SULFIDE TEXTURES AND ORE GRADES IN MINERALIZED CARBONATES DEPEND ON PÉCLET AND DAMKÖHLER NUMBERS

The sustainable development of mineral resources is one of the major challenges facing society and the transition toward a green, digital, and circular economy (Gonzalez-Alvarez et al., 2020). Our highest-value base metal (Zn, Pb) resources are associated with clastic-dominated (CD) deposits, which comprise thin (10¹ m), laterally extensive (10³ m), and stratiform to stratabound sulfide mineralized rocks within ancient sedimentary basins (Leach et al., 2005). The high grades, moderate to high tonnages, and potential for coproduction of critical raw materials (Ge, Ga, In) mean that CD-type deposits are particularly attractive exploration targets to meet increased resource demand.

In the Carpentaria zinc province of northern Australia (Fig. 1A), nearly 30% of the world’s Zn resources are hosted in a number of world-class CD-type deposits (Large et al., 2005). The orebodies of CD-type deposits have relatively simple sulfide assemblages (sphalerite, galena, and pyrite) and are hosted by mudstones, siltstones, or mixed siliciclastic and carbonate sequences. Mineralization textures in CD-type deposits can be highly variable, leading to debate over genetic models (Perkins and Bell, 1998; Large et al., 2005). For example, stratiform to stratabound relationships between sulfides and the host rock have been widely used to infer exhalative processes and a synsedimentary origin for ore-stage mineralization (e.g., Fig. 1B; Goodfellow and Franklin, 1993; Large et al., 1998). In contrast, other studies have described how host-rock replacement in the subsurface can produce a variety of textures, including stratiform to stratabound relationships (Williams, 1978; Chen et al., 2003; Kelley et al., 2004; Magnall et al., 2020; Spinks et al., 2021). Nevertheless, the physical flow parameters that controlled these types of mineralization textures and grade distribution are poorly known (Sheldon et al., 2023). Recent studies have now demonstrated how CD-type mineralization likely formed via carbonate replacement (Fig. 1C) at two localities in the Carpentaria province (Magnall et al., 2021; Spinks et al., 2021).

Carbonate dissolution, such as that required for the carbonate-replacement model, is a widely studied phenomenon because of its significance in radioactive waste storage, enhanced oil and gas recovery, and CO₂ sequestration (Schechter and Gidley, 1969; Fredd and Fogler, 1998; Ghommem et al., 2015). These studies provide a wealth of knowledge on fluid-rock interaction in sedimentary units that has not been widely applied to economic geology. In particular, the acidification of oil wells provides a useful analogue for better understanding the effects of reactive flow in carbonate replacement that form CD-type deposits. Mineral dissolution studies (Daccord, 1987; Fredd and Fogler, 1998; Steinwinder and Beckingham, 2019) have shown that carbonate dissolution patterns are controlled by velocity of the fluid (V), characteristic length of the flow path (i.e., nodal length of the model, L), reaction rate (K), diffusion coefficient (D), surface area (S), and discharge rate (Dc). Together, these parameters can be described by Péclet (Pe, eq. 1) and Damköhler (Da, eq. 2) numbers, as follows:

Importantly, Péclet and Damköhler numbers are unitless, meaning observations are scale independent and can be extrapolated from thin section to drill core to outcrop/deposit scales. When acid is injected into a carbonate unit, the following carbonate dissolution patterns have been described (Fig. 2; Golfier et al., 2002): (1) uniform dissolution in which dissolution occurs evenly throughout the domain; (2) wormholes that are formed through concentrated flow and dissolution along preferential flow paths; and (3) face dissolution, where a sharp dissolution front progresses out from an inlet. Here, we investigate how different types of reactive flow regimes might explain the variability in mineralization textures of carbonate-replacement CD-type systems and constrain the optimum conditions for high-grade ore formation. Although this work uses a CD-type deposit as a basis, the model is fairly general, and implications of the results transfer easily to other carbonate-replacement systems such as Mississippi Valley-type and Irish-type as well as ore deposition in noncarbonate units.

Recent work on the Teena and McArthur River (HYC) deposits in the Carpentaria zinc province (Fig. 1) has documented a range of styles of sphalerite mineralization (Fig. 3) that has been attributed to carbonate replacement (Magnall et al., 2021; Spinks et al., 2021). Hydrothermal experiments involving metal-bearing brines and carbonates have also been used to simulate ore formation and further constrain fluid temperature and composition (Liu et al., 2021). In this study, reactive flow models in which a metal-bearing brine is injected into carbonate were used to simulate different styles of host-rock dissolution in a simplified CD-type system analogous to Teena and McArthur River (HYC). We then examined the impact of varying the Damköhler and Péclet numbers on dissolution textures and sphalerite concentration/ore grade. Although the simulations represent a simplification of the natural systems, they are an entirely new and novel approach for interpreting mineralization textures in the context of fundamental physical flow parameters, which we hope will stimulate discussion and future research in this exciting field.

The Carpentaria zinc province comprises the Proterozoic Isa superbasin, which includes the Mt. Isa and McArthur basins (Fig. 1). The Mt. Isa basin has undergone extensive deformation and metamorphism, whereas the McArthur basin is relatively undeformed and not metamorphosed. As a result, the most well-preserved CD-type deposits are located in the McArthur basin (Large et al., 2005). The stratigraphy in the McArthur basin is dominated by a 5- to 15-km-thick succession of mixed carbonate and siliciclastic rocks that are intercalated with bimodal volcanic rocks, which are a potential metal source (Plumb, 1979; Jackson et al., 1987; Cooke et al., 1998; Southgate et al., 2000).

The McArthur basin consists of a series of intercontinental sediments, which were deposited during a period of extension in fault-bound depocenters between 1.82 and 1.58 Ga (Southgate et al., 2000). The Teena and McArthur River (HYC) deposits are hosted by the 1.64 Ga Barney Creek Formation. This formation varies from shallow-water dolostones to fine-grained silt- and mudstones (Kunzmann et al., 2022). Sulfide deposits are primarily hosted by dolomitic silt- and mudstones that are locally enriched in organic matter and diagenetic pyrite (McGoldrick et al., 2010; Hayward et al., 2021). The mineralization occurs at the base of the HYC Pyritic Shale Member and the top of the W-Fold Shale Member. The W-Fold Shale Member is dominated by thinly interbedded dolostones with cycles of arkosic sandstones and dolomitic siltstones (Hayward et al., 2021). This unit transitions into the overlying Lower HYC Pyritic Shale Member that consists of talus slope breccias and thinly bedded to laminated carbonaceous siltstones and mudstones (Spinks et al., 2021).

Total inferred resources at Teena are 58.0 million tonnes (Mt) at 11.1% Zn and 1.6% Pb (Rox Resources, 2016) and the premining resource at HYC was 227 Mt at 9.2% Zn and 4.1% Pb (Porter, 2017). The sulfide mineralization is interpreted to have formed when an oxidizing metal-bearing brine migrated up subbasin-bounding extensional faults before reacting and moving laterally through the host unit, dissolving carbonate (Magnall et al., 2021; Spinks et al., 2021).

Two separate models (1-D and 2-D; Fig. 4) were created for this study using Geochemists Workbench (GWB). In the resulting simulations, brine flowed from left to right and dissolved dolomite while depositing sphalerite and hematite:

Dolomite dissolution buffered pH and drove sphalerite precipitation, whereas the ferrous iron to hematite reaction provided a redox buffer and generated acidity. The thermodynamic data used was the thermo.com.V8.R6+ data set as included with GWB (Johnson et al., 1992). In both the 1-D and 2-D models, we varied the rate of dolomite dissolution as well as the flow rate to change the Damköhler and Péclet numbers between simulations. The amount of fluid flowing through the simulations was 500 pore volumes, which was enough fluid to significantly alter the model domain while not being such a high volume that it would cause model stability issues. This constant volume of flow was simulated to pass through the model at different rates over different lengths of time. So, although individual simulations had different flow rates, the amount of fluid passing through each was always the same.

Representative mineralogy of relatively unaltered rock at the Teena deposit is 56% dolomite, 7% feldspar, 13% quartz, 9% pyrite, and 15% chlorite, mica, and clays based on previous work (Magnall et al., 2021). We took an average of the dolomite concentration and excluded other minerals from the model to simplify the interpretation of the Damköhler number. Mineral reactions in the model were at equilibrium except for the dissolution of dolomite, which was kinetically constrained based on a published rate constant and activation energy (Palandri and Kharaka, 2004). The base model used a kinetic rate constant for dolomite of 1.7 × 10^–8 mol/cm²/s, but this was decreased by up to an order of magnitude to vary the Damköhler number. Because the altered minerals at Teena are dominated by fine silt and clays and lack fluid inclusions, the brine constraints used here were based on a mixture of assumptions about mineral equilibrium and extrapolations from other ore-forming brines at approximately the same salinities (Table 1; Hanor, 1996; Cooke et al., 2000; Kharaka and Hanor, 2003).

The 1-D model was 100 m long with 20 nodal grids and was run using X1t (Bethke et al., 2022). Using a 1-D model allowed us to maintain a constant Péclet number for the model while varying the Damköhler number so the effect of two values could be assessed separately. At the end of each simulation, we calculated the Damköhler number of each node.

We created 2-D models using X2t (Bethke et al., 2022) to simulate dissolution patterns. The dimensions of these models were 1 × 1 m with 25 nodes on each axis and matched the input conditions of the 1-D models. Again, these models had varying flow and dolomite dissolution rates in order to vary the Péclet and Damköhler numbers. This approach allowed for comparison of the simulations to the observed textures in the rocks. The initial porosity averaged 10% based on a reasonable estimate for a siltstone (Freeze and Cherry, 1979) with a standard deviation of 2.5%. The chosen permeability was dependent on the porosity based on a log-linear relationship (Bethke, 1985) with an average value of 3.1 × 10^–16 m². The random porosity and permeability field allowed for simulations to potentially generate preferential flow paths. Preferential flow paths are generated when mineral dissolution creates zones of high permeability, which then creates a feedback where more flow through the high-permeability zone leads to more dissolution and higher permeability along that path. At a constant flow rate, the proportion of fluid flowing along a flow path is dependent on the ratio of the permeability along that flow path relative to other potential flow paths; therefore, the ratio of permeabilities between nodes was significant, whereas their absolute permeability values were not. The complete input for both models is given in the Appendix.

Different ranges of values for the Péclet and Damköhler numbers describe types of feedback between fluid flow and carbonate dissolution (Fig. 2). For example, decreasing the Damköhler number means the reactants begin to flow into the system faster than the reactions can take place; decreasing the Péclet number shifts the balance between advection and diffusion in the system toward diffusion.

For the 1-D model, the Damköhler number varied depending on where a node is located in relationship to the reaction front. Ahead of the reaction front, fluid flowing through the system was much closer to equilibrium with the dolomite, leading to a slow rate of dolomite dissolution and a low Damköhler number. The area ahead of the front was, therefore, undergoing uniform dolomite dissolution and precipitation of disseminated sphalerite (Fig. 3A). The Damköhler number was higher at the reaction front, which placed the model within the wormhole dissolution domain (Fig. 2) and led to a gradual contact between ore lenses and host rock. The trend of ore grade vs. Damköhler number was dependent on the Péclet number of the system and vice versa. Mutual dependency of grade on the Damköhler and Péclet numbers is because the precipitation of sphalerite was driven by dolomite dissolution and as either number shifted, dissolution morphology was still partially dependent on the value of the other number.

In the 2-D, smaller-scale simulations, carbonate dissolution followed the predicted patterns for the different Péclet and Damköhler numbers with uniform dissolution in areas having a low Damköhler number. Where the Damköhler number was higher, face or wormhole dissolution was dependent on the Péclet number (Fig. 2) because the reaction fronts become unstable as the velocity increases (Chadam et al., 1986). As dolomite dissolution produced sphalerite precipitation, the model results yield ore grades that are essentially the inverse of dolomite concentrations. In summary, face dissolution created a massive deposit (Fig. 3A), wormholes created an interfingering deposit (Fig. 3C), and uniform dissolution created a disseminated deposit (Fig. 3E).

The subjective nature of mineralization textures in CD-type deposits has resulted in debate over paragenetic relationships and a poor understanding of paleofluid flow within these systems (Logan et al., 1990; Large et al., 1998; Magnall et al., 2021; Spinks et al., 2021). Importantly, results of the simulations are comparable to different dissolution styles observed in the Teena deposit, thereby providing a framework for interpreting styles carbonate-replacement deposits in the context of flow rate and mineral kinetics.

Disseminated mineralization, particularly where confined to specific beds (i.e., stratiform mineralization; see Fig. 3F), is commonly attributed to the earliest stages of ore formation and is sometimes used to support synsedimentary or exhalative models of mineralization (Sangster, 2018). However, it is also possible to produce such textures through uniform dissolution, which occurs at low Damköhler numbers where slow mineral kinetics or fast fluid velocities lead to uniform dissolution of carbonate minerals throughout a formation. Coarser mineralization (e.g., Fig. 3D) is commonly interpreted to occur in later stages of ore formation by telescoping of the system or as the result of multiple superimposed pulses of fluid flow and progressive host-rock dissolution (Hayward et al., 2021). However, these textures could also be generated through the development of wormholes, where dissolution is concentrated around preferential flow paths because of fast mineral dissolution and fluid velocities.

There is a correlation between Damköhler and Péclet numbers with grade in our models (Fig. 5). Low Péclet and high Damköhler numbers describe a system that is more likely to have concentrated sphalerite deposition, resulting in high ore grades. The smoothness of the correlation is due to the gradational boundaries between the different dissolution morphologies. In simulations with lower Péclet numbers, the dissolution of dolomite and precipitation of sphalerite was more focused, leading to higher ore grades. The low Péclet number models had face dissolution/massive precipitation and yielded the highest sphalerite concentrations (Fig. 3B). In systems with increasing Damköhler numbers, precipitation is increasingly concentrated around the inlet and preferential flow paths. The wormholes and the preferential flow paths they represent are the type of reaction-generated permeability that recent work suggests is important to accurately model ore deposition (Magnall et al., 2023; Sheldon et al., 2023). Correlation of higher grades with increasing Damköhler number was due to the shift away from disseminated mineralization. Therefore, systems with low Péclet and high Damköhler numbers are able to more efficiently strip metals out of brine flowing through the unit.

A key factor controlling dissolution efficiency and morphology is mineralogical heterogeneity (Wei et al., 2019). This issue raises the possibility of defining the mineralogical “Goldilocks zone” needed for the development of high Damköhler numbers and high grades, which would be a key parameter to include in exploration models. The kinetics of mineral dissolution and reactive surface area are main factors controlling the Damköhler number and, by extension, mineralization style and ore grades. For example, the rate of mineral dissolution increases at higher temperature; therefore, in a prograde Cu skarn more massive mineralization should form closer to the heat source—a pattern that is seen in nature (Atkinson and Einaudi, 1978). In some carbonate-hosted systems, however, this temperature dependence could be offset by the retrograde solubility of the host rock. The mineralogy of the host formation is also a control on the Damköhler number because of the faster reaction kinetics of, for example, calcite compared to dolomite (Hoefner and Fogler, 1988; Wei et al., 2019). The high surface area of certain types of diagenetic pyrite (e.g., framboidal morphologies) may also be an important factor in enhancing Damköhler numbers. Fast-reacting pyrite in silicate rich units, such as in Figure 3B, could enhance the reactivity of an otherwise slowly reacting unit and increase the Damköhler number.

The carbonate system shifts closer to equilibrium with increasing pH or the concentration of a dissolved ion such as calcium, which not only would decrease the amount of carbonate dissolution but also would lower the Damköhler number and shift the ore precipitation pattern toward a disseminated deposit. This is why there was disseminated mineralization ahead of the main reaction front (Fig. 3A). Each of the various styles of mineralization can occur at different points in the same flow regime, i.e., disseminated deposition ahead of a reaction front that could be characterized as interfingering or massive. Although the system presented here is relatively simple, most systems contain multiple reacting phases, each of which could be characterized by different Damköhler and Péclet numbers. Further modeling in conjunction with experimental results would help constrain the effects of different types of heterogeneity and the feedback between fluid flow and reaction rates.

The relationship between Péclet and Damköhler numbers and ore grade demonstrates how large volumes of metal-bearing brines could flow through a suitable host rock but not yield an economic deposit if the relationship between flow rate, diffusion, and mineral kinetics (Pe and Da) in the system is not optimal, such as is shown in Figure 3E where the kinetics were too slow. Specific flow rates cannot be given here, because the exact relationship between flow rate and grade is simulation dependent, owing to the different permeability/porosity relationships implemented in various models (Kühn, 2004). Incorporating the concept of flow rates into exploration strategies would likely require construction of 3-D basin models (Garven et al., 2001; Yang et al., 2006; Sheldon et al., 2021, 2023). These basin models could show which regions within a basin have experienced the optimal flow rate for ore deposition and, therefore, reduce the search area for exploration.

Using mineralization textures as a proxy for identifying low Péclet/high Damköhler domains and mapping how these textures correspond with other drill core information (e.g., host-rock facies) may help to further constrain highly reactive flow zones and parameterize numerical models. Incorporating facies variability into the reactive flow model would also help investigate the role of different types of heterogeneity on the system. For example, one testable hypothesis is that phases that enhance reactive surface area (e.g., diagenetic pyrite) are likely to result in Damköhler hot spots and further enhance ore grade. Another possibility is that the development of wormholes could explain some of the spectacular dissolution breccias formed by hydrothermal karstification processes (Corbella et al., 2006; Petrus and Szymczak, 2016; Wang et al., 2017) as well as other textures such as stylolites. Experimental work could also determine the thresholds between morphology domains and relevant Péclet and Damköhler numbers for particular systems; importantly, because these numbers are dimensionless, results from lab experiments could be used to inform deposit-scale processes.

Finally, these results are not limited to the systems in the Carpentaria province. Carbonate replacement is an important ore-forming process, not only in CD-type deposits but also in Mississippi Valley-type and Irish-type Pb and Zn deposits, stratiform Cu deposits, Au-mineralized Carlin systems, and higher-temperature systems such as carbonate-replacement deposits (mantos) and skarns. Similar reaction patterns could be seen in other mineral replacement systems or with redox-controlled deposition (Fig. 3B; Tranter et al., 2021; Jones et al., 2023). Across a range of deposit types, this modeling provides an entirely new and novel approach for interpreting mineralization textures in the context of fundamental physical flow parameters.

Laterally extensive orebodies are one of the most characteristic, but also enigmatic, features of sediment-hosted mineral systems. A number of important deposit types within this mineral system formed via carbonate replacement. Interpreted within the context of carbonate dissolution and matrix acidification, the variety of observed mineralization textures and geometries can be related to key properties of reacting flow systems that are described by Péclet and Damköhler numbers.

The variation of Péclet and Damköhler numbers can be used to explain the differences between deposit styles in carbonate-replacement ore systems. These results emphasize the need for understanding both the flow rates and the mineral kinetics of CD-type systems through incorporating the evolution of fluid flow rates over time into exploration strategies, which would likely require the expansion of basin-modeling efforts. Face dissolution leads to massive deposits, wormholes lead to interfingering and lenses, and uniform dissolution creates disseminated deposits. The relationship between Péclet and Damköhler numbers and deposit style leads to the correlation between these numbers and ore grade; more compact, higher-grade mineralization occurs at higher Damköhler and lower Péclet numbers.

Peter M. Berger, Joseph M. Magnall, and Sarah A. Gleeson are funded by a Helmholtz-Rekrutierungs initiative. We would like to thank reviewer Sam Spinks and Associate Editor John Slack for their valuable comments and suggestions.

Conflict of interest: The authors declare no competing interests.

Peter Berger is a postdoctoral researcher at the GFZ Hemholtz Centre for Geosciences (Germany). He obtained his B.S. degree at Bowling Green State University and M.S. degree at the University of Illinois Urbana-Champaign in geology. During his M.S. studies, he worked as a programmer on reactive transport modeling software. He then spent 10 years working at the Illinois State Geological Survey on carbon sequestration projects doing a mixture of modeling, field, and experimental work. Afterward, Peter got his Ph.D. degree at the University of Tasmania (Centre for Ore Deposit and Earth Sciences; CODES) and shifted to economic geology. Peter is currently working on numerical models of ore deposition in sediment-hosted systems.