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

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 (10³ 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.

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. 1995; Fu et al. 2003; Fan et al. 2006; Neumayr et al. 2008; Gessner et al. 2009; Leach et al. 2010; Lester et al. 2012; Hobbs & Ord, 2018), but is still not well understood (Ord et al. 2010; Lester et al. 2012). 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. 2012), where changes in redox potential, O₂ fugacity, pH and temperature can trigger mineral precipitation (Liu et al. 2021; Weihua et al. 2021). 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. 2006). Following Lester et al. (2012), 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. 2012). 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. 2013; Bons et al. 2014), where fluid is thought to come up in pulses and ores will precipitate close to an unconformity (Wagner et al. 2010). 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, 1999). 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. 2013; Bons et al. 2014). 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. 2019), 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, 2001; Fusswinkel et al. 2013). 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. 2010; Lester et al. 2012).

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. 2021) with the hydrogeochemical simulation module iphreeqc (Charlton & Parkhurst, 2011). Similar approaches in creating reactive transport modelling environments were reported decades ago (Steefel & Lasaga, 1994) and were performed by utilizing phreeqc as a reaction rate engine that is coupled with UTCHEM (Korrani et al. 2015), MATLAB (Shabani et al. 2019), COMSOL (Nardi et al. 2014; Wei & Cao, 2021) and POET (De Lucia et al. 2021), as well as other reactive transport codes (see Steefel et al. 2015 for a review). The evolution of permeability in porous media and rocks has been modelled extensively (Saripalli et al. 2001; Zhao et al. 2007; Jamtveit et al. 2009; Chen et al. 2014; Kang et al. 2014; Mostaghimi et al. 2016) with smooth particle hydrodynamics, lattice Boltzmann methods and computational fluid dynamic techniques (Manwart et al. 2002; Fredrich et al. 2006; Tartakovsky & Meakin, 2006; Shabro et al. 2012; Chen et al. 2013). 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. 2009, 2018; Hobbs & Ord, 2018; Moosavi et al. 2019). 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. 2008; Ghani et al. 2013, 2015; Piazolo et al. 2019; Koehn et al. 2020, 2021) with iphreeqc (Charlton & Parkhurst, 2011). 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, 1941). 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. 1992) 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 ) (m²) assuming the Carman–Kozeny relation (Carman, 1937; Phillips, 1991; Ghani et al. 2013) 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 (m²/N), P is fluid pressure (Pa), κ is permeability ( ${m^2})$ and μ is fluid viscosity (Pa s; Ghani et al. 2013). 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 (m² s^–1). The IMEX (IMplicit+Explicit; Ascher et al. 1997) 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, 2011) 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, 1999):

(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, 1992). 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, 1999). For simplicity we define alkalinity as HCO₃^– 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. 2007).

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 (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 HCO₃.

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⁻¹⁸ m², whereas the other two layers have a high porosity of 0.10 and a higher permeability of 2.75 × 10⁻¹⁵ m² (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⁻¹⁵ m². 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⁻¹⁸ m² (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 (1991) and a representative example of a brine (II) from the Smackover formation in SW Arkansas (Moldovanyi & Walter, 1992). 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, 1969; Bjørlykke, 1993). 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. 2013) 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. 2012).

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. 2013). 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, 1960) 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, 2017), and are important during mixing processes in ore deposits (Lester et al. 2012). One form of these structures are waves and spots of reactants that are stationary (e.g. Turing patterns) or travelling (Belousov, 1959; Zaikin & Zhabotinsky, 1970). As discussed by Lester et al. (2012), 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, 1959; Zaikin & Zhabotinsky, 1970) (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 5–7) 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. (2012) 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. 1987; Bons et al. 2013; Koehn et al. 2021; 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, 2009); the Damköhler numbers also influence reactive transport for the alteration of pore space (Mostaghimi et al. 2016), as well as the roughness of reaction fronts (Koehn et al. 2021).

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. 2013). 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⁻¹⁰ m s^–1 for barite growth (Kuwahara et al. 2016) the Damköhler I number is c. 10⁻⁴ with a fluid velocity of c. 10⁻⁶ 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.