Hostname: page-component-7c8c6479df-xxrs7 Total loading time: 0 Render date: 2024-03-28T20:05:46.874Z Has data issue: false hasContentIssue false

Outcrop scale mixing enhanced by permeability variations: the role of stationary and travelling waves of high saturation indices

Published online by Cambridge University Press:  13 October 2022

Daniel Koehn*
Affiliation:
GeoZentrum Nordbayern, Friedrich Alexander University (FAU) Erlangen-Nuremberg, Schlossgarten 5, Erlangen 91054, Germany
Ulrich Kelka
Affiliation:
CSIRO – Mineral Resources, 26 Dick Perry Ave, Kensington, WA 6151, Australia
Renaud Toussaint
Affiliation:
University of Strasbourg, CNRS, Institut Terre et Environnement de Strasbourg, UMR 7063, 5 rue René Descartes, Strasbourg Cedex 67084, France SFF PoreLab, the Njord Centre, Department of Physics, University of Oslo, PO Box 1048 Blindern, NO-0316 Oslo, Norway
Coralie Siegel
Affiliation:
CSIRO – Mineral Resources, 26 Dick Perry Ave, Kensington, WA 6151, Australia
Gary Mullen
Affiliation:
School of Geographical and Earth Sciences, Glasgow University, G12 8QQ, UK Scottish Universities Environmental Research Centre, East Kilbride, Glasgow, G75 0QF, UK
Adrian Boyce
Affiliation:
Scottish Universities Environmental Research Centre, East Kilbride, Glasgow, G75 0QF, UK
Sandra Piazolo
Affiliation:
School of Earth and Environment, University of Leeds, LS2 9JT, UK
*
Author for correspondence: Daniel Koehn, Email: daniel.koehn@fau.de
Rights & Permissions [Opens in a new window]

Abstract

To study the ore mineralization at the outcrop scale we merge an advection–diffusion simulation with the geochemical software iphreeqc to model the mixing of two realistic fluids. We simulate the infiltration of a metal-rich fluid into a rock that is saturated with pore fluid. We test the feedback effects with a number of scenarios based on an outcrop-scale 5 × 5 m model consisting of two high-permeable vertical faults within a low-permeable host rock that lead into a permeable layer. The hot metal-rich fluid enters the model through the faults from below. We solve the advection–diffusion equation for 12 chemical species and temperature, and use iphreeqc to determine the resulting properties of local fluid domains as well as related saturation indices for minerals. The faults in the model act as pathways for the metal-rich fluid, with the infiltrating fluid displacing the pore fluid. Mixing in the model takes place as a function of advection along permeable faults coupled with diffusion of chemical species at the interface between two fluids, while heat diffusion is fast enough (103 times faster) to equilibrate temperature. Simulations show a high saturation index of mixing-derived minerals such as barite at the interface between the two fluids as a result of fluid mixing. Fast fluid pathways (i.e. faults) show travelling waves of high saturation indices of barite, while low-permeability zones such as fault walls and areas below less permeable layers experience stationary waves of high saturation indices. Our results show that, depending on the dominant transport process (advection or diffusion), mineralization will localize next to permeability contrasts in zones where local diffusion dominates.

Type
FLUID FLOW AND MINERALIZATION IN FAULTS AND FRACTURES
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© The Author(s), 2022. Published by Cambridge University Press

1. Introduction

Fluid mixing or unmixing is suggested to be one of the dominant ore-forming processes that lead to the emplacement of economic hydrothermal precious and base metal mineralization (Haynes et al. Reference Haynes, Cross, Bills and Reed1995; Fu et al. Reference Fu, Williams, Oliver, Dong, Pollard and Mark2003; Fan et al. Reference Fan, Hu, Yang and Wang2006; Neumayr et al. Reference Neumayr, Walshe, Hagemann, Petersen, Roache, Frikken, Horn and Halley2008; Gessner et al. Reference Gessner, Kühn, Rath, Kosack, Blumenthal and Clauser2009; Leach et al. Reference Leach, Taylor, Fey, Diehl and Saltus2010; Lester et al. Reference Lester, Ord and Hobbs2012; Hobbs & Ord, Reference Hobbs, Ord, Gessner, Blenkinsop and Sorjonen-Ward2018), but is still not well understood (Ord et al. Reference Ord, Hobbs, Walshe, Metcalfe, Trefry, Zhao, Lester, Regenaur-Lieb, Rudman, Schmidt and Brugger2010; Lester et al. Reference Lester, Ord and Hobbs2012). One important characteristic of hydrothermal ore-forming systems is the development of pronounced spatial gradients in chemical potential (e.g. pH/redox fronts) that result from the mixing of two or more fluids of different hydrochemical character and origin (Lester et al. Reference Lester, Ord and Hobbs2012), where changes in redox potential, O2 fugacity, pH and temperature can trigger mineral precipitation (Liu et al. Reference Liu, Spinks, Glenn, Macrai and Pearce2021; Weihua et al. Reference Weihua, Spinks, Glenn, MaCrae and Pearce2021). Mixing of fluids with different chemical composition alone can result in supersaturation with respect to the gangue and ore minerals at conditions far from equilibrium with the host formation (Schwinn et al. Reference Schwinn, Wagner, Baatartsogt and Markl2006). Following Lester et al. (Reference Lester, Ord and Hobbs2012), fluid mixing itself can take place by (1) pure diffusion, (2) physical differences and thus instabilities between two fluids or (3) a combination of advection along local geological pathways and diffusion. For larger deposits diffusion alone as a mixing mechanism is thought to be too slow, so effective mixing should contain an advective part that leads to folding or fingering of fluid streams so that diffusive length scales are reduced significantly (Lester et al. Reference Lester, Ord and Hobbs2012). Here we perform coupled outcrop-scale numerical simulations with geological features that include sedimentary layers and faults for two realistic fluids to show that permeability variations can profoundly enhance mixing and thus mineralization.

Several studies highlight the importance of permeability contrasts and thus geological geometries on ore localization. One example is the unconformities between basement and overlying sedimentary basins, such as the German Schwarzwald of the European Variscan Orogen (Fusswinkel et al. Reference Fusswinkel, Wagner, Wälle, Wenzel, Heinrich and Markl2013; Bons et al. Reference Bons, Fusswinkel, Gomez-Rivas, Markl, Wagner and Walter2014), where fluid is thought to come up in pulses and ores will precipitate close to an unconformity (Wagner et al. Reference Wagner, Okrusch, Weyer, Lorenz, Lahaye, Taubald and Schmitt2010). At the Finnish Outokumpu Cu-Zn-Co deposit, black shales form a cap rock above the ore-body and mineralizing fluids are suggested to have circulated below this cap forming economic deposits (Loukola-Ruskeeniemi, Reference Loukola-Ruskeeniemi1999). Fluids involved in ore deposits are often interpreted as being a basement-derived metal-bearing brine and metal-poor formation water (e.g. Fusswinkel et al. Reference Fusswinkel, Wagner, Wälle, Wenzel, Heinrich and Markl2013; Bons et al. Reference Bons, Fusswinkel, Gomez-Rivas, Markl, Wagner and Walter2014). Structures that may represent fluid mixing and associated mineral precipitation can be observed across scales from the wall rock of faults (e.g. Hollinsworth et al. Reference Hollinsworth, Koehn and Dempster2019), inside of veins (Fig. 1a, b) and at the scale of hand-specimen inside small-scale fracture networks. In addition, mineral precipitation adjacent to high-permeability fractures can be observed in laboratory experiments (Fig. 1c). Determining the fluid compositions that lead to hydrothermal mineralization is often challenging as only limited data can be obtained, for instance from fluid inclusions (Wilkinson, Reference Wilkinson2001; Fusswinkel et al. Reference Fusswinkel, Wagner, Wälle, Wenzel, Heinrich and Markl2013). Quantifying the effects of fluid mixing on ore deposition represents an additional major challenge as the impact and nature of permeability structures that lead to localization of economic mineralization remains poorly understood (Ord et al. Reference Ord, Hobbs, Walshe, Metcalfe, Trefry, Zhao, Lester, Regenaur-Lieb, Rudman, Schmidt and Brugger2010; Lester et al. Reference Lester, Ord and Hobbs2012).

Fig. 1. Example of mineralization at different scales within and adjacent to high-permeability sites. (a) Mineralization associated with faulting, Trollers Gill Barite–Fluorspar–Galena mine, Yorkshire Dales, UK. (b) MicroXRF geochemical mapping image acquired with a Bruker M4 Tornado from the hole NDIBK04 drilled by the National Drilling Initiative through the MinEx CRC. NDIBK04 is the most prospective well in the East Tennant Creek (Northern Territory) NDI region, an emerging mineral province for base metal mineral exploration (Zn, Cu and As). This sample is from the top of a c. 150-m-thick banded graphitic schist of the Palaeoproterozoic Alroy Formation intersected at depth 236.2–237 m. The nearest geochronology available is from a sample at 301–307 m depth, giving a maximum depositional age of 1969 ± 10 Ma (dated by SHRIMP at Geoscience Australia with zircons). The main vein has gypsum at its centre with Fe–Mn-rich carbonate to the edge. Such vein geometries can arise from the interaction with different fluids or by the evolution of a single infiltrating fluid with the host rock. The presence of gypsum in the centre of the vein points towards an open channel that prevailed after the diffusion-dominated reaction between fluid and rock that led to the development of the reaction rims surrounding the veins. (c) Experiment of calcite precipitation in low-porosity dolomite directly adjacent to high-porosity fracture in response to carbonate-rich fluid infiltration (courtesy of R. Morgan and Q. Fisher).

To address the knowledge gap in the effect of permeability variations, fluid mixing and transport mechanisms on the localization of ore mineralization, we investigate the effect of permeability structures on fluid mixing and mineral saturation states in simple two-dimensional (2D) numerical models on the scale of the geological geometries, that is, the outcrop scale, with two realistic fluids that contain all the important chemical components. We present a numerical framework that couples fluid transport in porous media (elle, Koehn et al. Reference Koehn, Piazolo, Beaudoin, Kelka, Spruženiece, Putnis and Toussaint2021) with the hydrogeochemical simulation module iphreeqc (Charlton & Parkhurst, Reference Charlton and Parkhurst2011). Similar approaches in creating reactive transport modelling environments were reported decades ago (Steefel & Lasaga, Reference Steefel and Lasaga1994) and were performed by utilizing phreeqc as a reaction rate engine that is coupled with UTCHEM (Korrani et al. Reference Korrani, Sepehrnoori and Delshad2015), MATLAB (Shabani et al. Reference Shabani, Kalantariasl, Abbasi, Shahrabadi and Aghaei2019), COMSOL (Nardi et al. Reference Nardi, Idiart, Trinchero, de Vries and Molinero2014; Wei & Cao, Reference Wei and Cao2021) and POET (De Lucia et al. Reference De Lucia, Kühn, Lindemann, Lübke and Schnor2021), as well as other reactive transport codes (see Steefel et al. Reference Steefel, Appelo, Arora, Jacques, Kalbacher, Kolditz, Lagneau, Lichtner, Mayer, Meeussen and Molins2015 for a review). The evolution of permeability in porous media and rocks has been modelled extensively (Saripalli et al. Reference Saripalli, Meyer, Bacon and Freedman2001; Zhao et al. Reference Zhao, Hobbs, Hornby, Ord, Peng and Liu2007; Jamtveit et al. Reference Jamtveit, Putnis and Malthe-Sørenssen2009; Chen et al. Reference Chen, Kang, Carey and Tao2014; Kang et al. Reference Kang, Chen, Valocci and Viswanathan2014; Mostaghimi et al. Reference Mostaghimi, Liu and Arns2016) with smooth particle hydrodynamics, lattice Boltzmann methods and computational fluid dynamic techniques (Manwart et al. Reference Manwart, Aaltosalmi, Koponen, Hilfer and Timonen2002; Fredrich et al. Reference Fredrich, Digiovanni and Noble2006; Tartakovsky & Meakin, Reference Tartakovsky and Meakin2006; Shabro et al. Reference Shabro, Torres-Verdin, Javadpour and Sepehrnoori2012; Chen et al. Reference Chen, Kang, Robinson, He and Tao2013). There has been significant development in numerical models of ore deposits in the last 30 years, mainly focusing on large-scale simulations coupling geological geometries with fluid flow, heat transfer, species transport and reactions as well as deformation (e.g. review in Gessner et al. Reference Gessner, Kühn, Rath, Kosack, Blumenthal and Clauser2009, Reference Gessner, Blenkinsop and Sorjonen-Ward2018; Hobbs & Ord, Reference Hobbs, Ord, Gessner, Blenkinsop and Sorjonen-Ward2018; Moosavi et al. Reference Moosavi, Kumar, De Wit and Schroter2019). In this contribution we aim to study the details of the outlined processes on the metre scale to understand the local influence of geological structures on fluid mixing. To do this we use geologically realistic outcrop-scale scenarios including a formation fluid and an incoming metal-rich fluid. While the location of high permeability areas remains static, the presented approach can be readily adapted to include fracture dynamics.

2. Methodology

2.a. Overview

In this contribution we investigate mixing between formation equilibrated pore fluid and an incoming metal-rich fluid pulse in an outcrop-scale modelling domain of 5 × 5 m with vertical faults and horizontal layers. We couple the hydrodynamic advection–diffusion model “latte” within the microstructural modelling software package elle for the basic simulations of fluid transport and mixing (Bons et al. Reference Bons, Koehn and Jessell2008; Ghani et al. Reference Ghani, Koehn, Toussaint and Passchier2013, Reference Ghani, Koehn, Toussaint and Passchier2015; Piazolo et al. Reference Piazolo, Bons, Griera, Lloens, Gomez-Rivas, Koehn, Wheeler, Gardner, Godinho, Evans, Lebensohn and Jessell2019; Koehn et al. Reference Koehn, Piazolo, Sachau and Toussaint2020, Reference Koehn, Piazolo, Beaudoin, Kelka, Spruženiece, Putnis and Toussaint2021) with iphreeqc (Charlton & Parkhurst, Reference Charlton and Parkhurst2011). In the simulations, the solid that represents a geological zone with its specific porosity and permeability, that is, rock type (chemical composition), is represented by a lattice spring network (Hrennikoff, Reference Hrennikoff1941). It should be noted that in the presented scenarios the elastic properties of the lattice spring network is not utilized; only the lattice geometry is utilized, allowing for porous flow. The lattice network is coupled to a finite difference grid (Press et al. Reference Press, Flannery, Teukolsky and Vetterling1992) that is used to determine the evolution of fluid pressure as well as advection and diffusion of temperature and the chosen necessary 12 chemical species present in the two fluids. Each node in the finite difference fluid grid can have a different chemical composition depending on the mixing fluids. This information is then passed to iphreeqc, where the pH and density of each fluid node is calculated, as well as saturation indices of minerals to understand where minerals precipitate.

2.b. General model set-up

2.b.1. Solid geometry and fluid movement

Table 1 provides a list of all model parameters with notation and units. The solid is represented by a lattice spring network with a triangular mesh of cells that stores the rock properties, that is, porous media (porosity, mineral composition) and the geometry (faults, layers). This network is then connected to a square fluid grid where the transport is simulated assuming porous flow along a pressure gradient. The solid passes a local porosity ϕ x,y to the fluid grid, which is used to calculate the permeability κ(ϕ x,y ) (m2) assuming the Carman–Kozeny relation (Carman, Reference Carman1937; Phillips, Reference Phillips1991; Ghani et al. Reference Ghani, Koehn, Toussaint and Passchier2013) as

(1) $$\kappa \left( {{\phi _{x,y}}} \right) = {{{r^2}{{\left( {{\phi _{x,y}}} \right)}^3}} \over {45{{\left( {1 - {\phi _{x,y}}} \right)}^2}}}$$

Table 1. Parameters used in the simulations.

where r (m) is a fixed grain size. A hot fluid is thought to infiltrate the model from below at a higher temperature and pressure. If we assume that gravitational effects do not play a significant role on a scale of 5 m, then the fluid pressure evolution can be determined as

(2) $$\phi \beta \left[ {{{\partial P} \over {\partial t}}} \right] = \nabla \cdot \left[ {\left( {1 + \beta P} \right){\kappa \over \mu }\nabla {\rm{P}}} \right]$$

where ϕ is porosity, β is fluid compressibility (m2/N), P is fluid pressure (Pa), κ is permeability ( ${m^2})$ and μ is fluid viscosity (Pa s; Ghani et al. Reference Ghani, Koehn, Toussaint and Passchier2013). In the present simulations the porosity is constant and there is no associated volume change, since the reaction itself is not included. The fluid pressure gradients are then used to derive the fluid velocity v (m s–1) according to

(3) $$\vec v = - {{{{\rm{\kappa }} \over \mu }\nabla P} \over \phi }\;$$

The fluid velocity is then used in the general transport equation to derive the advective part for temperature and the different chemical species

(4) $${{\partial \left( {C,T} \right)} \over {\partial t}} + \vec v\;\nabla \left( {C,T} \right) - \;D\Delta \left( {C,T} \right) = 0$$

where the physical effects of advection and diffusion are separated and added after each time step as:

(5) $${\left( {C,T} \right)^t} = {\left( {C,T} \right)^{t - 1}} + \delta \left( {C,T} \right)_{adv}^t \,+ \,\delta \left( {C,T} \right)_{dif\,f}^t$$

where T is temperature (°C), C is concentration of species (mol kg–1) and D is the diffusion coefficient for temperature and matter diffusion respectively (m2 s–1). The IMEX (IMplicit+Explicit; Ascher et al. Reference Ascher, Ruuth and Spiteri1997) approach is used, where the advection is treated in an explicit way and the diffusion in an implicit way, with internal time loops in the advection to increase stability. The diffusion coefficient varies for temperature and the chemical species, and we assume that it is linearly related to the porosity. The model time loop is given in Figure 2 with an initial starting configuration and a given geometry, followed by an application of boundary conditions. The right- and left-hand side of the model are closed, while fluid is entering through the lower boundary and exits at the upper boundary. This information is used to calculate the fluid pressure gradients, before solving the transport equation for temperature and chemical species. In a final step, iphreeqc is used to calculate both fluid compositions and saturation indices.

Fig. 2. Schematic diagram illustrating the simulation model set-up. An initial structure with certain starting conditions builds the base of the system with a background porosity, temperature, fluid pressure and pore fluid with concentration of species and density. The model then applies boundary conditions fluid pressure (P f), temperature (T) and concentration of different species of the incoming fluid (c). The next step in the model is the evaluation of the fluid pressure evolution per time step followed by transport of temperature and concentration by advection and diffusion. The local concentration of species and temperature are then used in iphreeqc to calculate the local fluid properties including density and saturation indices of minerals.

2.b.2. Fluid composition: coupling with hydrochemical calculations

We added the iphreeqc class (Charlton & Parkhurst, Reference Charlton and Parkhurst2011) to the existing simulation framework in elle. The equilibrations with the dominant host-rock mineral are performed during the initialization step prior to the start of the simulations. The fluid compositions are then stored in sequence containers, the sizes of which are based on the species present in the pore- and infiltrating fluid.

The coupling between transport and hydrochemical calculations is performed in a sequential manner: we first solve for the transport by obtaining pressure (Equation (2)) and fluid velocities (Equation (3)). After determining the concentration changes at every particle, the fluid compositions stored at the particles are updated and predefined saturation states are reported. The calculation of a mineral saturation index SI p in iphreeqc is defined as (Parkhurst & Appelo, Reference Parkhurst and Appelo1999):

(6) $$S{I_P} = log\mathop \prod \limits_m^{{M_{aq}}} a_m^{{c_{m,p}}}$$

where M aq is the total number of aqueous master species, a m is the activity of the master species and c m,i is the stoichiometric coefficient of master species m. The saturation index is one that shows whether a fluid will tend to dissolve or precipitate a particular mineral, and is a dimensionless measure. Positive values indicate supersaturation, a value of zero corresponds to equilibrium, and negative values indicate undersaturation of the solution with respect to the master species. Because of uncertainties (e.g. in fluid temperature, ionic activity, analytical precision and thermodynamic data), a log SI of ± 0.2 is often considered to be in the range of equilibrium (Moldovanyi & Walter, Reference Moldovanyi and Walter1992). The function f p used by iphreeqc for phase equilibrium to reach the target SI p,target is:

(7) $${f_p} = \left( {ln{K_p} + \left[ {\ln \left( {10} \right)} \right]S{I_{p,{\rm{target}}}}} \right) - \mathop \sum \limits_m^{{M_{aq}}} {c_{m,p}}\ln \left( {{a_m}} \right),$$

where K p is the pure phase equilibrium (Parkhurst & Appelo, Reference Parkhurst and Appelo1999). For simplicity we define alkalinity as HCO3 and let iphreeqc calculate the pH based on alkalinity and inorganic carbon. By following this approach, we ensure that we always obtain the input pH of the solution for 100% of pore fluid or 100% of infiltrating fluid (Fig. 2). The computational efficiency is enhanced by parallelizing the hydrochemical calculations via multithreading using the software openMP (Chapman et al. Reference Chapman, Jost and Van der Pas2007).

2.c. Simulation set-up

2.c.1. Model geometry and boundary conditions

For the geometry of the model, we use a simple but geologically realistic scenario to avoid complex interferences and to study fluid mixing and developing saturation indices in detail (Fig. 3). In particular, we model a sedimentary sequence with layers of different chemical composition and permeability. In addition, one layer exhibits faults with distinctly lower permeabilities. In the model a vertical 2D simulation box is chosen with a width and height of 5 m. The side boundaries are confined for a solid and impermeable for the fluid. The upper boundary is open and can move as a function of gravity assuming a depth of 200 m, a density of 2700 g m–3 for the overlying rock and acceleration due to gravity of 9.81 m s–2. The fluid pressure at the upper boundary is hydrostatic and is kept constant. The lower boundary is confined for the solid so that the lowermost elements are not allowed to move vertically. The metal-rich fluid will enter along the whole lower boundary with a fixed concentration of species (see Table 2 for concentrations), fixed temperature (90 °C) and increase in fluid pressure (above hydrostatic) for the initial 2000 time steps until the lower boundary pressure is 3.3 MPa, 1.3 MPa above hydrostatic pressure. The 2000 model time steps are chosen such that the fluid pressure gradients within the model are low enough for numerical stability.

Fig. 3. Set-up of the simulations. (a) Mineralogy of the layers with two sandstone layers and one limestone layer in the centre. Initial boundary conditions are defined by equilibrated pore fluid in the different layers. Side walls are confined so that no fluid escapes sideways. The upper boundary is open and contains a low fluid pressure. The metal fluid enters through the lower boundary with an initially increasing and then constant fluid pressure, temperature and incoming fluid composition at the lower boundary. (b) Porosity/permeability scenario I with two permeable faults in the otherwise low-permeability lower layer. (c) Porosity/permeability scenario II with an additional low-permeability horizontally aligned area in the centre of the model above the two faults. (d) Flow lines for scenario II showing how flow concentrates in the two faults and then spreads out in the high-permeability layers and goes around the low-permeability zone in the centre.

Table 2. Fluid compositions used during simulations. Seawater composition is adapted from Ball & Nordstrom (Reference Ball and Nordstrom1991) and brine composition from Moldovanyi & Walter (Reference Moldovanyi and Walter1992) (Smackover Formation brine 55). Element concentrations are given as molarities (per kilogram seawater or kgw) as this is the standard concentration unit used by phreeqc. Alkalinity is defined as HCO3.

The model geometry has three horizontal layers: a lower sandstone layer overlain by a limestone and another sandstone layer (Fig. 3a). We test two scenarios to investigate the effect of local low- and high-permeability zones. For all scenarios the lower sandstone layer has a low porosity of 0.01 and permeability of 2.37 × 10−18 m2, whereas the other two layers have a high porosity of 0.10 and a higher permeability of 2.75 × 10−15 m2 (Fig. 3b). The low-permeability sandstone layer is transected by two high-permeability vertical faults with two different widths. The faults have a porosity of 0.10 and permeability of 2.75 × 10−15 m2. In scenario II an additional low-permeability zone is present in the middle limestone layer, defined by a low porosity of 0.01 and a permeability of 2.37 × 10−18 m2 (Fig. 3c, d) relative to the limestone layer (see also Table 2).

A total of 61 simulations were performed, with each simulation lasting 5 days on a normal desktop computer. Each time step includes the solution of the pressure diffusion equation, the advection and diffusion for temperature and 12 different species, and the following iphreeqc calculation. Up to 100 000 modelling steps were used for the simulations, with iphreeqc being solved every step for every node of the model (13 000 fluid nodes). The pore fluid is initially equilibrated with the respective layers, either sandstone or limestone.

2.c.2. Simulation specific information: fluid compositions

For the simulations we chose two simple but realistic fluid compositions: seawater (I) published by Ball & Nordstrom (Reference Ball and Nordstrom1991) and a representative example of a brine (II) from the Smackover formation in SW Arkansas (Moldovanyi & Walter, Reference Moldovanyi and Walter1992). The brine represents a highly saline formation water, and we chose well 55 as this exhibits the highest salt content (total dissolved solids: 336 492 mg L–1). The composition of these fluids is provided in Table 2. The error reported in Table 2 represents the charge balance error (CBE); a CBE of ± 5% is usually acceptable. We assume that the pore fluid is derived from seawater and is equilibrated with the dominant mineral species of the host formation, that is, it is a formation fluid. In our model we chose to equilibrate the seawater with calcite and quartz, representing formation water that is stagnating in the two chosen lithologies (limestone and sandstone, respectively) (Fig. 3). For the metal-bearing fluid we chose the composition of a saline formation water. Consequently, we need to account for the following 12 chemical species: Ba, Ca, Cl, Fe, K, Mg, Mn, Na, Pb, S, Si and Zn. We chose this scenario as mixing of formation versus fluxing fluid is frequently advocated to trigger mineralization in sedimentary basins (e.g. Hitchon & Friedman, Reference Hitchon and Friedman1969; Bjørlykke, Reference Bjørlykke1993). The initial mineral saturation states of the pore fluids and the infiltrating fluid calculated with the standard phreeqc database are shown in online Supplementary Figure S1 (available at http://journals.cambridge.org/geo).

2.d. Assumptions

For simplicity, the model contains several assumptions. We assume that we can approach the permeability–porosity relation by the Carman–Kozeny equation (Equation (1)) and that the solid can be approximated as a porous medium. The driving force for the fluid flux is assumed as a pressure gradient representing an incoming pulse of fluid that is at pressure between hydrostatic and lithostatic. The pressure gradient initially builds up in the model and, at later stages, flow is assumed to proceed along the pressure gradient. On the scale of the model the temperature diffusion is very fast (a factor of 1000 faster than matter diffusion; see also Bons et al. Reference Bons, Van Milligen and Blum2013) so that the hot incoming fluid increases the temperature in the whole box and does not lead to the development of density gradients resulting from temperature differences that are sustainable over the model time. Matter diffusion is a lot slower meaning that, on the scale of our simulation, density differences because of concentration differences of species could matter. We calculated the local density of every fluid node with iphreeqc and added it to the pressure evolution in several test runs; however, the effects of density differences were very low and did not change the pattern. In addition, our hot incoming fluid has a higher density than the formation water because of the metal content, so that flow into the model was only slowed down. We therefore argue that, on the scale of the model, density differences because of matter or temperature variations do not influence the flow and no instabilities develop. This could change once reactions are included in the model and material precipitates, changing the physical properties of the fluid (Lester et al. Reference Lester, Ord and Hobbs2012).

For both the chemical species and temperature diffusion we assume that the diffusion coefficient is a scalar that is linearly dependent on the porosity. In the current version of the code, the diffusion coefficients for temperature and matter for the same porosity vary by three orders of magnitude (Bons et al. Reference Bons, Van Milligen and Blum2013). The diffusion coefficients for the different chemical species do not vary, even though in reality they show a variation. In addition, diffusion is thought to be independent of the number of other species present in the fluid. The advection is calculated using the Lax–Wendroff scheme (Lax & Wendroff, Reference Lax and Wendroff1960) that adds the necessary numerical diffusive term to the advection to keep the calculation stable.

In iphreeqc we assume that the saturation index provides us with a proxy to understand which minerals will precipitate. In order to follow the effect of fluid mixing we track the saturation index of barite as an indicator mineral, as its saturation index is only high in areas where the fluids mix (Fig. 4). We sum up the saturation indices for different time steps to explore where in the model they remain high for a longer time period.

Fig. 4. Ideal mixing of pore fluid (seawater) and infiltrating fluid (oil-field brine) calculated using MIX in the iphreeqc module. pH, barite saturation state and CBE are obtained for three databases: phreeqc (solid lines), wateqf4 (dotted lines) and llnl (dashed lines). The blue lines represent the mixing with unaltered seawater and the orange lines show the resulting fluid if mixed with seawater equilibrated with calcite. Equilibration with quartz did not alter the results compared with the non-equilibrated case; only the experiments with non-equilibrated and calcite-equilibrated pore fluid are shown.

3. Results

3.a. Testing the robustness of the fluid mixing algorithm

We first performed pure mixing experiments using the iphreeqc module to assess the effect of gradual mixing on the evolution of pH and saturation states (Fig. 4). The simulations start with 100% pore fluid that is gradually mixed with the infiltrating fluid, adding 1% during each step of the simulation. The tests are performed with three different hydrochemical databases: phreeqc.dat, llnl.dat and wateqf4.dat.

Figure 4 shows that the maximum of the saturation index of barite is reached when 80–90% of pore fluid is still present followed by a slowly descending saturation index with increasing metal-rich fluid mixed with the pore fluid. Barite saturation exhibits the same trend for all three databases with only minor divergence between the respective values. The CBE of the phreeqc and wateqf4 databases are superimposed, while the CBE for the llnl database diverges above 10% of pore fluid. While the error obtained by using phreeqc and wateqf4 are above −4% the error, the llnl database that is used for the simulations is about −6%. The temperature and ionic strength evolution during ideal mixing of pore and infiltrating fluid is shown in online Supplementary Figure S2 (available at http://journals.cambridge.org/geo). In summary, the phreeqc and wateqf4 databases yield very similar results in terms of pH, saturation states, error and ionic strength. The trend of the barite saturation index is similar for all three databases. For the coupled simulations we chose the standard phreeqc database as it is expected that the CBE will remain in an acceptable regime, and general trends for barite saturation indices, the species investigated in this study, are similar for all tested databases.

3.b. Scenario I: fluid infiltration, fluid mixing and saturation indices within faults

Scenario I is designed to show the infiltration of the metal-rich fluid into the two faults within a layered sequence of variable chemistry and permeability (Fig. 5). Due to the high porosity and thus permeability of the faults, the fluid pressure is higher within the permeable faults leading to fluid infiltration into the low-permeability wall rock (Fig. 5a). The infiltrating fluid is visualized using Zn concentration that is only present in the incoming fluid, whereas the pore fluid is visualized by Mg concentration only present in the pore fluid. The metal-rich fluid infiltrates the faults with very steep concentration gradients at the fault walls (i.e. host-rock side of the interface between fault and wall rock). In contrast, there are more spread-out gradients in the infiltration front within the faults (Fig. 5b). While the metal-rich fluid infiltrates, the pore fluid is displaced from the faults. As with the metal fluid, the concentration gradients are very steep in the fault walls but more spread out in the actual faults (Fig. 5c).

Fig. 5. Evolution of fluid in two vertical faults with blue low and red high numbers; width of box is 5 m. (a) Evolution of the fluid pressure into the two faults. (b) Influx of metal-rich fluid into the faults with a Péclet number (see Section 4b for definition) of around 8000 for concentration in the fault. (c) Pore-fluid is displaced in the faults by the incoming fluid. (d) Evolution of barite saturation index in the faults. (e) Barite saturation index summed up over time to show how long the saturation index is high. (f) Temperature after 2 min of model run. Temperature is equilibrating as a result of fast diffusion (T-Péclet number of 100 at the start, and 8.3 after 11 days in the fault).

The barite saturation index (Fig. 5d) is depicted as a snapshot in time, showing high indices at the lower boundary (this is a model–boundary artefact), in the fault walls and in the faults themselves where the two fluids meet. The positive index is larger within the fault than in the fault walls. When considering the accumulated index over time (Fig. 5e), the positive index is very sharp and small at the lower boundary and the fault walls but smeared out and relatively wide within the faults. The area with the highest index is found within the fault walls and in front of the infiltrating metal-rich fluid in the faults, because only 10–20% of incoming fluid is needed to reach the highest index. Finally, the temperature advection is shown for a very small time period of only 2 minutes relative to the 11 days for the other values (Fig. 5f). The temperature field is already smeared out after this short time, with diffusion being relatively fast with respect to the advection. After 50 days the whole box has an almost equal temperature.

To investigate the mixing behaviour of the fluid we produce horizontal and vertical profiles through one of the faults (Fig. 6). A horizontal profile through the wider fault is shown in Figure 6a, b with the concentration of Zn for the metal fluid, Mg for the pore fluid and the saturation index of barite. The metal fluid is displacing almost all the pore fluid in the fault. Mixing is concentrated on the fault walls where the metal-rich fluid stops infiltrating. Here the barite saturation index is highest in a zone of c. 10 cm width into the wall rock. In the fault itself the index is relatively low where the pore fluid is displaced. The high saturation index of barite in the fault wall geometrically resembles a relative stationary wave. A vertical cross-section through the fault is shown in Figure 6c from bottom (left) to top (right). Here the mixing is more spread out with a steep advection front where the metal-rich fluid is displacing, followed by a dragged-out area where 10–20% of the pore fluid is still present. The barite saturation index is more complicated with a frontal part that increases up to 1 m before the actual metal fluid displaces the pore fluid. This frontal zone is followed by a zone that is very steep (c. 20 cm wide) and with a very high index (c. 10–70%) of metal fluid followed by a descending index behind the advection front. This zone of high saturation index of barite resembles a geometrical wave that travels through the fault. The evolution in time is shown in Figure 6d, e for two time steps in a vertical cut through the fault, where diffusion leads to a rounding of the concentration pattern and a more gradual change into the fault walls. This is reflected by the saturation index wave that progresses into the fault walls, growing in wavelength from c. 20 to 40 cm over time. Figure 6f shows the time evolution of the saturation index for one single point in the fault with a very steep increase when the advection front passes, followed by a fast decay to about half the maximum index and a slow increase towards the end of the model run.

Fig. 6. Evolution of pore fluid, metal fluid and barite saturation index in one fault. (a, b) Cross-section through a fault where the metal fluid has infiltrated: (a) fluids and barite solution index in three plots and (b) all confined to one plot; mixing best in the fault walls. (c) Vertical section through the incoming metal fluid showing how the pore fluid is pushed away and the barite saturation index is highest within the area where 80% of pore fluid is still present. This wave of high barite saturation index is moving with the influxing fluid through the fault. (d) Two evolutions for the incoming fluid where the concentration curve flattens out into the wall rock as a result of diffusion. (e) Saturation index of barite for (d), showing how the maxima move with diffusion into the wall rock of the fault. (f) Time evolution of a point in the model with a very fast increase in barite saturation index when the fluids mix, followed by a relaxation when the metal fluid dominates (full model run is 231 days).

The time evolution of the cumulative saturation index of barite (Fig. 7) shows how the metal fluid enters the faults with a stationary wave of high cumulative index developing in the fault walls (and at the lower boundary of the model) and a travelling wave with a larger wavelength within the fault. The zone within the fault is moving, is slowing down with time and is spreading out progressively.

Fig. 7. Evolution of summed saturation index of barite in the faults: (a–l) progressive time steps. The saturation index shows a relatively stationary and diffusion-controlled wave of high accumulated SI of barite in the wall rock, but a travelling wave of high SI of barite in the faults that grow in width.

3.c. Scenarios I and II: effects of porous layers and seals

In this section we compare scenarios I and II to explore the patterns that develop when the fluid leaves the faults and enters the porous limestone layer, as well as how an additional low-porosity zone within the layer influences the patterns. Once the metal-rich fluid enters the porous limestone layer it spreads sideways and builds a large plume-like structure (Fig. 8a, b). In this case, the wider fault leads to further infiltration of the fluid than the smaller fault. Above the faults in the model, the concentration pattern of Zn shows relatively wide gradients into the porous layer, whereas it is relatively steep in the fault walls. The accumulated saturation index of barite shows very high saturation in the fault walls. In addition, zones of high accumulated saturation indices develop at the sides of the faults where the metal-rich fluid meets the pore fluid in the porous layer. In contrast, the saturation indices are low within and above the large fault and high above the smaller fault.

Fig. 8. Effects of horizontal permeable layer and a horizontal low-permeability layer on flux and saturation index of barite. (a) Accumulated barite saturation index in a simulation with a permeable layer on top of the two faults but no seal. High saturation index above the fault is present on the fault sides and where fluids between the two faults mix as a result of diffusion. (b) Concentration of incoming metal fluid for scenario I. (c) Accumulated barite saturation index in a simulation with a horizontal permeable layer covered by a central low-permeability layer (scenario II). (d) Concentration of metal fluid in scenario II. Fluid is pushed against the barrier and flows around it at the left- and right-hand side. Below the barrier the accumulated saturation index of barite is high, indicating that barite would precipitate at this location. All pictures shown are after 115 days.

Once a low-porosity zone is added to the limestone layer it has a profound effect on the infiltration of the metal-rich fluid and on the accumulated saturation index of barite (Fig. 8c, d). Above the large fault, the fluid flows against the low-porosity zone and concentration gradients become steep, similar to those in the fault walls. This is indicated by the accumulated saturation index of barite, which is becoming high below the seal, again similar to what is happening in the fault walls.

The evolution of the accumulated saturation index of barite over time is shown in Figure 9 from when the metal-rich fluid leaves the fault. The saturation index moves sideways and builds a table-like structure below the seal. In addition, the high indices at the entry points of the faults to the permeable layer merge between the faults.

Fig. 9. Time evolution of the accumulated saturation index of barite in simulations with a permeable layer and a barrier (scenario II). The highest indices are in the fault walls, between and at the exit of faults, in the permeable layer and below the barrier.

4. Discussion

4.a. Stationary and travelling waves in fault networks

It is well known that dissipative structures will develop in reaction–diffusion and reaction–diffusion–advection systems, forming structures in space and time (Glansdorff & Prigogine, 1971; Budroni & De Wit, Reference Budroni and De Wit2017), and are important during mixing processes in ore deposits (Lester et al. Reference Lester, Ord and Hobbs2012). One form of these structures are waves and spots of reactants that are stationary (e.g. Turing patterns) or travelling (Belousov, Reference Belousov1959; Zaikin & Zhabotinsky, Reference Zaikin and Zhabotinsky1970). As discussed by Lester et al. (Reference Lester, Ord and Hobbs2012), an additional form of enhanced mixing in ore deposits can be induced by geometrical structures.

In accordance with nomenclature from chemistry and physics for chemical waves, and because these zones represent waves geometrically, we use the term ‘travelling waves’ for zones of maximum instantaneous saturation index that move in space and time, and ‘stationary waves’ for zones of high saturation index that do not move (Faraday, 1828; Belousov, Reference Belousov1959; Zaikin & Zhabotinsky, Reference Zaikin and Zhabotinsky1970) (Fig. 6a–c). The highest saturation indices develop in front of the incoming fluid with 80–90% pore fluid still present. Once the incoming fluid dominates, the saturation index goes down again. This leads to moving waves of high saturation indices in areas of high fluid velocity where advection dominates. In areas where diffusion is dominant, the mixing zone moves very slowly and the waves of high saturation indices are almost stationary within the model time. Moving waves are most pronounced in high-permeability zones, so that faults as well as porous layers act as advection dominating systems. These advection-dominated systems only lead to a minor amount of fluid mixing, where fluids do not travel towards each other to mix but rather displace each other. As long as these travelling waves of saturation indices move faster than a reaction can progress, they will be dissipative and disappear. What leads to the main mixing of fluids on the local scale is diffusion of species, and this is mainly happening at the interfaces between the fluids. There are typically two main interfaces in the system: (1) between the fluids in high-permeability zones, where the metal fluids displace the pore fluid; and (2) at low-permeability zones, where advection is not dominating. High-permeability zones are the faults and the permeable layer, whereas low-permeability zones are represented by the fault walls (Figs 57) and low-permeability barriers (Figs 8, 9).

In considering mineralization, we assume that a positive saturation index leads to precipitation; however, this assumption has a time factor involved, such that the index needs to be positive for a given time or else minerals will not precipitate. We therefore also summed up the saturation index over time, which gives an indication of locations where the index is positive over a long time and should lead to precipitation. These zones of high saturation index travel as waves and typically develop where the velocity of the incoming fluid is high, whereas they are relatively stationary where the velocity is low. Travelling waves are therefore present in the permeable faults where we do not expect precipitation of minerals (e.g. barite) to develop because of mixing. The high advection in the faults potentially leads to them remaining permeable. In contrast, where the velocity is low (in the fault walls, below low-permeability zones and also between the two faults) the waves of high saturation index are stable and minerals such as barite have enough time to precipitate. Feeders like the faults would remain open in a mineral system, whereas areas of low fluid velocity represent sinks and can lead to mineralization.

4.b. Péclet and Damköhler numbers

In a system where advection and diffusion are important, the dimensionless Péclet number (Pe) is used to describe the relationship between advection rate and diffusion rate for chemical and temperature transport:

(8) $$Pe = {{vL} \over D}$$

where v is the fluid velocity, L the characteristic length scale of the system and D the diffusion coefficient. Advection dominates for high values of Pe, whereas diffusion takes over for low values of Pe. Lester et al. (Reference Lester, Ord and Hobbs2012) suggested different regimes in hydrothermal deposits that can be distinguished based on the ratio of diffusion and advection/convection, that is, the Péclet number. In addition, in any given system the local Péclet number may vary significantly both spatially and temporally (Ortoleva et al. Reference Ortoleva, Chadam, Merino and Sen1987; Bons et al. Reference Bons, Van Milligen and Blum2013; Koehn et al. Reference Koehn, Piazolo, Beaudoin, Kelka, Spruženiece, Putnis and Toussaint2021; Fig. 10). Permeable systems in our model such as the faults have a high Péclet number and are advection dominated, whereas fault walls and areas below the seal have a low Péclet number and are diffusion dominated. In a system where reactions occur along with advection and diffusion, two additional dimensionless numbers are used to assess the influence of the relative rates of these processes. These two numbers are: (1) Damköhler I for reaction rate relative to advection rate

(9) $$D{a_I} = {R \over v},$$

Fig. 10. (a) Initial evolution of the Péclet number (Pe) with time in the simulations. Initially, the Péclet number for concentration of species is decaying fast. Inset shows a log-log plot of Péclet numbers for concentration versus temperature, illustrating that temperature is a lot more diffusion controlled than concentration. (b) Schematic illustration of a high-permeability fluid channel or fault with a travelling wave of high saturation index of barite. Stable waves develop in the wall rock of the fault where mineralization can occur. A plot of SI barite along the fault on the left-hand side of the figure shows a snapshot of the travelling wave of high-SI barite, whereas the plot at the bottom of the figure shows the stationary waves of high-SI barite in the wall rock of the faults. The index is only high where the mixing between the fluids is optimal, with c. 80–90% of pore fluid and only a small amount of incoming fluid of 10–20%. The index therefore decays once the incoming fluid dominates and displaces the pore fluid. The geometries are similar to the experiment shown in Figure 1c, even though the experiment did not directly involve fluid mixing.

and (2) Damköhler II, defined as reaction rate R relative to diffusion rate

(10) $$D{a_{II}} = {{RL} \over D}$$

Changes in the Damköhler numbers lead to changes in the shape and spacing of wormholes in karst systems (Szymczak & Ladd, Reference Szymczak and Ladd2009); the Damköhler numbers also influence reactive transport for the alteration of pore space (Mostaghimi et al. Reference Mostaghimi, Liu and Arns2016), as well as the roughness of reaction fronts (Koehn et al. Reference Koehn, Piazolo, Beaudoin, Kelka, Spruženiece, Putnis and Toussaint2021).

In terms of the Péclet number in our simulations, Figure 10 illustrates that it is decaying rapidly at the beginning of the simulations where it is initially very high. The inset in Figure 10 also shows the extreme difference between temperature and concentration Péclet numbers, because of the difference in diffusion coefficients (three orders of magnitude; Bons et al. Reference Bons, Van Milligen and Blum2013). The concentration remains within the permeable zones in the advection-dominated regime as for the faults and permeable layers, whereas the temperature reaches a diffusion-controlled regime early in the simulations (Pe of 10 after 100 hours). The Damköhler I number that relates fluid velocity to reaction rate in the system is also evolving, since the Péclet number changes. If we assume a reaction rate of roughly 10−10 m s–1 for barite growth (Kuwahara et al. Reference Kuwahara, Liu, Makio and Otsuka2016) the Damköhler I number is c. 10−4 with a fluid velocity of c. 10−6 m s–1. This is a relatively low number and means that, within the faults where the wave with high saturation index is travelling, the reaction rate is probably too low to lead to precipitation. However, the Damköhler number II is c. 1–10 for diffusion of chemical species, 1 for diffusion in the wall rock and low-permeability rocks, and 10 within the high-permeability faults and layers. The reaction and diffusion rate are then on the same length and time scales so that minerals can probably precipitate once diffusion is dominant, especially within the low-permeability wall rock of faults and low-permeability zones in general.

4.c. Model development: a future perspective

We have successfully merged a simple transport model that deals with two mixing fluids with iphreeqc to study the development of saturation indices advecting and diffusing 12 different species and temperature. We note that the utilization of different databases might yield different results (Fig. 4, Figure S2). Uncertainties must therefore be considered for the results presented here that arise from the underlying uncertainties related to the empirically determined parameters present in the databases.

The next steps in the development of the models should include the actual reactions and potentially also changes in the porosity structure of the solid, as well as more realistic diffusion models, hydromechanical feedbacks and the creation of fracture networks. Reactions would of course change the local concentration in the fluid (an additional term in Equation (5)) and therefore influence the time evolution of the saturation indices of minerals. This already leads to the question of how the metal-rich fluid would initially behave in the faults, since it is not in equilibrium with the solid. It would therefore lead directly to the precipitation of minerals from the oversaturated solution, even without mixing. Depending on the volume change and the possible dissolution of material versus precipitation, this could change the permeability of the fault and potentially clog it. It would be interesting to explore the effects of direct reactions in a new fluid relative to effects of mixing of the two fluids. Under what circumstances will the system remain open to incoming fluids and when will it close? In terms of the diffusion, the different diffusion coefficients of the 12 species used as well as the effect of overall concentration on the coefficients also needs to be explored.

However, care must be taken in these approaches since the variability of parameters becomes increasingly complex while the potential error grows, meaning that it will be challenging to separate real physical effects from numerical model artefacts. This is the main advocate for simplicity within the current model.

5. Conclusion

In this contribution we have introduced a reactive transport environment that is open-source and may include the mechanics of a porous media via a lattice-spring model. In the current study we focus on the coupling of the multicomponent fluid transport with simple hydrochemical calculations. We show that, depending on the dominant transport process (i.e. advection or diffusion), supersaturated fronts will localize at permeability contrasts. We applied our framework to a simple fluid-mixing scenario to demonstrate several mechanisms that are active at the spatial and temporal scale of the investigated models (outcrop scale, 5 × 5 m): (1) infiltrating fluid displaces the pore fluid; (2) heat diffusion is fast enough to quickly reach temperature equilibrium, and the temperature dissipates fast; (3) local diffusion coupled with advection along faults/fractures is the dominant mixing process, and permeability differences are important; (4) stationary waves of high saturation indices lead to mineralization in quiet zones; (5) travelling waves of high saturation indices can leave faults or fractures permeable; and (6) a high density of fault or fracture networks increases the surface area where fluids meet, leading to enhanced mixing and mineralization.

These mechanisms are paramount in the context of hydrothermal mineral systems. The displacement of pore fluid is most pronounced in high-permeability regions such as faults and fractures, thus preventing efficient mixing and subsequent precipitation in these zones. The travelling waves of high saturation indices that are observable within the fault zones are advection dominated and do not lead to considerable mineral precipitation because of the timescale between reaction rate and transport. Stationary waves of high saturation indices prevail in the wall rock of faults, below low-permeability zones and between permeable faults. These stationary waves will lead to the precipitation of minerals as they retain high saturation states over a long time.

In summary, ore deposits that develop as a result of fluid mixing may be hosted in zones that are diffusion dominated with low fluid velocities, whereas high-permeability feeders such as faults and fractures may stay open even though waves of high saturation indices pass through them. The best mixing scenario would encounter a fracture or fault network that has about twice the spacing of the diffusion distance of 10–20 cm of the system, leading potentially to complete mineralization of the host rock. While the simulations presented here are simple, comprising only a pore fluid equilibrated with the host rock and an infiltrating metal-bearing fluid, the results allow us to gain new insights into the processes that lead to localized reactions (including mineral precipitation) that are active during the formation of economic hydrothermal mineralization on the outcrop scale, and demonstrate that permeability variations such as fracture networks and faults significantly enhance mixing and mineralization.

Supplementary material

To view supplementary material for this article, please visit https://doi.org/10.1017/S001675682200070X

Acknowledgements

This research was part of the project “Integrated geophysical-structural-kinematic analysis of the fault patterns in Northern Bavaria” funded by the Bayerisches Landesamt für Umwelt (LfU) and the Bavarian Ministry of Science and Art (StMWK) of the Geothermal Alliance Bavaria (GAB). GM, AB and DK acknowledge funding by the Natural Environment Research Council (Case Studentship grant no. NE/R007772/1, title: “Predicting ore mineralization using quantitative forward modelling and high resolution analytical data”). CS acknowledges MinEx CRC phase 1 for funding the acquisition of microXRF maps. RT acknowledges the support of the CNRS INSU CESSUR programmes and the Research Council of Norway through its Centre of Excellence funding scheme (project no. 262644). The authors also acknowledge the support of ITN FlowTrans, the European Union’s Seventh Framework Programme for research (grant no. 316889). We thank reviewers and the associated editor for very useful comments that enhanced the earlier version of the manuscript.

References

Ascher, UM, Ruuth, SJ and Spiteri, RJ (1997) Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics 25, 151–67.CrossRefGoogle Scholar
Ball, W and Nordstrom, DK (1991) User’s manual for WATEQ4F, with revised thermodynamic data base and text cases for calculating speciation of major, trace, and redox elements in natural waters. USGS Open-File Report 91–183, 193.Google Scholar
Belousov, BP (1959) Periodicheski deistvuyushchaya reaktsia i ee mekhanism [Periodically acting reaction and its mechanism]. Sbornik referatov po radiotsionnoi meditsine, 1958 [Collection of abstracts on radiation medicine, 1958]. Moscow: Medgiz, 145–47.Google Scholar
Bjørlykke, K (1993) Fluid flow in sedimentary basins. Sedimentary Geology 86, 137–58.CrossRefGoogle Scholar
Bons, PD, Fusswinkel, T, Gomez-Rivas, E, Markl, G, Wagner, T and Walter, B (2014) Fluid mixing from below in unconformity-related hydrothermal ore deposits. Geology 42, 1035–8. doi: 10.1130/G35708.1.CrossRefGoogle Scholar
Bons, PD, Koehn, D and Jessell, MW (2008) Microdynamics Simulation. New York: Springer. Lecture Series in Earth Sciences, p. 406.CrossRefGoogle Scholar
Bons, PD, Van Milligen, BP and Blum, P (2013) A general unified expression for solute and heat dispersion in homogeneous porous media. Water Resources Research 49. doi: 10.1002/wrcr.20488.CrossRefGoogle Scholar
Budroni, MA and De Wit, A (2017) Dissipative structures: from reaction-diffusion to chemo-hydrodynamic patterns. Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 104617.CrossRefGoogle ScholarPubMed
Carman, PC (1937) Fluid flow through granular beds. Transactions of the Institution of Chemical Engineers 15, 150.Google Scholar
Chapman, B, Jost, G and Van der Pas, R (2007) Using OpenMP: Portable Shared Memory Parallel Programming. Cambridge, MA: MIT Press.Google Scholar
Charlton, SR and Parkhurst, DL (2011) Modules based on the geochemical model PHREEQC for use in scripting and programming languages. Computers & Geosciences 37, 1653–63.CrossRefGoogle Scholar
Chen, L, Kang, Q, Carey, B and Tao, W-Q (2014) Pore-scale study of diffusion-reaction processes involving dissolution and precipitation using the lattice Boltzmann method. International Journal of Heat and Mass Transfer 75, 483–96.CrossRefGoogle Scholar
Chen, L, Kang, Q, Robinson, BA, He, Y-L and Tao, W-Q (2013) Pore-scale modeling of multiphase reactive transport with phase transitions and dissolution-precipitation processes in closed systems. Physical Review E 87, 043306.CrossRefGoogle ScholarPubMed
De Lucia, M, Kühn, M, Lindemann, A, Lübke, M and Schnor, B (2021) POET (v0.1): speedup of many-core parallel reactive transport simulations with fast DHT lookups. Geoscientific Model Development 14, 7391–409.CrossRefGoogle Scholar
Fan, HR, Hu, FF, Yang, KF and Wang, KY (2006) Fluid unmixing/immiscibility as an ore-forming process in the giant REE–Nb–Fe deposit, Inner Mongolian, China: evidence from fluid inclusions. Journal of Geochemical Exploration 89, 104–7.CrossRefGoogle Scholar
Faraday, M (1827) Chemical Manipulation, Being Instructions to Students in Chemistry. London: W. Phillips, 656 p.Google Scholar
Fredrich, JT, Digiovanni, AA and Noble, DR (2006) Predicting macroscopic transport properties using microscopic image data. Journal of Geophysical Research-Solid Earth 111. doi:10.1029/2005JB003774.CrossRefGoogle Scholar
Fu, B, Williams, PJ, Oliver, NH, Dong, G, Pollard, PJ and Mark, GM (2003) Fluid mixing versus unmixing as an ore-forming process in the Cloncurry Fe-oxide-CuAu District, NW Queensland, Australia: evidence from fluid inclusions. Journal of Geochemical Exploration 78, 617–22.CrossRefGoogle Scholar
Fusswinkel, T, Wagner, T, Wälle, M, Wenzel, T, Heinrich, CA and Markl, G (2013) Fluid mixing forms basement-hosted Pb-Zn deposits: insight from metal and halogen geochemistry of individual fluid inclusions. Geology 41, 679–82.CrossRefGoogle Scholar
Gessner, K, Blenkinsop, TG and Sorjonen-Ward, P eds. (2018) Characterization of Ore-Forming Systems from Geological, Geochemical and Geophysical Studies. London: Geological Society of London.Google Scholar
Gessner, K, Kühn, M, Rath, V, Kosack, C, Blumenthal, M and Clauser, C (2009) Coupled process models as a tool for analysing hydrothermal systems. Surveys in Geophysics 30, 133–62.CrossRefGoogle Scholar
Ghani, I, Koehn, D, Toussaint, R and Passchier, CW (2013) Dynamic development of hydrofracture. Pure and Applied Geophysics 170, 1685–703.CrossRefGoogle Scholar
Ghani, I, Koehn, D, Toussaint, R and Passchier, C (2015) Dynamics of hydrofracturing and permeability evolution in layered reservoirs. Frontiers in Physics 3, 67.CrossRefGoogle Scholar
Glansdorff & Prigogine (1971) Thermodynamic Theory of Structure, Stability and Fluctuations. Hoboken: Wiley-Interscience, 306p.Google Scholar
Haynes, DW, Cross, KC, Bills, RT and Reed, MH (1995) Olympic Dam ore genesis; a fluid-mixing model. Economic Geology, 90, 281307.CrossRefGoogle Scholar
Hitchon, B and Friedman, I (1969) Geochemistry and origin of formation waters in the western Canada sedimentary basin—I. Stable isotopes of hydrogen and oxygen. Geochimica et Cosmochimica Acta 33, 1321–49.CrossRefGoogle Scholar
Hobbs, BE and Ord, A (2018) Episodic modes of operation in hydrothermal gold systems: Part II. A model for gold deposition. In Characterization of Ore-forming Systems from Geological, Geochemical and Geophysical Studies (eds Gessner, K, Blenkinsop, TG and Sorjonen-Ward, P), pp. 147–64. Geological Society of London, Special Publication no. 453.Google Scholar
Hollinsworth, AD, Koehn, D and Dempster, AK (2019) Structural controls on the interaction between basin fluids and a rift flank fault: constraints from the Bwamba Fault, East African Rift. Journal of Structural Geology 118, 236–49.CrossRefGoogle Scholar
Hrennikoff, A (1941) Solution of problems of elasticity by the framework method. Journal of Applied Mechanics 8, A619715.CrossRefGoogle Scholar
Jamtveit, B, Putnis, CV and Malthe-Sørenssen, A (2009) Reaction induced fracturing during replacement processes. Contributions to Mineralogy and Petrology 157, 127–33.CrossRefGoogle Scholar
Kang, Q, Chen, L, Valocci, AJ and Viswanathan, HS (2014) Pore-scale study of dissolution-induced changes in permeability and porosity of porous media. Journal of Hydrology 517, 1049–55.CrossRefGoogle Scholar
Koehn, D, Piazolo, S, Beaudoin, NE, Kelka, U, Spruženiece, L, Putnis, CV and Toussaint, R (2021) Relative rates of fluid advection, elemental diffusion and replacement govern reaction front patterns. Earth and Planetary Science Letters 565, 116950.CrossRefGoogle Scholar
Koehn, D, Piazolo, S, Sachau, T and Toussaint, R (2020) Fracturing and porosity channelling in fluid overpressure zones in the shallow earth’s crust. Geofluids, article no. 7621759. doi: 10.1155/2020/7621759.CrossRefGoogle Scholar
Korrani, AKN, Sepehrnoori, K and Delshad, M (2015) Coupling IPhreeqc with UTCHEM to model reactive flow and transport. Computers & Geosciences 82, 152–69.CrossRefGoogle Scholar
Kuwahara, Y, Liu, W, Makio, M and Otsuka, K (2016) In situ AFM study of crystal growth on a barite (001) surface in BaSO4 solutions at 30°C. Minerals 6, 117.CrossRefGoogle Scholar
Lax, PD and Wendroff, B (1960) Systems of conservation laws. Communications on Pure and Applied Mathematics 13, 217–37.CrossRefGoogle Scholar
Leach, DL, Taylor, RD, Fey, DL, Diehl, SF and Saltus, RW (2010) A deposit model for Mississippi Valley-type lead-zinc ores. In Mineral Deposit Models for Resource Assessment. Reston, VA: United States Geological Survey. Scientific Investigations Report.CrossRefGoogle Scholar
Lester, DR, Ord, A and Hobbs, BE (2012) The mechanics of hydrothermal systems: II. Fluid mixing and chemical reactions. Ore Geology Reviews 49, 4571.CrossRefGoogle Scholar
Liu, W, Spinks, S, Glenn, M, Macrai, C and Pearce, M (2021) How carbonate dissolution facilitates sediment-hosted Zn-Pb mineralisation. Geology 49, 1363–8.CrossRefGoogle Scholar
Loukola-Ruskeeniemi, K (1999) Origin of black shales and the serpentinite-associated Cu-Zn-Co ores at Outokumpu, Finland. Economic Geology 94, 1007–28.CrossRefGoogle Scholar
Manwart, C, Aaltosalmi, U, Koponen, A, Hilfer, R and Timonen, J (2002) Lattice-Boltzmann and finite-difference simulations for the permeability for three-dimensional porous media. Physical Review E 66, 016702.CrossRefGoogle ScholarPubMed
Moldovanyi, EP and Walter, LM (1992) Regional trends in water chemistry, Smackover Formation, southwest Arkansas: geochemical and physical controls. AAPG Bulletin 76, 864–94.Google Scholar
Moosavi, R, Kumar, A, De Wit, A and Schroter, M (2019) Influence of mineralization and injection flow rate on flow patterns in three-dimensional porous media. Physical Chemistry Chemical Physics 21, 14605–11.CrossRefGoogle ScholarPubMed
Mostaghimi, P, Liu, M and Arns, CH (2016) Numerical simulation of reactive transport on micro-CT images. Mathematical Geosciences 48, 963–83.CrossRefGoogle Scholar
Nardi, A, Idiart, A, Trinchero, P, de Vries, LM and Molinero, J (2014) Interface COMSOL-PHREEQC (iCP), an efficient numerical framework for the solution of coupled multiphysics and geochemistry. Computers & Geosciences 69, 1021.CrossRefGoogle Scholar
Neumayr, P, Walshe, J, Hagemann, S, Petersen, K, Roache, A, Frikken, P, Horn, L and Halley, S (2008) Oxidized and reduced mineral assemblages in greenstone belt rocks of the St Ives gold camp, Western Australia: vectors to high-grade ore bodies in Archaean gold deposits? Mineralium Deposita 43, 363–71.CrossRefGoogle Scholar
Ord, A, Hobbs, B, Walshe, J, Metcalfe, G, Trefry, M, Zhao, CB, Lester, D, Regenaur-Lieb, K and Rudman, M (2010) Architecture and mixing for crustal plumbing systems: St Ives as a type example. In Proceedings of the Tenth Biennial SGA Meeting of the Society for Geology Applied to Mineral Deposits (eds, Schmidt, A, Brugger, J), pp. 833–35. Geneva: Society for Geology Applied to Mineral Deposits.Google Scholar
Ortoleva, P, Chadam, J, Merino, E and Sen, Z (1987) Geochemical self-organization ii: the reactive-infiltration instability. American Journal of Science 287, 1008.CrossRefGoogle Scholar
Parkhurst, DL and Appelo, C (1999) User’s guide to PHREEQC (Version 2). Water-Resources Investigations Report 99, 4259.Google Scholar
Phillips, OM (1991) Flow and Reactions in Permeable Rocks. Cambridge: Cambridge University Press.Google Scholar
Piazolo, S, Bons, PD, Griera, A, Lloens, M-G, Gomez-Rivas, E, Koehn, D, Wheeler, J, Gardner, R, Godinho, JRA, Evans, L, Lebensohn, RA and Jessell, MW (2019) A review of numerical modelling of the dynamics of microstructural development in rocks and ice: past, present and future. Journal of Structural Geology 125, 111–23.CrossRefGoogle Scholar
Press, WH, Flannery, BP, Teukolsky, SA and Vetterling, WT (1992) Numerical Recipes in C. The Art of Scientific Computing. Cambridge: Cambridge University Press, 1020 p.Google Scholar
Saripalli, KP, Meyer, PD, Bacon, DH and Freedman, VL (2001) Changes in hydrologic media due to chemical reactions: a review critical reviews. Environmental Science & Technology 31, 311–49.Google Scholar
Schwinn, G, Wagner, T, Baatartsogt, B and Markl, G (2006) Quantification of mixing processes in ore-forming hydrothermal systems by combination of stable isotope and fluid inclusion analyses. Geochimica et Cosmochimica Acta 70, 965–82.CrossRefGoogle Scholar
Shabani, A, Kalantariasl, A, Abbasi, S, Shahrabadi, A and Aghaei, H (2019) A coupled geochemical and fluid flow model to simulate permeability decline resulting from scale formation in porous media. Applied Geochemistry 107, 131–41.CrossRefGoogle Scholar
Shabro, V, Torres-Verdin, C, Javadpour, F and Sepehrnoori, K (2012) Finite-difference approximation for fluid flow simulation and calculation of permeability in porous media. Transport in Porous Media 94, 775–93.CrossRefGoogle Scholar
Steefel, CI, Appelo, CAJ, Arora, B, Jacques, D, Kalbacher, T, Kolditz, O, Lagneau, V, Lichtner, PC, Mayer, KU, Meeussen, JCL and Molins, S (2015) Reactive transport codes for subsurface environmental simulation. Computational Geosciences 19, 445–78.CrossRefGoogle Scholar
Steefel, CI and Lasaga, AC (1994) A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems. American Journal of Science 294, 529–92. doi: 10.2475/ajs.294.5.529.CrossRefGoogle Scholar
Szymczak, P and Ladd, AJC (2009) Wormhole formation in dissolving fractures. Journal of Geophysical Research 114. doi: 10.1029/2008JB006122.CrossRefGoogle Scholar
Tartakovsky, AM and Meakin, P (2006) Pore scale modeling of immiscible and miscible smoothed particle hydrodynamics. Advances in Water Resources 29, 1464–78.CrossRefGoogle Scholar
Wagner, T, Okrusch, M, Weyer, S, Lorenz, J, Lahaye, Y, Taubald, H and Schmitt, RT (2010) The role of the Kupferschiefer in the formation of hydrothermal base metal mineralization in the Spessart ore district, Germany: insight from detailed sulfur isotope studies. Mineralium Deposita 45, 217–39.CrossRefGoogle Scholar
Wei, Y and Cao, X (2021) A COMSOL-PHREEQC coupled python framework for reactive transport modeling in soil and groundwater. Groundwater 60(2), 284–94. doi: 10.1111/gwat.13144.CrossRefGoogle ScholarPubMed
Weihua, L, Spinks, SC, Glenn, M, MaCrae, C and Pearce, MA (2021) How carbonate dissolution facilitates sediment-hosted Zn-Pb mineralization. Geology 49, 1363–8.Google Scholar
Wilkinson, J (2001) Fluid inclusions in hydrothermal ore deposits, Lithos 55, 229–72.CrossRefGoogle Scholar
Zaikin, AN and Zhabotinsky, AM (1970) Concentration wave propagation in two-dimensional liquid-phase self-oscillating system. Nature 225, 535–7.CrossRefGoogle ScholarPubMed
Zhao, C, Hobbs, BE, Hornby, P, Ord, A, Peng, S and Liu, L (2007) Theoretical and numerical analyses of chemical-dissolution front instability in fluid-saturated porous rocks. Numerical and Analytical Methods in Geomechanics 32, 1107–30.CrossRefGoogle Scholar
Figure 0

Fig. 1. Example of mineralization at different scales within and adjacent to high-permeability sites. (a) Mineralization associated with faulting, Trollers Gill Barite–Fluorspar–Galena mine, Yorkshire Dales, UK. (b) MicroXRF geochemical mapping image acquired with a Bruker M4 Tornado from the hole NDIBK04 drilled by the National Drilling Initiative through the MinEx CRC. NDIBK04 is the most prospective well in the East Tennant Creek (Northern Territory) NDI region, an emerging mineral province for base metal mineral exploration (Zn, Cu and As). This sample is from the top of a c. 150-m-thick banded graphitic schist of the Palaeoproterozoic Alroy Formation intersected at depth 236.2–237 m. The nearest geochronology available is from a sample at 301–307 m depth, giving a maximum depositional age of 1969 ± 10 Ma (dated by SHRIMP at Geoscience Australia with zircons). The main vein has gypsum at its centre with Fe–Mn-rich carbonate to the edge. Such vein geometries can arise from the interaction with different fluids or by the evolution of a single infiltrating fluid with the host rock. The presence of gypsum in the centre of the vein points towards an open channel that prevailed after the diffusion-dominated reaction between fluid and rock that led to the development of the reaction rims surrounding the veins. (c) Experiment of calcite precipitation in low-porosity dolomite directly adjacent to high-porosity fracture in response to carbonate-rich fluid infiltration (courtesy of R. Morgan and Q. Fisher).

Figure 1

Table 1. Parameters used in the simulations.

Figure 2

Fig. 2. Schematic diagram illustrating the simulation model set-up. An initial structure with certain starting conditions builds the base of the system with a background porosity, temperature, fluid pressure and pore fluid with concentration of species and density. The model then applies boundary conditions fluid pressure (Pf), temperature (T) and concentration of different species of the incoming fluid (c). The next step in the model is the evaluation of the fluid pressure evolution per time step followed by transport of temperature and concentration by advection and diffusion. The local concentration of species and temperature are then used in iphreeqc to calculate the local fluid properties including density and saturation indices of minerals.

Figure 3

Fig. 3. Set-up of the simulations. (a) Mineralogy of the layers with two sandstone layers and one limestone layer in the centre. Initial boundary conditions are defined by equilibrated pore fluid in the different layers. Side walls are confined so that no fluid escapes sideways. The upper boundary is open and contains a low fluid pressure. The metal fluid enters through the lower boundary with an initially increasing and then constant fluid pressure, temperature and incoming fluid composition at the lower boundary. (b) Porosity/permeability scenario I with two permeable faults in the otherwise low-permeability lower layer. (c) Porosity/permeability scenario II with an additional low-permeability horizontally aligned area in the centre of the model above the two faults. (d) Flow lines for scenario II showing how flow concentrates in the two faults and then spreads out in the high-permeability layers and goes around the low-permeability zone in the centre.

Figure 4

Table 2. Fluid compositions used during simulations. Seawater composition is adapted from Ball & Nordstrom (1991) and brine composition from Moldovanyi & Walter (1992) (Smackover Formation brine 55). Element concentrations are given as molarities (per kilogram seawater or kgw) as this is the standard concentration unit used by phreeqc. Alkalinity is defined as HCO3.

Figure 5

Fig. 4. Ideal mixing of pore fluid (seawater) and infiltrating fluid (oil-field brine) calculated using MIX in the iphreeqc module. pH, barite saturation state and CBE are obtained for three databases: phreeqc (solid lines), wateqf4 (dotted lines) and llnl (dashed lines). The blue lines represent the mixing with unaltered seawater and the orange lines show the resulting fluid if mixed with seawater equilibrated with calcite. Equilibration with quartz did not alter the results compared with the non-equilibrated case; only the experiments with non-equilibrated and calcite-equilibrated pore fluid are shown.

Figure 6

Fig. 5. Evolution of fluid in two vertical faults with blue low and red high numbers; width of box is 5 m. (a) Evolution of the fluid pressure into the two faults. (b) Influx of metal-rich fluid into the faults with a Péclet number (see Section 4b for definition) of around 8000 for concentration in the fault. (c) Pore-fluid is displaced in the faults by the incoming fluid. (d) Evolution of barite saturation index in the faults. (e) Barite saturation index summed up over time to show how long the saturation index is high. (f) Temperature after 2 min of model run. Temperature is equilibrating as a result of fast diffusion (T-Péclet number of 100 at the start, and 8.3 after 11 days in the fault).

Figure 7

Fig. 6. Evolution of pore fluid, metal fluid and barite saturation index in one fault. (a, b) Cross-section through a fault where the metal fluid has infiltrated: (a) fluids and barite solution index in three plots and (b) all confined to one plot; mixing best in the fault walls. (c) Vertical section through the incoming metal fluid showing how the pore fluid is pushed away and the barite saturation index is highest within the area where 80% of pore fluid is still present. This wave of high barite saturation index is moving with the influxing fluid through the fault. (d) Two evolutions for the incoming fluid where the concentration curve flattens out into the wall rock as a result of diffusion. (e) Saturation index of barite for (d), showing how the maxima move with diffusion into the wall rock of the fault. (f) Time evolution of a point in the model with a very fast increase in barite saturation index when the fluids mix, followed by a relaxation when the metal fluid dominates (full model run is 231 days).

Figure 8

Fig. 7. Evolution of summed saturation index of barite in the faults: (a–l) progressive time steps. The saturation index shows a relatively stationary and diffusion-controlled wave of high accumulated SI of barite in the wall rock, but a travelling wave of high SI of barite in the faults that grow in width.

Figure 9

Fig. 8. Effects of horizontal permeable layer and a horizontal low-permeability layer on flux and saturation index of barite. (a) Accumulated barite saturation index in a simulation with a permeable layer on top of the two faults but no seal. High saturation index above the fault is present on the fault sides and where fluids between the two faults mix as a result of diffusion. (b) Concentration of incoming metal fluid for scenario I. (c) Accumulated barite saturation index in a simulation with a horizontal permeable layer covered by a central low-permeability layer (scenario II). (d) Concentration of metal fluid in scenario II. Fluid is pushed against the barrier and flows around it at the left- and right-hand side. Below the barrier the accumulated saturation index of barite is high, indicating that barite would precipitate at this location. All pictures shown are after 115 days.

Figure 10

Fig. 9. Time evolution of the accumulated saturation index of barite in simulations with a permeable layer and a barrier (scenario II). The highest indices are in the fault walls, between and at the exit of faults, in the permeable layer and below the barrier.

Figure 11

Fig. 10. (a) Initial evolution of the Péclet number (Pe) with time in the simulations. Initially, the Péclet number for concentration of species is decaying fast. Inset shows a log-log plot of Péclet numbers for concentration versus temperature, illustrating that temperature is a lot more diffusion controlled than concentration. (b) Schematic illustration of a high-permeability fluid channel or fault with a travelling wave of high saturation index of barite. Stable waves develop in the wall rock of the fault where mineralization can occur. A plot of SI barite along the fault on the left-hand side of the figure shows a snapshot of the travelling wave of high-SI barite, whereas the plot at the bottom of the figure shows the stationary waves of high-SI barite in the wall rock of the faults. The index is only high where the mixing between the fluids is optimal, with c. 80–90% of pore fluid and only a small amount of incoming fluid of 10–20%. The index therefore decays once the incoming fluid dominates and displaces the pore fluid. The geometries are similar to the experiment shown in Figure 1c, even though the experiment did not directly involve fluid mixing.

Supplementary material: Image

Koehn et al. supplementary material

Koehn et al. supplementary material 1

Download Koehn et al. supplementary material(Image)
Image 2.1 MB
Supplementary material: Image

Koehn et al. supplementary material

Koehn et al. supplementary material 2

Download Koehn et al. supplementary material(Image)
Image 5 MB