Geodynamic controls on clastic-dominated base metal deposits

To meet the growing global demand for metal resources, new ore deposit discoveries are required. However, finding new high-grade deposits, particularly those not exposed at the Earth's surface, is very challenging. Therefore, understanding the geodynamic controls on the mineralizing processes can help identify new areas for exploration. Here we focus on clastic-dominated Zn–Pb deposits, the largest global resource of zinc and lead, which formed in sedimentary basins of extensional systems. Using numerical modelling of lithospheric extension coupled with surface erosion and sedimentation, we determine the geodynamic conditions required to generate the rare spatiotemporal window where potential metal source rocks, transport pathways, and host sequences are present. We show that the largest potential metal endowment can be expected in narrow asymmetric rifts, where the mineralization window spans about 1–3 Myr in the upper ∼ 4 km of the sedimentary infill close to shore. The narrow asymmetric rift type is characterized by rift migration, a process that successively generates hyper-extended crust through sequential faulting, resulting in one wide and one narrow conjugate margin. Rift migration also leads to (1) a sufficient life span of the migration-side border fault to accommodate a thick submarine package of sediments, including coarse (permeable) continental sediments that can act as source rock; (2) rising asthenosphere beneath the thinned lithosphere and crust, resulting in elevated temperatures in these overlying sediments that are favourable for leaching metals from the source rock; (3) the deposition of organic-rich sediments that form the host rock at shallower burial depths and lower temperatures; and (4) the generation of smaller faults that cut the major basin created by the border fault and provide additional pathways for focused fluid flow from source to host rock. Wide rifts with rift migration can have similarly favourable configurations, but these occur less frequently and less potential source rock is produced, thereby limiting potential metal endowment. In simulations of narrow symmetric rifts, the conditions to form ore deposits are rarely fulfilled. Based on these insights, exploration programmes should prioritize the narrow margins formed in asymmetric rift systems, in particular regions within several tens of kilometres from the paleo-shoreline, where we predict the highest-value deposits to have formed.

1 Introduction

Decarbonization, population increase, urbanization, and a growing middle class mean that current metal reserves will not meet projected future demands (IRP, 2020). Improved recycling will only go partway to meeting this increased demand (e.g. Elshkaki et al., 2018; IEA, 2021) and new ore deposit discoveries will be needed for metals such as copper, zinc, nickel, and lead (Lawley et al., 2022; Ali et al., 2017, and references therein). Approximately 50 % of global zinc and lead resources are found in sediment-hosted Zn–Pb deposits (Goodfellow et al., 1993; Singer, 1995; Mudd et al., 2017). The highest-value mineral systems within this broad family of deposits are clastic-dominated (CD-type), which are hosted by mixed siliciclastic–carbonate successions in Proterozoic and Paleozoic marine sedimentary basins (Wilkinson, 2014) formed in continental rifts, passive margins, and back-arcs (Leach et al., 2005, 2010). Large CD-type deposits can be found in the Mt. Isa Superbasin (Australia), the Kuna and Selwyn basins (USA and Canada), and other ancient rifted margins in Russia and China (Hoggard et al., 2020). However, there have been very few major CD-type discoveries in the last 25 years (Mudd et al., 2017), and developing exploration targets for CD-type deposits in ancient sedimentary basins 10² to 10⁵ km² in size remains a considerable challenge. This challenge is compounded by the exhaustion of resources located close to Earth's surface (< 100 m; Schodde, 2020, 2014), meaning future discoveries will be located at greater depths. Exploration under several hundred metres of cover will require new exploration methods on all spatial scales, from geodynamic to deposit scale, with an increasing contribution of geophysical, geochemical, and conceptional or geological methods (Schodde, 2020).

Mineral system models describe the key processes (often across different scales) involved in the genetic evolution of ore deposits (Hagemann et al., 2016; McCuaig et al., 2010). In sedimentary basins, key components of the mineral system model are the metal source, flow conduits for metal-bearing fluids, and a trap for concentrating metals at the deposit site (e.g. Wilkinson, 2014; Lawley et al., 2022). In CD-type deposits, the primary ore fluid was most likely derived from evaporated seawater, which when heated and circulated at depth was capable of leaching metals from permeable source (aquifer) units (Leach et al., 2005; Banks et al., 2002) during fluid–rock interaction with syn-rift siliciclastic source units or possibly with basement rocks or mafic volcanics (e.g. Wilkinson, 2023; Ayuso et al., 2004). Compartmentalization of the aquifer units by overlying sag-phase sediments could have resulted in locally overpressured aquifer units that were then periodically breached during periods of extension (e.g. Rodríguez et al., 2021). Syn-extensional faults also provide the main flow conduits for fluids between aquifer units and the site of metal deposition (e.g. Walsh et al., 2018; Hayward et al., 2021). The deposition of sulfide mineral CD-type ores is generally concentrated in units undergoing early to burial stages of diagenesis where there is a strong redox gradient, controlled either by a change in stratigraphic unit or by fluid mixing (Wilkinson, 2014). Sedimentary facies that are enriched in organic matter (e.g. mudstones) form particularly effective host rock units in which the required reduced sulfur may be generated via a number of biotic and abiotic processes (e.g. Kelley et al., 2004; Magnall et al., 2018).

Many of the world-class CD-type districts are known to be characterized by long periods of sedimentary basin evolution (> 100 Myr; e.g. Betts et al., 2006; Beranek, 2017; Gibson and Edwards, 2020). Nevertheless, the underlying geodynamic controls on transient ore-forming events remain unclear. Numerical modelling of hydrothermal fluid flow has focused on the Mt. Isa Superbasin (Yang et al., 2004a, b; Oliver et al., 2006; Sheldon et al., 2019, 2023) and Selwyn Basin (Schardt et al., 2008; Rodríguez et al., 2021). These basin-scale (< 70 km) hydrothermal studies demonstrate the interplay between the permeability and geometry of the hydrostratigraphic units, heating events, and opening of faults during tectonic events. However, all drivers and facilitators of fluid flow in these simulations greatly rely on assumptions and simplifications of basin, fault, and sedimentary layer geometry and other constraints such as thermal architecture. Moreover, deformation of the solid rock medium is often ignored. These conditions are, however, greatly influenced by the dynamics of rifting and the resulting margin architecture and sedimentary record.

In the field of geodynamics, numerical models have greatly advanced our understanding of sedimentary rift basin formation. The lithospheric strength distribution – dependent on composition, temperature, pressure, strain rate and/or extension velocity, and weakening mechanisms – controls deformation patterns and the resulting margin architecture (Pérez-Gussinyé et al., 2023; Brune et al., 2023). For example, coupling of upper crust and mantle lithosphere determines whether deformation occurs over a broad area (for a weak lower crust) or narrow area (strong lower crust; Buck, 1991). Strain weakening breaks the symmetry of deformation, leading to asymmetric extension (Huismans and Beaumont, 2003). High temperatures and viscous strain softening at the tip of the dominant asymmetric fault and cooling in its footwall generate a strength gradient that sustains rift migration (Brune et al., 2014). New faults continuously form in the direction of migration, accreting material of the opposite margin to the wide margin left behind. This produces highly asymmetric margins, with the wide margin characterized by hyperextension. Through varying extension velocity, crustal layer thickness, strain weakening, and other parameters affecting the lithospheric strength distribution, a range of rifting styles and conjugate margin geometries has been observed in numerical simulations (e.g. Dyksterhuis et al., 2007; Gueydan et al., 2008; Sharples et al., 2015; Tetreault and Buiter, 2018), in particular a metamorphic core complex mode, wide rifting, and narrow rifting leading to symmetric or asymmetric margins.

The more recent coupling of geodynamic modelling software with surface processes codes (Andrés-Martínez et al., 2019; Wu et al., 2019; Beucher and Huismans, 2020; Balázs et al., 2021; Wolf et al., 2021; Neuharth et al., 2022a) has also helped unravel the interaction of sediment (un)loading and rift tectonics. For one, sediment erosion and deposition increase the longevity and offset of faults in the proximal margin, prolong the different phases of rifting, and therefore delay the transition to oceanic spreading (Olive et al., 2014; Theunissen and Huismans, 2019; Neuharth et al., 2022a). Sediment loading can also enhance ductile crustal flow in weak crust, promoting sag basin formation (Theunissen and Huismans, 2019). Sediment blanketing maintains higher temperatures in the thinned crust, favouring ductile crustal deformation in the distal margins during late rifting (Andrés-Martínez et al., 2019). Beucher and Huismans (2020) found, however, that for a weak lithosphere, distributed deformation reduces the accommodation space and the spread-out sediment load will have little feedback on tectonics. Understanding how different geodynamic and erosional parameters might influence the spatial and temporal alignment of CD-type mineral system components thus has the potential to significantly improve basin prospectivity mapping (e.g. Hoggard et al., 2020; Lawley et al., 2022).

In this contribution, we provide the first numerical modelling study employing a geodynamic model coupled to a surface processes model to understand ore formation in sedimentary rift basins. Having identified the main components (i.e. source, fault fluid conduit, trap and host) that characterize the CD-type mineral system, we simulate the spatial and temporal evolution of these components in order to investigate the geodynamic controls on prospectivity. We focus on the role of rift type (narrow vs. wide, asymmetric vs. symmetric), erosion and sedimentation efficiency, and plate velocities in the generation of basins with potential source and host rocks connected by faults. Importantly, we identify three potential ore-deposit-forming mechanisms and find that asymmetric narrow rifting favours all of these mechanisms. We discuss these results in light of known major CD-type deposits as well as their implications for mineral exploration.

2 Methods

To investigate the geodynamic controls on CD-type deposit formation, we use a coupled numerical modelling approach that combines the dynamics of continental rifting (Sect. 2.1) with surface processes (Sect. 2.2). In Sect. 2.3 we describe how the surface processes information is applied to track the sedimentary infill of the tectonically created basins. The sediment types combined with the thermal and mechanical structure of the basins are used to identify ore formation mechanisms in a post-processing step described in Sect. 2.4. Finally, we list the model assumptions (Sect. 2.5) and define the model setups (Sect. 2.6) before presenting the modelling results in Sect. 3.

2.1 Governing equations for continental rifting

Using the ASPECT software (Kronbichler et al., 2012; Heister et al., 2017; Rose et al., 2017; Glerum et al., 2018; Bangerth et al., 2021a, b; Glerum, 2023), we solve for the conservation equations of momentum (1), mass (2), and energy (3) using the extended Boussinesq approximation.

$$\begin{array}{}\text{(1)}& {\displaystyle }-\mathrm{\nabla }\cdot \left[\mathrm{2}\mathit{\eta }\dot{\mathit{ϵ}}\right(\mathit{u}\left)\right]+\mathrm{\nabla }P=\mathit{\rho }\mathit{g}\text{(2)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0}\text{(3)}& {\displaystyle }\begin{array}{rl}\stackrel{\mathrm{‾}}{\mathit{\rho }}{C}_P({\displaystyle \frac{\partial T}{\partial t}}& +\mathit{u}\cdot \mathrm{\nabla }T)-\mathrm{\nabla }\cdot k\mathrm{\nabla }T=\stackrel{\mathrm{‾}}{\mathit{\rho }}H\\ & +\mathrm{2}\mathit{\eta }\dot{\mathit{ϵ}}\left(\mathit{u}\right):\dot{\mathit{ϵ}}\left(\mathit{u}\right)+\mathit{\alpha }\stackrel{\mathrm{‾}}{\mathit{\rho }}T(\mathit{u}\cdot \mathit{g})\end{array}\end{array}$$

Here η is the effective viscosity, $\dot{\mathit{ϵ}}\left(\mathit{u}\right)$ the strain rate depending on velocity u,P the total pressure, g the gravity vector, $\mathit{\rho }={\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{1}-\mathit{\alpha }\left(T-T_{\mathrm{ref}}\right)\right)$, T_ref the adiabatic reference temperature, T the temperature, $\stackrel{\mathrm{‾}}{\mathit{\rho }}$ the adiabatic reference density, C_P the specific isobaric heat capacity, t the time, k the thermal conductivity, α the thermal expansivity, and H the intrinsic specific heat production. ρ₀, C_P, k, H, and α are all constant per compositional material.

To track such materials and other fields (e.g. upper crust and plastic strain), we solve an advection equation for each field c_i:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial c_i}{\partial t}}+\mathit{u}\cdot \mathrm{\nabla }{c}_i=q_i\left(\dot{\mathit{ϵ}}\left(\mathit{u}\right)\right).\end{array}$$

Here the source term q_i is used to update the plastic and viscous strain as $\mathrm{\Delta }t\dot{\mathit{ϵ}}$ in the plastic or viscous domain, respectively, with Δt the time step size.

The effective viscosity in Eqs. (1) and (3) is equal to a composite (η_comp) of diffusion (diff) and dislocation (disl) creep flow laws when not on yield (Karato and Wu, 1993; van den Berg et al., 1993; Karato, 2008) and equals the plastic viscosity η_p based on a Drucker–Prager yield criterion when on yield (Willett, 1992; Davis and Selvadurai, 2002; Glerum et al., 2018):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathit{\eta }}_{\mathrm{diff}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}_{\mathrm{w}\mathrm{diff}}^{-\mathrm{1}}{d}^m\exp \left({\displaystyle \frac{E_{\mathrm{diff}}+P{V}_{\mathrm{diff}}}{RT}}\right),\\ {\mathit{\eta }}_{\mathrm{disl}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}_{\mathrm{w}\mathrm{disl}}^{\frac{-\mathrm{1}}{n}}{\mathit{ϵ}}_{\mathrm{eff}}^{\frac{\mathrm{1}-n}{n}}\exp \left({\displaystyle \frac{E_{\mathrm{disl}}+P{V}_{\mathrm{disl}}}{nRT}}\right),\\ {\mathit{\eta }}_{\mathrm{comp}}& ={\displaystyle \frac{{\mathit{\eta }}_{\mathrm{diff}}{\mathit{\eta }}_{\mathrm{disl}}}{{\mathit{\eta }}_{\mathrm{diff}}+{\mathit{\eta }}_{\mathrm{disl}}}},\\ {\mathit{\eta }}_p& ={\displaystyle \frac{C\cos \left({\mathit{\varphi }}_w\right)+P\sin \left({\mathit{\varphi }}_w\right)}{\mathrm{2}{\mathit{ϵ}}_{\mathrm{eff}}}},\\ {\mathit{\varphi }}_w& =F_{\mathit{\varphi }}\left({\displaystyle \frac{\max \left({\mathit{ϵ}}_{\min },\min \left({\mathit{ϵ}}_p,{\mathit{ϵ}}_{\max }\right)\right)-{\mathit{ϵ}}_{\min }}{{\mathit{ϵ}}_{\max }-{\mathit{ϵ}}_{\min }}}\right)\mathit{\varphi }\\ & \left(\mathrm{same}\mathrm{for}{A}_{\mathrm{w}}\right).\end{array}\end{array}$$

The internal friction angle ϕ_w and the viscous prefactor A_w of each material are linearly reduced (weakened) up to a factor F of 0.25 or 0.75 over the plastic and viscous strain intervals of ϵ_min=0 to ϵ_max=2 (e.g. Huismans and Beaumont, 2003; Naliboff and Buiter, 2015; Le Pourhiet et al., 2017). This strain weakening represents frictional softening and viscous strain softening, respectively, two mechanical weakening mechanisms observed in laboratory experiments and in nature (Bürgmann and Dresen, 2008). The former leads to fault localization and the latter to localization in ductile shear zones, mimicking the drop in viscosity upon grain size reduction by large strains (White et al., 1980; Handy, 1989; Rybacki and Dresen, 2004; Austin and Evans, 2009). E is the activation energy, V the activation volume, R the gas constant, C the cohesion, and ${\mathit{ϵ}}_{\mathrm{eff}}=\sqrt{\left|{\dot{\mathit{ϵ}}}_{\mathrm{II}}\right|}$ the effective strain rate (scalar), with ${\dot{\mathit{ϵ}}}_{\mathrm{II}}$ the second moment invariant of the deviatoric strain rate tensor. All parameter values can be found in Table 1.

Table 1 Material and other model parameters and their symbols, values, and units. See Yuan et al. (2019a, b), Guerit et al. (2019), Neuharth et al. (2022a), and references therein for FastScape parameter values.

Figure 1 Schematic illustrations of the sedimentary processes in our coupled ASPECT–FastScape workflow and the parameters used in Eqs. (6)–(8). (a) Endmember case of purely continental sediment erosion, transport, and deposition in the continental and marine domain. (b) Endmember case of purely marine sedimentation. (c) Combined continental and marine erosion, transport, and deposition, including reworking, as used in this work. (d) Definitions of continental source and marine host rock units. Abbreviations: DD – deposition depth, MF – fraction of marine sediments, SF – silt fraction.

2.2 Governing equations for surface processes

As described in Neuharth et al. (2022a, b), two-way coupling with the surface process code FastScape (Braun and Willett, 2013; Yuan et al., 2019a, b) handles the erosion, transport, and deposition of sediment in both the marine and continental parts of the ASPECT model domain (see Figs. 1 and 2). These processes are governed by the stream power law with a sediment transport and deposition term, hillslope diffusion, and marine transport and deposition. For topography h above sea level (SL) height h_SL, the following equation is thus used (Yuan et al., 2019a):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle \frac{\partial h}{\partial t}}& =U-K_f{A}^m{S}^n+{\displaystyle \frac{G}{A}}\underset{A}{\int }\left(U-{\displaystyle \frac{\partial h}{\partial t}}\right)\mathrm{d}A+K_{\mathrm{d}}{\mathrm{\nabla }}^{\mathrm{2}}h\\ & +\mathit{v}\cdot \mathrm{\nabla }h.\end{array}\end{array}$$

Here U is the tectonic uplift velocity, K_f the incision rate parameter (for the stream power law, SPL), A the drainage area (SPL), m the drainage area exponent (SPL), S the slope (SPL), n the slope exponent (SPL), G the dimensionless deposition and transport coefficient (enriched SPL), K_d the transport coefficient (hillslope diffusion), and v the horizontal velocities used to advect the topography. U and v are supplied by ASPECT at each time step; all values used for the other parameters can be found in Table 1.

In the marine domain, where h is smaller than h_SL, the following equation is solved (Yuan et al., 2019b):

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial h}{\partial t}}=K_{\mathrm{M}}{\mathrm{\nabla }}^{\mathrm{2}}h+\mathit{v}\cdot \mathrm{\nabla }h+Q_s+R_{\mathrm{M}}.\end{array}$$

Here K_M is the marine sediment transport coefficient, and $Q_{\mathrm{s}}={\int }_A\dot{e}\mathrm{d}A/\left(\mathrm{\Delta }x\mathrm{\Delta }y\right)$ is the sediment flux coming from the continental domain at points along the shoreline (Q_s is nonzero only along the shoreline), with the erosion rate $\dot{e}$ and the surface area of the whole domain ΔxΔy. R_M is the user-set background marine sedimentation rate that deposits a homogenously thick layer of sediments in the marine domain.

2.3 Formation of basin stratigraphy

The ASPECT side of the coupling to FastScape tracks the deposited sediment and its age and deposition depth (DD) on compositional fields that are advected according to Eq. (4) (Neuharth et al., 2022a). We have extended the interface to also track (1) the dominant grain size of the continental sediment (coarse or fine, using FastScape's “silt fraction”, hereafter called SF) within the marine domain and (2) the fraction of marine sediments (MF) of the sum of marine and continental sediments deposited within one ASPECT time step at a particular location along the submarine surface. On land, FastScape's silt fraction SF is 0 everywhere; upon entering the marine domain sediments receive the default value of 0.5 (Fig. 1a). Due to the different transport properties of silt and sand within the marine domain (K_M, see Table 1), this fraction then evolves with transport distance from the shoreline (staying between 0 and 1). Marine sediments are deposited everywhere in the sub-sea-level domain at a user-specified rate R_M by increasing the ASPECT surface height in the marine domain before passing it to FastScape (Fig. 1b). The increase is equal to the time step size times the deposition rate, but the new height cannot exceed the sea level height. The fraction of marine sediments (MF) out of the total thickness of sediments deposited in one location in one ASPECT time step is computed as

$$\begin{array}{}\text{(8)}& \mathrm{MF}={\displaystyle \frac{h_{\mathrm{marine}}-h_{\mathrm{0}}}{h_{\mathrm{total}}-h_{\mathrm{0}}-U\cdot \mathrm{\Delta }t}}.\end{array}$$

Here h₀ is the height above the initial ASPECT surface at the beginning of the ASPECT time step, h_marine the height after marine sediments have been deposited, h_total the height that has been subsequently returned by FastScape (including uplift, erosion and sedimentation effects), U the uplift velocity provided to FastScape, and Δt the ASPECT time step size. If material has been eroded (i.e. $h_{\mathrm{total}}-U\cdot \mathrm{\Delta }t-h_{\mathrm{0}}<\mathrm{0})$, MF is set to zero.

2.4 Proxies for mineral system components and metal endowment

The general model for CD-type metallogenesis describes the leaching of metals from a source rock by fluids, transport of these fluids upwards along faults, and deposition of the metals they carry in units undergoing early to burial stages of diagenesis in marine sedimentary basins (e.g. Wilkinson, 2003; Heinrich and Candela, 2014; Magnall et al., 2021). As base metal solubilities increase with temperature and salinity (Yardley, 2005), evaporated seawater that has been heated up has the potential to leach and transport high concentrations of metals (e.g. Banks et al., 2002; Wilkinson et al., 2009). Leaching occurs at depth through fluid–rock interaction with the intra-basinal syn-rift siliciclastic rocks such as sandstones (Ayuso et al., 2004; Leach et al., 2005). The bottom sedimentary sequences and the fluids flowing through them can be heated up by temperature anomalies generated by rifting (e.g. through proximity to hot asthenosphere after thinning of the lithosphere), driving thermal fluid convection and allowing for ore-forming fluid temperatures between 150 and 300° (e.g. Cooke et al., 2000; Magnall et al., 2016), which can result in Zn and Pb solubilities of hundreds to thousands of parts per million (ppm) (e.g. Wilkinson et al., 2009). Fault activity (e.g. seismic pumping and pulsing as well as fault valving) and passive upward flow then bring the fluids to the surface (Walsh et al., 2018), mechanisms that focus the fluids and concentrate the metals in the host rock. Host rocks to the mineralization are typically carbonaceous marine mudstones and siltstones with some carbonate units (e.g. Ridley, 2013; Wilkinson, 2014). In this syn-diagenetic environment, base metal solubilities decrease along steep thermal and chemical gradients. When the fluids encounter reduced sulfur, sulfide minerals are formed, which are only stable under reducing conditions (Cooke et al., 2000; Ridley, 2013). Sedimentary facies that are enriched in organic matter (e.g. mudstones) therefore form particularly effective host rock units as they can produce large volumes of reduced sulfur via a number of biotic and abiotic pathways (e.g. Kelley et al., 2004; Magnall et al., 2018, 2020). For bacteriogenic (biotic) processes to occur, maximum temperatures of the host rock should be below ∼ 110° (Takai et al., 2008; Ridley, 2013). At temperatures over 110°, abiotic pathways of sulfate reduction prevail, e.g. thermochemical sulfate reduction (Machel, 2001).

We translate the above description of the CD-type mineral system into a series of simplified conditions that identify the main components of the system in our simulations using the predicted dominant sediment type, its deposition depth, and the temperature field. These conditions are evaluated in a post-processing step, in which the ASPECT results are analysed using ParaView (Ahrens et al., 2005; Ayachit, 2015) Python scripts as follows.

1. Based on the tracked fields described in Sect. 2.3, we first separate sediments from basement materials (e.g. upper and lower crust). We then further divide and colour-code the sediments into five types: (1) sandstone (where silt fraction SF is smaller than 0.5 and marine fraction MF is smaller than 0.5; yellow colours in the figures), (2) siltstone (where SF > = 0.5 and MF < 0.5; orange), (3) limestone (where MF > = 0.5 and deposition depth DD below sea level < 200 m; light purple), (4) carbonates (where MF > = 0.5 and 200 m < = DD < 4500 m; medium purple), and (5) silicious marine sediments (MF > = 0.5 and DD > = 4500 m, which represents the carbonate compensation depth (e.g. Middelburg, 2019; dark purple). Note that sediment types are classified according to their dominant component (Fig. 1c). For example, when sediment is labelled as sandstone, it means that the dominant (i.e. more than 50 %) grain size is coarse and from a continental source region, not that no continental fine-grained or marine sediment is present. When sediments in the marine domain are eroded after being deposited, their previous type is lost and they are classified as continental in the computation of MF even if they were labelled as marine before.

2. Combining these sediment types and the temperature field, we identify and colour-code potential source and host rock (Fig. 1d). Where temperatures in the sandstone exceed 250 °C, it is coloured dark green, indicating source rock conditions favourable for enhanced metal leaching and the generation of ore-forming solutions. Where marine sediments (types 3–5) are below a temperature of 150 °C, they are coloured light green, indicating potential host rock, i.e. where redox becomes the predominant factor controlling metal solubilities.

3. Potential deposit-forming mechanisms are then defined as events where an active fault (high strain rate; approximately $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{15}}$ s⁻¹ and above, see greyscale in figures) or an inactive fault (high plastic strain of ∼ 3 and above, low strain rate) connects source and host rock within the same or neighbouring basins. Ore-forming fluids will precipitate metals when they ascend along a fault and first encounter prospective host rock units (the thickness of the host rock is thus less important than the area of source rock). We assume that active faults are always more permeable than the surrounding rock and therefore act as preferred upward pathways for fluids. However, as summarized by Walsh et al. (2018), faults also facilitate replenishing downward flow of fluids, and their permeability can be greatly decreased by, for example, the generation of clay smears. Shearing of unlithified sediments with high clay content (like CD-type host rock) can lead to lower permeabilities than in the surrounding shallow (< 3 km) sediments (e.g. Barnicoat et al., 2009; Neuzil, 2019; Sheldon et al., 2023). For a discussion of fault permeability, see also Bense et al. (2013) and Gleeson and Ingebritsen (2016).

Within ASPECT we also compute at each time step the total area of potential source and host rock using the same definitions as above and integrating over each mesh element. Based on the total source rock area A, we estimate the possible metal endowment E as

$$\begin{array}{}\text{(9)}& E=A\cdot \mathrm{d}z\cdot \mathit{\rho }\cdot F_{\mathrm{leaching}}\cdot C_{\mathrm{initial}},\end{array}$$

where ρ is the rock density, dz=10 km is the assumed in-plane extent of the basins to create a source rock volume (based on basin extent in Blaikie and Kunzmann, 2020, and on the assumptions of Manning and Emsbo, 2018), F_leaching is the leaching capacity factor assumed to be 0.65 (from Windermere Supergroup leaching experiments; Lydon, 2015), and C_initial is the initial metal concentration of the source rock, taken to be 100 ppm for Zn and 23 ppm for Pb (Lydon, 2015).

Figure 2 Numerical model setup. (a) ASPECT and FastScape model domains for reference model NA-2, showing initial composition (beige and grey areas), initial temperature (white isotherms), initial plastic strain (blue), boundary conditions, and initial FastScape topography with respect to sea level at −200 m from the initial ASPECT model surface at y= 300 km. Rifting driven by the prescribed tensional horizontal velocities is initially localized in the centre of the domain due to the higher-amplitude plastic strain and thicker upper crust. The bottom boundary of the ASPECT domain is open to flow, and its normal traction is prescribed the initial lithostatic pressure at x= 10 km and y= 0 km. FastScape is initialized with a random topography of ±5 m. (b) Temperature and strength profiles in unperturbed and perturbed (centre of domain) locations for each of the three modelled rift types: narrow asymmetric (NA), narrow symmetric, and wide rifting assuming an effective strain rate of 10⁻¹⁵ s⁻¹. (c) Strain rate patterns (black) in the crust at the horizontal centre of the domain after 1.5 Myr for three out of the nine different initial plastic strain configurations of the narrow asymmetric rift type. The effective strain rate is defined as $\sqrt{\left|{\dot{\mathit{ϵ}}}_{\mathrm{II}}\right|}$, the square root of the second moment invariant of the deviatoric strain rate.

2.5 Model assumptions

The coupled geodynamic–surface process model cannot explicitly resolve hydrothermal fluid flow and fluid–rock interactions; therefore, the proxies for metal endowment provide maximum estimates under the simplifying assumptions that (1) the source rocks are permeable enough to support metal transport by fluid flow, (2) the hydrology provides an efficient focusing mechanism along fault structures, (3) the thermal structure remains close to the conduction-dominated profile calculated by the geodynamic model during most of the ore-forming event, and (4) the availability of reduced sulfur from the host rock is sufficient to provide an active metal trap.

Although our workflow does not directly account for the above sub-fault-scale processes, it nevertheless captures the essence of CD-type deposit formation by identifying source and host rocks and fluid migration pathways. As such, this workflow uses the same solution to the cross-scale problem as petroleum and mineral system analysis (Allen and Allen, 2013, Chap. 11; McCuaig et al., 2010). Similar workflows to deduce the first-order impact of fluid flow and associated geochemical alteration from numerical forward models have been shown to be successful avenues for estimating hydrogen generation at mid-ocean ridges (Rüpke and Hasenclever, 2017) and rifted margins (Liu et al., 2023), as well as carbon dioxide release at extensional plate boundaries (Hasenclever et al., 2017; Gorczyk and Gonzalez, 2019). To truly cross the scales from large-scale dynamics to mineralization, the inclusion of hydrothermal flow and reactive transport modelling (e.g. Yapparova et al., 2019) into the combined modelling of geodynamic and surface processes as presented in this paper is required.

While our simulations account for self-consistent interactions between thermal, mechanical, and sedimentary processes, they also involve several geodynamic limitations. For one, the governing equations approximate the Earth as a highly viscous fluid on geological timescales, modifying the viscosity to allow for non-viscous behaviour such as brittle deformation. Second, our simulations are 2D and therefore ignore potential along-strike effects on the structural heterogeneity of the rift and margins (e.g. Naliboff et al., 2020; Gouiza and Naliboff, 2021), effects of oblique extension (Brune et al., 2012, 2018), and, on a smaller scale, the possible compartmentalization of subbasins into a complex, multi-generational, multi-directional 3D fault geometry as observed in the McArthur Basin (Australia) that hosts several world-class deposits (Blaikie and Kunzmann, 2020). Furthermore, as the initial Proterozoic state of the lithosphere is mostly unknown, we start off our simulations with a relatively smooth initial lithosphere geometry and approximate all possible sources of tectonic inheritance with initial plastic strain (see Sect. 2.6). Finally, we do not include melting and magma flow and thus also ignore their effect on heat flow, faulting, the overall stress state, the generation of oceanic crust, and its potential role as metal source rock (e.g. Murphy et al., 2011).

2.6 Model setup

The model setup represents a 2D extensional lithosphere–asthenosphere system (for ASPECT input files and source code, see Glerum, 2023). The model domain spans 700×300 km and is discretized with a resolution that varies from 5 km at the bottom of the domain to 312.5 m in the upper 42 km of the domain. The reference case includes a 20 km thick upper crust, a 15 km thick lower crust, and 85 km of mantle lithosphere; see Fig. 2. The rheologies assigned to these materials are wet quartzite (Rutter and Brodie, 2004), wet anorthite (Rybacki et al., 2006), dry olivine (Hirth and Kohlstedt, 2004), and wet olivine for the asthenosphere (Hirth and Kohlstedt, 2004), respectively. See Table 1 for the properties of each material layer. The upper crust and the mantle lithosphere are locally thickened or thinned around x= 350 km using a Gaussian shape with a standard deviation of 60 km. To further localize deformation, the plastic strain field is initialized with random noise with a maximum amplitude of 0.25 centred at x= 350 km (see Fig. 2a).

The initial temperature follows a 1D continental steady-state geotherm (Chapman, 1986) at each x coordinate, constrained by a lithosphere–asthenosphere boundary (LAB) temperature (coinciding with the bottom of the mantle lithosphere field at the start of the model run) of 1350 °C and a surface temperature of 20 °C (see Table 1 for other relevant parameters). Below the lithosphere, the initial temperature follows an adiabat with a surface potential temperature of 1303 °C such that at the LAB the geotherm and adiabat temperatures match (at 1350 °C). As such, we assume that the tectonic settings of the Paleoproterozoic in which most CD-type deposits formed can be represented by a setup representative of present-day continental rifting.

Extension is prescribed through outward horizontal velocities on the left and right vertical boundaries (see Fig. 2a). The reference total prescribed horizontal velocity is 10 mm yr⁻¹, which resembles the velocity of the initial slow rift phase (Brune et al., 2016) and falls in the ultra-slow range that produces magma-poor margins (Pérez-Gussinyé et al., 2023). On the bottom boundary, the initial lithostatic pressure below reference lithosphere at 300 km depth (at x= 1 km) is prescribed as the normal traction component such that material is free to enter or leave the domain in response to the internal stress state. Temperature and composition are fixed to their initial values at this boundary (i.e. the inflowing material is asthenosphere). The surface boundary is deformed using FastScape (see Sect. 2.2), with full stabilization according to Kaus et al. (2010) and Rose et al. (2017). The FastScape sea level is set to 200 m below the initial surface of the model domain, so at y=299.8 km, where y is the vertical upward coordinate. Along the deforming surface boundary, the fixed boundary condition for the deposition depth (DD) field is set to DD $=\mathrm{299.8}-y$ (in km), the sediment age field to time in millions of years, and the sediment field to 1. Other material fields are kept fixed at a value of zero along the surface boundary, while SF and MF are set to the values computed as described in Sect. 2.3.

The simulations presented in the following section differ in the strain weakening factor F (see Eq. 5), the thickness of the lithospheric layers, the initial perturbation of these layers to focus deformation, surface process parameters, and the initial plastic strain. By varying the first three parameters, we obtain three endmember rifting styles: narrow asymmetric, narrow symmetric, and wide rifting (e.g. Huismans and Beaumont, 2003; Brune et al., 2014; Svartman Dias et al., 2015). For the narrow asymmetric cases, the upper crust is locally thickened to 25 km and the mantle lithosphere locally thinned by 15 km at the domain's centre (see Fig. 2b). For narrow symmetric rifting, strain weakening is reduced by increasing the weakening factor F from 0.25 to 0.75 with respect to the narrow asymmetric cases. Wide rifting is obtained for an upper crust of 25 km, a lower crust of 20 km, and a mantle lithosphere layer of 75 km of which the upper crust is locally thickened to 30 km and the mantle lithosphere thinned to 70 km. The weakening factor F is again 0.25. For the narrow asymmetric simulations, we subsequently vary the surface process parameters R_M, K_M, and K_f and the extension velocity.

For each of the above cases, nine simulations are run with different initial plastic strain configurations. The initial strain configuration determines the early fault pattern and coalescence, as shown for three simulations in Fig. 2c, and therefore the subsequent co-evolution of thinning and faulting. The initial strain represents possible tectonic inheritance, such as foliation, lattice-preferred-orientation-induced mechanical anisotropy (Tommasi and Vauchez, 2001), and other pre-existing structures where continental rifts prefer to localize. The nine initial configurations account for possible along-strike variability in inherited heterogeneity (Richter et al., 2021) and allow us to assess the robustness of the modelling results. To this end, we both discuss results of individual simulations and report averages of, for example, the number of basins favourable to ore formation and area of source rock.

Each simulation runs for 25 Myr of model time, which is enough to reach continental breakup, after which the syn-rift, syngenetic, and syn-diagenetic CD-type deposits are no longer expected to form. Although studies have suggested a post-rift or syn-inversion timing for some CD-type deposits (e.g. Gibson and Edwards, 2020), syn-rift deposits have been most comprehensively studied and are the focus of our study.

Figure 3 Model evolution of narrow asymmetric rift simulation NA-4. (a–f) Time series of the whole domain showing different materials, isotherms (white and black), strain rate (greyscale), sediment type (orange, yellow, and purple shades), and velocity (brown vectors). Black rectangles indicate zoomed-in areas in (g)–(l). (g–l) Time series of zoom-ins focusing on FastScape topography (top) and sedimentary basins (bottom), additionally showing plastic strain (blue scale) and sediment age (brown contours in steps of 5 Myr). Coarse continental sediments (sandstone) that exhibit temperatures higher than 250 °C are classified as source rock and coloured dark green. Marine sediments at temperatures below 150 °C are potential host rocks for metal deposition and are coloured in light green shades. The black rectangle in (k) indicates the location of a further zoom-in in Fig. 4.

3 Results

We first present the overall rift evolution and describe the favourable conditions for ore formation in the simulated sedimentary basins for each rifting type (narrow asymmetric, narrow symmetric, and wide rifting). Then we investigate the effect of varying surface process parameters and extension velocity on source and host rock creation in narrow asymmetric rifts.

3.1 Narrow asymmetric rifts

Out of the nine simulations of narrow asymmetric rifting (simulations NA-1 to NA-9), the representative temporal evolution of simulation NA-4 is presented in Fig. 3. After an initial phase of distributed deformation (Fig. 3a and g), by 2.5 Myr two conjugate normal faults have formed. By 5 Myr (Fig. 3b and h), the currently active left border fault (with black indicating high strain rates) extends to the base of the lithosphere and new faults are forming in the rift centre. Rift shoulder topography (maximally 2139 m) has reached its highest point and sedimentary basins are being filled with mostly coarse continental sediments (sandstone; yellow areas). At 10 Myr, the rift centre has started migrating towards the right, and the old left shear zone is abandoned (Fig. 3c and i; blue indicates high plastic strain). An exhumation channel runs from the uprising asthenosphere to the active faults, along which lower crust flows up from the narrow margin. This counterbalances the thinning of the crust by the new faults, postponing continental breakup (e.g. Brune et al., 2014). Marine sediments (shallow limestones) deposited in the rightmost basin (Fig. 3i; shown and labelled as light-green potential host rock) are subsequently buried by new continental sediments (Fig. 3j). Rift shoulder topography is eroded and the sequential faulting of rift migration continues. By 15 Myr, the rightmost basin has heated up enough for the lowermost sandstone layers to be classified as source rock (dark-green area in Fig. 3j). They are connected by the now inactive border fault to the marine sediments higher up in the stratigraphy, which are cool enough to act as host rock (light green). This constellation aligns potential source rock, fluid pathway, and host rock and could promote ore formation by hydrothermal fluid flow. As the border fault is inactive, the constellation is not the most favourable constellation encountered in the NA simulations, and we therefore term these conditions “ore formation mechanism 3”; see Fig. 5.

5 million years later (Fig. 3k), the border fault has been reactivated, but the marine sediments along the fault are no longer in the host temperature window (represented by the colour change from light green to light purple in Fig. 3j and k). However, the sequential faulting associated with rift migration now leads to smaller synthetic faults fracturing the major basin adjacent to the border fault, tapping the potential source rock, which has grown in size during basin evolution, and providing a pathway to the potential host rock lying above. We consider these conditions more favourable because a larger source rock region is now connected to an overlying host rock by active faults instead of an inactive fault, so we label them as “ore formation mechanism 1” (Fig. 4). At the same time, one of the deeper basins on the wide margin develops and maintains a small pocket of source rock (Fig. 3e–k) because the rising asthenosphere that has replaced the thinning lithosphere maintains higher temperatures in the extended crust above. This source rock is connected to overlying host rock by an inactive fault; this constitutes another instance of ore formation mechanism 3. Rift migration then ceases after 20 Myr and the right border fault becomes inactive. By 22 Myr, oceanic spreading has started and highly asymmetric margins have formed (Fig. 3f and l).

Figure 4 The most favourable ore formation mechanism observed in simulations of narrow asymmetric rifting, i.e. a scenario that most favourably combines source and host rock as well as fault activity as required for metallogenesis. The upper- and lowermost rows show all rift basins and the zoom-ins in between only the relevant subbasin over time. Note that the zoom-in location changes over time to stay centred on the basin. White arrows and annotations indicate the predicted depth of mineralization at the bottom of host rock sequences. Blue arrows and annotations indicate the minimum distance of the mineralization to the shoreline. Dotted white arrows delineate the interpreted fluid flow paths, with potential deposit locations indicated by black stars.

Source rock can still be found at the sea floor of the rightmost basin, in the spreading centre, and in one basin on the wider left margin, while host rock covers large parts of the wider margin. However, at this stage, no more potential ore-forming mechanisms can be identified.

The nine narrow asymmetric rifting simulations with varying initial plastic strain configurations evolve very similarly. All simulations produce asymmetric margins, with five simulations showing rift migration towards the left and four towards the right (i.e. the direction of migration is random). The quantity of sediments deposited increases approximately linear over time and totals on average 962 km² after 25 Myr. Only a fraction of those sediments could actually serve as source (≤ 21.9 km²; Fig. 6a) or host (≤ 271.7 km²) rock because it is either not the right type of sediment or not in the right temperature window. Whereas the amount of host rock increases steadily over time, source rock only starts to form after 8 Myr (in or close to the rift centre, Figs. 3i and 6a) when the rift shoulders have provided significant sediment supply to the marine domain. Source rock area massively increases from around 15 Myr to its maximum around 20 Myr (for seven out of the nine simulations) as source rock forms in the major basin on the narrow margin, before slowly declining during the last ≤ 5 Myr (with the exception of NA-7, Fig. 6a). The decline is temporally related to the stabilization of the rift centre (starting between 17.5 and 21.5 Myr except for NA-7) and its transition to oceanic spreading (starting from 20 Myr, grey area in Fig. 6a, except for NA-7).

The model results contain three resolvable parameters (see Sect. 2.4) that are critical for the overall mineral system: (1) source rock for fluids to leach metals from, (2) host rocks for fluids to deposit the metals upon reaching them, and (3) active and/or inactive faults that act as fluid pathways connecting source and host rock. Figures 4 and 5 zoom in on the three ore formation mechanisms combining these ingredients that occur in the suite of nine narrow asymmetric rift simulations. Ore formation mechanism 1 (Fig. 4, also Fig. 3k) shows a scenario in which the faults of the rift centre migrate into the major active basin on the narrow margin. This basin contains both source rock and host rock, which are subsequently connected by the new faults providing upward fluid pathways. Ore formation mechanism 2 (Fig. 5a) shows how the major basin on the left narrow margin develops both source rock and host rock that are connected by the reactivated (or continuously active) border fault. In ore formation mechanism 3 (Fig. 5b), the source rock forming in a basin is connected to host rock by an inactive fault. This can occur in both the wider margin and the major basin (Fig. 3k and j, resp.).

To summarize the occurrence of components required for ore formation, we count the number of basins in which (1) source rock occurs, (2) source rock and host rock co-occur, and (3) source rock and host rock are connected by an inactive (ore formation mechanism 3) or (4) an active fault (mechanisms 1 and 2). This analysis is performed at 2.5 Myr intervals. The maximum counts are averaged over the nine simulations with different initial strain configurations and plotted in Fig. 7. On average, narrow asymmetric rifts contain two basins that align with the most prolific ore formation mechanisms (1 and 2); i.e. they contain source and host rock connected by active faults. Note that there are also basins with host rock and active faults but no source rocks (e.g. Fig. 5a at 15 Myr) and basins with only host rock and inactive faults (e.g. border basin in Fig. 3i and in all wide margins). These cases are, however, not represented in Fig. 7, as they do not have a potential metal source and therefore cannot be a site of mineralization.

Figure 5 Ore formation mechanisms 2 and 3 in narrow asymmetric rift simulations. Note that the NA-9 snapshots have been mirrored such that rift migration is directed towards the right instead of the left, as for ore formation mechanism 1 in Fig. 4. See Fig. 4 for more information.

3.2 Narrow symmetric rifts

A representative narrow symmetric simulation is shown in Fig. 8 (NS-8 out of NS-1 to NS-9). Initially, the symmetric rift evolution resembles the asymmetric one. However, the phase of distributed faulting is prolonged compared to narrow asymmetric rifting (compare Figs. 8b and 3b) and a symmetric sedimentary infill of sand, silt, and shallow marine sediments is deposited in the area of deformation. Similarly high rift shoulders develop, with a maximum topography of 2070 m. By 10 Myr, two major faults dipping away from the rift centre accommodate most of the deformation and the lithospheric mantle has broken up. Along both large shear zones, a small pocket of source rock forms that is connected to host rock deposited at great water depth through antithetic faults (ore formation mechanism 1). By 15 Myr, continental breakup has occurred and oceanic spreading starts. Marine sediments continue to be deposited away from the spreading centre, and a small quantity of coarse continental sediments in the oceanic spreading centre is hot enough to be labelled as source rock (Fig. 8d). However, no favourable conditions for ore formation occur after breakup. On average, less than one mechanism 1 or 2 basin occurs in the narrow symmetric rift evolutions, and the average maximum amount of source rock is only about one-third of the narrow asymmetric source rock (Fig. 7). Source rock area soon returns to zero after oceanic spreading has started (Fig. 6).

Figure 6 (a–c) Source area over time for each simulation plus their arithmetic average for each rift type. Grey areas indicate the period of time when oceanic spreading occurs in the individual simulations. (d) The narrow asymmetric (NA) average is used to compute zinc and lead endowment over time.

3.3 Wide rifts

A typical simulation (W-2) of the initially wide rifting simulations (W-1 to W-9) is presented in Fig. 9. A broader area of deformation leads to lower rift shoulder topography (maximally 868 m) and a shallower but broader marine domain than in the previous simulations. By 10 Myr, five larger faults shoal into the shear zone along the interface of the upper and lower crust. The basins are predominantly filled with sand and shallow marine sediments. 5 million years later, the mantle lithosphere has almost thinned to breakup (Fig. 9d) and the development of source rock starts (Fig. 6c). Soon after, the rift migrates towards the left (Fig. 9e). By this time, two major basins on the deserted margin have heated enough to generate source rock connected to host rock by inactive faults (Fig. 9i). Ore formation mechanism 3 persists at least until 25 Myr, when the rift centre has stabilized again (Fig. 9f and j). The asymmetric margins also exhibit asymmetry in terms of sediment type, with the narrower left margin hosting more fine-grained continental sediments and the wider margin hosting extensive younger marine sediments. The average maximum amount of source rock (∼ 9.4 km²) is again smaller than in the narrow asymmetric rift (Figs. 6 and 7). Additionally, Fig. 7 shows that fewer active faults connect source and host rock than for narrow asymmetric rifts.

3.4 Varying sedimentation and erosion efficiency

Surface process efficiencies can vary significantly as a function of lithologic, environmental, and climatic factors and are often not well constrained. By changing key parameters, we therefore assess their relative impact on model results. Starting from the model setup for narrow asymmetric rifting, which leads to the most prolific ore formation mechanisms (as shown by Fig. 7), we first vary the marine sedimentation rate, then the efficiency of marine and of continental erosion, and finally the extension velocity. We show changes in model results by comparing the total source and host area over time averaged over the nine simulations with different initial plastic strain configurations (Fig. 10).

Figure 7 The average number of basins that contain source rock (red, labelled “S”), source and host rock (orange, S + H), and source and host rock connected by inactive faults (light green, S + H + IF) or active faults (dark green, S + H + AF) per rift type (NA – narrow asymmetric, NS – narrow symmetric, W – wide). The blue line plots the average maximum source area per rift type, which eventually determines the size of the potential mineral resource. The narrow asymmetric rift type produces the most source rock and has the most basins with active faults that provide fluid pathways between source and host rock. All averages are taken over the nine initial strain configurations per rift type.

3.4.1 Varying marine sedimentation rate

Both source and host rock total area increases when the background marine sedimentation rate is increased (Fig. 10a). Source rock starts forming at around 10 Myr in each case, but the period of fastest increase starts earlier and lasts longer the higher the marine sedimentation rate. The approximately linear increase in host rock has a larger slope for higher rates, as would be expected. The higher source rock area for a sedimentation rate of $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{4}}$ m yr⁻¹ is reflected in the occurrence of source rock in every basin at 20 Myr (Fig. 10a top), not just in the basin in the narrow margin. However, these basins only experience ore formation mechanism 3; no active faulting occurs in the wide margin. Ore formation mechanisms 1 and 2 still occur in the major basin on the narrow margin. The overall potential to form deposits has therefore increased.

Figure 8 Model evolution of narrow symmetric rift simulation NS-8 showing continental breakup after 15 Myr and the development of equally wide margins. (a–e) Time series of zoom-ins focusing on FastScape topography (top) and sedimentary basins (bottom). The figures are always horizontally centred around x= 350 km. Black rectangles indicate the locations of further zoom-ins in (f)–(h). (f–h) Zoom-ins of specific basins with and without ore formation mechanisms. (g) Ore formation mechanism 1, where new active faults connect source and host rock in neighbouring basins.

3.4.2 Varying marine erosion efficiency (diffusivity)

Increasing the marine diffusion coefficient has a very limited effect on the amount of host rock, but it leads to an increased generation of source rock in the major basin on the narrow margin (Fig. 10b). A larger diffusivity leads to a faster distribution of sediments away from the continent and therefore to a larger distal sediment load. Source rock area therefore occurs more often in the wide margin (see example in Fig. 10b top). This means that ore formation mechanism 3 could potentially produce larger CD-type deposits than predicted in the reference narrow asymmetric simulations.

3.4.3 Varying continental erosion efficiency

An increased continental erosion efficiency (Fig. 10c) translates to a larger supply of continental sediments to the marine domain. As expected, this is mirrored in a larger quantity of total sediments and a larger amount of source rock. There is also an increase in host rock, but this is relatively smaller. The high sediment supply of very efficient erosion ($K_f=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{5}}$ m^1–2 m yr⁻¹) seems to prolong rift migration: where for the other erosion efficiencies the rift stabilizes in seven and eight of the nine simulations, respectively, in the high-K_f case, only four do within 25 Myr. In the high-K_f case, continental sediment supply during the erosional phase that removes rift shoulder topography is so large that close to the shore no marine sediments are deposited (the accommodation space up to sea level is completely filled by continental sediments). Therefore, the ore formation mechanism 2 phase is shortened or absent in all simulations (Fig. 10c bottom). Mechanism 1, on the other hand, occurs in all five prolonged rift migration simulations around 22.5–25 Myr.

3.5 Varying extension velocity

In a final variation of the reference model of narrow asymmetric rifting, we vary the total applied extension velocity and therefore plot source and host area against total extension instead of time (Fig. 10d). Both source rock and host rock are more abundant when the extension velocity is lower. The lower the velocity, the slower the rift migration, the longer faults are active, and the more marine and continental sediments can be deposited for a certain amount of extension. Thus, the thermal effects of greater burial depths and slower rift migration lead to source rock forming in most basins for an extension velocity of 5 mm yr⁻¹. Total sediment area is reduced for faster velocities, and the faster-migrating rift also prohibits the formation of thick sediment packages. Therefore, no ore formation mechanisms occur for an extension velocity of 20 mm yr⁻¹ (Fig. 10d bottom).

Figure 9 Model evolution of wide rift simulation W-2 showing widespread deformation followed by rift migration and the transition to continental breakup. (a–f) Time series of zoom-ins focusing on FastScape topography (top) and sedimentary basins (bottom). The figures are always centred around x= 350 km. Black rectangles indicate the locations of further zoom-ins in (g)–(j). (g–j) Zoom-ins of specific basins with and without ore formation mechanisms. (i–j) Ore formation mechanism 3, where inactive faults connect source and host rock in the same basins. Note that the host rock in the left basin consists of deep marine sediments (silicates).

4 Discussion

4.1 Predicted spatial and temporal mineralization windows

Our numerical models of continental rifting indicate that asymmetric narrow rifting is most likely to produce favourable conditions for CD-type metallogenesis. This rifting type creates the largest area of source rock and the most basins where ore formation mechanisms 1 and 2 can occur (Figs. 4, 5, and 7). These scenarios occur in the major basin on the narrow margin, where the border fault is in turn active, inactive, and reactivated until rift migration ceases. Fluids can migrate either along the border fault itself or along smaller synthetic and antithetic faults that develop in the major basin as the rift migrates. Prolonged activity of the border fault is important as it creates accommodation space for the deposition of a thick package of both continental and marine sedimentary rock so that the lowermost units of that package reach the temperatures required to form a source rock as they are heated up by rising asthenosphere and are insulated by overlying sediments.

The distributed faulting and shorter-lived narrow faulting in the narrow symmetric rifting case do not allow thick sedimentary sequences to form before oceanic spreading starts around 15 Myr, and only two small patches of source rocks form in one of the nine simulations (Figs. 7 and 8). In five of the wide rifting simulations, a large basin forms once rift migration starts, but soon after source rocks begin to form, activity on its border fault ceases as the rift stabilizes. More common is ore formation mechanism 3, which occurs in the wide margin for eight out of the nine simulations (e.g. Fig. 9e, f, i, j). Our simulations further show that the ore-forming potential depends on a favourable balance between rifting velocities, sedimentation rates, and erosion efficiencies (Fig. 10).

From the above, it is clear that the relative timing of border or synthetic and antithetic fault activity, source rock generation, and the deposition and burial of marine sediments determines the temporal window for mineralization. For example, in the narrow asymmetric simulation in Figs. 3j and k and 4, marine sediments along the border fault have been heated up too much to be able to act as host rock by the time the fault is reactivated. Meanwhile the distal part of the major basin has received less sediment and rotated less along the border fault such that the marine sediments stayed cool enough to act as host rock for the fluids potentially tapped by the new faults. Here, the temporal mineralization window for ore formation mechanism 1 is 3 Myr, and mineralization would occur in the upper several kilometres of sediment within 38–50 km from the shoreline. In the example of ore formation mechanism 2 (Fig. 5a), the border fault is reactivated just in time to connect source and host rock for ∼ 3 Myr before all marine sediments along the fault have heated too much and the fault becomes inactive again. In this case, mineralization would occur closer to shore (∼ 14 km) in the upper few kilometres of sediment. The example simulations for mechanism 1 and 2 each also contain instances of the other mechanism (see Figs. 4 and 5), with a temporal window of 2 and < 1 Myr, respectively.

As discussed in Sect. 2.4, we assume that faults are more permeable than their surroundings and thus form upward-discharging fluid pathways. Recharge of fluids in the ore-producing basins is, however, also important; this can occur downwards along faults. In subaerial parts of the system, a hydraulic head can drive fluid flow and generally follows the water table that mimics the topography; as such, recharge along the main border fault is likely (e.g. Ridley, 2013). However, along-strike, both upward and downward flow may occur along the same fault and therefore both ore formation mechanism 1, with multiple smaller faults cross-cutting the major basin, and mechanism 2, which relies on one larger fault, are plausible mechanisms. Moreover, during both mechanisms other faults are also active in the same basin on the narrow margin (see Figs. 4 and 5), providing additional recharge pathways.

4.2 Comparison to observations and current mineral system concepts

Our physics-based simulations of coupled geodynamic and surface processes can be used to better understand the formation of CD-type deposits in well-studied provinces, which are often complicated by later tectonic overprints. Modelled ore formation mechanisms 1 and 2 reproduce some important characteristics of major Zn–Pb provinces. For one, in the Selwyn Basin (Canada) and the Carpentaria Zinc Belt (Australia), up to several kilometres of sediments can vertically separate the source rocks from the host rocks (Blaikie and Kunzmann, 2020; Hayward et al., 2021; Rodríguez et al., 2021). This agrees well with the 1–3 km between source and host rock at the time of mineralization (20 Myr) for ore formation mechanisms 1 and 2 in Figs. 4 and 5. Additionally, the source sequences themselves are 2–5 km thick in the McArthur (Australia) and Selwyn basins (Cooke et al., 2000; Rodríguez et al., 2021, and references therein), as we find for mechanisms 1 and 2. In the McArthur Basin, Blaikie and Kunzmann (2020) describe the Umbolooga Subgroup as an alternating stratigraphic sequence of siliciclastics, mixed siliciclastics and carbonates, and carbonates, which agrees with the level of stratigraphic detail we can resolve in our simulations. Lastly, many CD-type deposits are considered to have formed as the host unit was undergoing early to mid-burial diagenesis, i.e. at < 2 km depth (e.g. Hnatyshin et al., 2015; Hayward et al., 2021). The predicted depths of mineralization of 1–3.5 km in Figs. 4 and 5 are of similar magnitude.

Figure 10 Total source and host area over time for variations of the narrow asymmetric (NA) rift type, varying (a) marine sedimentation rate R_M, (b) marine diffusion coefficients K_M sand and K_M silt, (c) river incision rate K_f, and (d) extension velocity v_x. Averages of each set of nine initial plastic strain configurations are shown per setup for easier comparison. Note the different ranges on the vertical axes and the total extension that is plotted instead of time in (d). Each graph is accompanied by a zoom-in of the rift basins for a simulation with the same initial strain configuration as the reference simulation in Fig. 3 but with a higher value of the varied parameter instead of the reference value in bold. Increasing the marine sedimentation rate, the marine diffusion coefficients, and the river incision rate increases the total area of source rock at the end of a model run. Final host rock area increases for higher marine sedimentation rates and higher river incision rates, while it decreases for faster extension velocities. No ore formation mechanisms occur for the highest extension velocity.

Fault reactivation or episodic activity along fault zones has been suggested for some districts (e.g. the Selwyn Basin, Lund, 2008). Some known Zn–Pb mineralization is identified as post-faulting, so deposition occurred after the associated fault had become inactive, but syn-rift, as the locus of deformation shifts throughout the rift's evolution (Walsh et al., 2018). In this case, faults still provide fluid pathways, either passively or through fault-valve behaviour driven by fluid overpressure (Walsh et al., 2018). These observations are mirrored by some of our results: the major border fault that acts as a fluid pathway in ore formation mechanism 2 (Fig. 5) shows multiple episodes of inactivity and reactivation, which could provide a trigger for fault-valve behaviour by tectonic events. Furthermore, ore formation mechanism 3 in the wide margin would correspond to the post-faulting mineralization.

There are relatively few reliable temperature constraints for the fluids that formed CD-type deposits. In some localities, homogenization temperatures of fluid inclusions trapped in vein assemblages from the fault-bound feeder zone record a minimum temperature of fluid entrapment between 250 and 300 °C (Cooke et al., 2000; Rajabi et al., 2015; Magnall et al., 2016). These temperatures provide the lowest potential temperature constraints for fluids in the source rock. Within the temporal and spatial mineralization windows of mechanisms 1 and 2 in Figs. 4 and 5, the maximum source rock temperature is 508 and 455 °C, respectively. These temperatures would provide enough heat to the fluids to enter the mineralization site at > 250 °C, where they would quickly cool due to mixing with cooler fluids.

Using Eq. (9), which includes a conversion of source rock area to volume assuming an in-plane extent of the basins of 10 km, the potential endowment for the average of all nine narrow asymmetric simulations can be estimated (Fig. 6d). As the majority of source rock in these simulations is concentrated in the major basin adjacent to the border fault in the narrow margin (see, for example, the source rock distribution in Fig. 3k), the total source rock area is representative of the source rock tapped by faults through mechanisms 1 and 2. These mechanisms were active between 13.5 and 21.0 Myr in Figs. 4 and 5; this time frame corresponds to a possible endowment of up to 35 Mt zinc and 8 Mt lead. This is of the same order of magnitude as the maximum currently estimated zinc tonnage of 27.4 Mt (Red Dog, Alaska) and similar to the maximum lead tonnage of 28.0 Mt (Broken Hill, Australia) reported in the database of Hoggard et al. (2020). Endowment estimates would improve from resolving the relevant processes in 3D and from the inclusion of fluid flow. Future modelling should also investigate the role of thick, cold cratonic lithosphere (e.g. Raghuram et al., 2023; Gouiza and Naliboff, 2021) in promoting the formation and/or preservation of CD-type deposits, as Hoggard et al. (2020) (and Groves and Santosh, 2021; Lawley et al., 2022; Huston et al., 2022) found a correlation between deposit location and the edges of thick and/or cratonic lithosphere.

Our results show the co-occurrence of source rock, host rock, and active faults predominantly in the major basin in the narrow margin of asymmetric rift systems (ore formation mechanism 1 and 2). Possible sites of mineralization are located within several tens of kilometres of the shoreline at the basin margin (Figs. 4 and 5), which supports a number of geological observations. For example, high-salinity fluids derived from seawater evaporation along the basin margin are preserved in fluid inclusions in some locations (Banks et al., 2002; Leach et al., 2004), inspiring the hypothesis of a platform “brine factory” (e.g. Leach et al., 2010). In some of the Paleoproterozoic Australian deposits, mineralization is hosted by units that are transitional between shallow water, evaporitic facies, and deep-water facies (Kunzmann et al., 2019). In terms of reduced sulfur, the metal trap in some of the largest and highest-grade deposits is derived primarily from sulfate-reducing bacteria (e.g. Fallick et al., 2001), which are most active in regions of high carbon burial along basin margins (Kallmeyer et al., 2012). Some of the largest deposits in the North American Cordillera are also hosted in organic-rich units that were deposited in highly productive settings, similar to modern-day upwelling zones (Dumoulin et al., 2014). Trace fossil evidence from the Cordilleran Red Dog system also suggests that host rocks were deposited on the middle to outer shelf (Reynolds et al., 2015). Notably, numerical modelling of mineralizing fluid flow in the Selwyn Basin by Rodríguez et al. (2021) also shows larger areas of metal enrichment towards the higher topography and shallower water of the carbonate platform. Here, cooling by recharging fluids and along the sloping source-layer topography is most efficient. Modelling of the brine reflux system by Manning and Emsbo (2018) shows that evaporated seawater generated on a carbonate platform (the brine factory) can form ore fluids of the required salinity and temperature that sink into the platform, flow oceanward, and discharge along faults to the nearby ocean floor. In short, different lines of evidence point to the optimum alignment of ingredients for CD-type mineralization in the nearshore marine domain.

Our numerical modelling results also indicate that rift migration is key for large CD-type systems. Migration generates a large, deep, hot sedimentary basin on the narrow margin, in which episodic faulting connects source and host rocks, leading to favourable conditions for CD-type metallogenesis. Identifying narrow passive margins of ancient asymmetric rifts could therefore help constrain the exploration space for CD-type deposits. For example, asymmetry of present-day conjugate margins can be found in the Iberia–Newfoundland and central South Atlantic conjugates (Brune et al., 2017). The upper-plate and lower-plate passive margin concept (Lister et al., 1986) has been used to describe the structural asymmetry of such conjugates. In the Brazil and Angola–Gabon margins, the upper-plate margin is narrow with a sharp crustal taper and little hyper-extension, while lower-plate margins exhibit hyper-extended and exhumed domains and a distal supra-detachment (Péron-Pinvidic et al., 2017). Our model predicts that the potential for large endowments via ore formation mechanisms 1 and 2 is highest on the upper-plate side where the continental margin is narrow (i.e. the Brazilian Camamu basin and the Namibe and Benguela basins offshore of Africa), although there are additional secular factors (e.g. seawater chemistry) that might limit margin fertility in modern settings (Wilkinson, 2014). There are not many ancient margins where upper- or lower-plate sections have been identified, reflecting a need for additional geophysical data to characterize the structure of ancient margins. In the ancient western Laurentian margin (North America), however, both plate types have been recognized, and both types host CD-type deposits (Lund, 2008; Hayward and Paradis, 2021). This occurrence of mineralization on both margin types could reflect a combination of all three ore formation mechanisms occurring, with mechanism 3 dominant on the lower-plate side.

5 Conclusions

This study is the first to propose a framework for the first-order identification of favourable conditions for sediment-hosted ore formation systems using coupled numerical models of rift dynamics and surface processes. Using this framework, we identified three mechanisms that combine the main components of the CD-type mineral system (metal source, fluid conduit, and metal trap and host) in continental rifts: (1) smaller, younger faults deform a sedimentary basin and connect the potential source rock and host rock that have formed in the active basin; (2) an active fault provides a pathway from source to host rock that formed in the basin created by this fault; and (3) source rock and host rock are generated in a basin that is no longer active but whose inactive fault could still act as a fluid conduit either through free fluid convection or during fault reactivation.

The time during which and space in which these three ore formation mechanisms are active – the temporal and spatial mineralization window – are small. Mineralization windows occur most often in the narrow continental margin of narrow asymmetric rift systems, where over a time period of less than 3 Myr active faults provide a pathway for fluid to migrate from potential source rocks to potential host rocks. The predicted mineralization then occurs within several tens of kilometres from the shoreline within the top several kilometres of sediment. Temporal constraints on the mineralization window arise from conflicting results of heating and fault activity. For example, the heating of basal coarse continental sediments by rising asthenosphere and continuous deposition of new sediments that blanket the units below generates source rocks but simultaneously heats up shallower marine sediments too much to act as host rocks.

Rift migration, especially during narrow asymmetric rifting, promotes the formation of one deep basin along the persistent border fault on the narrow margin. Its thick sedimentary infill is continuously heated from below and generates large amounts of source rock that could produce metal endowment of similar magnitudes as the world's largest deposits. Moreover, all three ore formation mechanisms occur in the narrow margin during rift migration, with the border fault acting as a fluid conduit in mechanisms 2 and 3 and the sequential faults of the migrating rift centre providing fluid pathways according to mechanism 1. The spatiotemporal window of mineralization during rift migration is modulated by the rate of sediment erosion and deposition as well as the extension velocity, stressing the importance of coupled geodynamics and surface processes for a holistic understanding of sediment-hosted mineral systems. The amount of source rock produced in narrow symmetric rifts and wide rifts is only one-third of that generated in narrow asymmetric rifts, and favourable conditions for ore formation occur less often. Therefore, exploration for CD-type deposits could benefit from the identification of ancient narrow margins.