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

Daniel Koehn, Ulrich Kelka, Renaud Toussaint, Coralie Siegel, Gary Mullen, Adrian Boyce, Sandra Piazolo

<figure><img src="image_1.png" /></figure>

## To cite this version:

Daniel Koehn, Ulrich Kelka, Renaud Toussaint, Coralie Siegel, Gary Mullen, et al.. Outcrop scale mixing enhanced by permeability variations: the role of stationary and travelling waves of high saturation indices. Geological Magazine, 2022, Faults and fractures in rocks: mechanics, occurrence, dating, stress history and fluid flow, 159 (11-12), pp.2279-2292. 10.1017/S001675682200070X. hal-04311857

HAL Id: hal-04311857

https://hal.science/hal-04311857v1

Submitted on 28 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Outcrop scale mixing enhanced by permeability variations: the role of stationary and travelling waves of high saturation indices

Daniel Koehn¹ ID, Ulrich Kelka², Renaud Toussaint³,⁴ ID, Coralie Siegel², Gary Mullen⁵,⁶, Adrian Boyce⁶ and Sandra Piazolo⁷

## Original Article

**Cite this article:** Koehn D, Kelka U, Toussaint R, Siegel C, Mullen G, Boyce A, and Piazolo S (2023) Outcrop scale mixing enhanced by permeability variations: the role of stationary and travelling waves of high saturation indices. *Geological Magazine* **159**: 2279–2292. https://doi.org/10.1017/S001675682200070X

Received: 30 January 2022
Revised: 30 May 2022
Accepted: 20 June 2022
First published online: 13 October 2022

**Keywords:**
diffusion; advection; saturation index; phreeqc; modelling

**Author for correspondence:**
Daniel Koehn,
Email: daniel.koehn@fau.de

¹GeoZentrum Nordbayern, Friedrich Alexander University (FAU) Erlangen-Nuremberg, Schlossgarten 5, Erlangen 91054, Germany; ²CSIRO – Mineral Resources, 26 Dick Perry Ave, Kensington, WA 6151, Australia; ³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; ⁵School of Geographical and Earth Sciences, Glasgow University, G12 8QQ, UK; ⁶Scottish Universities Environmental Research Centre, East Kilbride, Glasgow, G75 0QF, UK and ⁷School of Earth and Environment, University of Leeds, LS2 9JT, UK

## 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^3$ 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_2$ 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

© The Author(s), 2022. Published by Cambridge University Press. 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.

<figure><img src="image_3.png" /></figure>

CAMBRIDGE
UNIVERSITY PRESS
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).

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 $\phi_{x,y}$ to the fluid grid, which is used to calculate the permeability $\kappa(\phi_{x,y})$ (m$^2$) assuming the Carman-Kozeny relation (Carman, 1937; Phillips, 1991; Ghani et al. 2013) as

$$ \kappa(\phi_{x,y}) = \frac{r^2(\phi_{x,y})^3}{45(1 - \phi_{x,y})^2} \quad (1) $$

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

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

where $\phi$ is porosity, $\beta$ is fluid compressibility (m$^2$/N), $P$ is fluid pressure (Pa), $\kappa$ is permeability (m$^2$) and $\mu$ 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

$$ \vec{v} = - \frac{\kappa \nabla P}{\mu \phi} \quad (3) $$

The fluid velocity is then used in the general transport equation to derive the advective part for temperature and the different chemical species
<figure><img src="image_4.png" /><figcaption>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).</figcaption></figure>

$$ \frac{\partial(C, T)}{\partial t} + \vec{v} \nabla(C, T) - D\Delta(C, T) = 0 \quad (4) $$

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

$$ (C, T)^t = (C, T)^{t-1} + \delta(C, T)_{adv}^t + \delta(C, T)_{diff}^t \quad (5) $$

where T is temperature (°C), C is concentration of species (mol kg⁻¹) and D is the diffusion coefficient for temperature and matter diffusion respectively (m² s⁻¹). 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.

## 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 SIp in iphreeqc is defined as (Parkhurst & Appelo, 1999):

$$ SI_p = \log \prod_m^{M_{aq}} a_m^{c_{m,p}} \quad (6) $$

where Maq is the total number of aqueous master species, am is the activity of the master species and cm,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 fp used by iphreeqc for phase equilibrium to reach the target SIp target is:

$$ f_p = (\ln K_p + [\ln(10)]SI_{p, target}) - \sum_m^{M_{aq}} c_{m,p} \ln(a_m), \quad (7) $$

where Kp 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
<table><caption>**Table 1.** Parameters used in the simulations.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Model size</td><td>S</td><td>5</td><td>m</td></tr><tr><td>Time step</td><td>t</td><td>200</td><td>s</td></tr><tr><td>Grid size solid</td><td>R<sub>s</sub></td><td>200 × 230</td><td>Model units</td></tr><tr><td>Grid size fluid</td><td>R<sub>f</sub></td><td>100 × 100</td><td>Model units</td></tr><tr><td>Temperature lower boundary</td><td>T</td><td>90</td><td>°C</td></tr><tr><td>Temperature background model and upper boundary</td><td>T</td><td>25</td><td>°C</td></tr><tr><td>Fluid pressure boundary</td><td>P<sub>f</sub></td><td>(2.001 to 3.3) × 10<sup>6</sup></td><td>Pa</td></tr><tr><td>Diffusion constant T fault/wall rock</td><td>D<sub>T</sub></td><td>10<sup>-6</sup> to 10<sup>-7</sup></td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Diffusion constant fault/wall-rock Ba, Ca, Cl, Fe, K, Mg, Mn, Na, Pb, S(6), Si, Zn</td><td>D<sub>c</sub></td><td>10<sup>-10</sup> to 10<sup>-11</sup></td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Porosity</td><td>φ</td><td>0.01–0.10</td><td></td></tr><tr><td>Fluid viscosity</td><td>μ</td><td>1.0 × 10<sup>-3</sup></td><td>Pa s</td></tr><tr><td>Fluid compressibility</td><td>β</td><td>4.5 × 10<sup>-10</sup></td><td>m<sup>2</sup> N<sup>-1</sup></td></tr><tr><td>Fluid density (pore fluids)</td><td>ρ</td><td>1.02323</td><td>kg m<sup>-3</sup></td></tr><tr><td>Fluid density (infiltrating fluid)</td><td>ρ</td><td>1.15907</td><td>kg m<sup>-3</sup></td></tr><tr><td>Carman–Kozeny grain size</td><td>r</td><td>0.00001</td><td>M</td></tr><tr><td>Permeability</td><td>κ</td><td>(2.37 × 10<sup>-18</sup>) to (2.74 × 10<sup>-15</sup>)</td><td>m<sup>2</sup></td></tr><tr><td>Density of rocks</td><td>ρ<sub>r</sub></td><td>2.700</td><td>kg m<sup>-3</sup></td></tr><tr><td>Gravitational acceleration</td><td>g</td><td>9.81</td><td>m s<sup>-2</sup></td></tr></tbody></table>

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⁻³ for the overlying rock and acceleration due to gravity of 9.81 m s⁻². 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.

The model geometry has three horizontal layers: a lower sandstone layer overlain by a limestone and another sandstone layer

<figure><img src="image_5.png" /><figcaption>**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<sub>f</sub>), 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.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

(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 \times 10^{-18} \text{ m}^2$, whereas the other two layers have a high porosity of 0.10 and a higher permeability of $2.75 \times 10^{-15} \text{ m}^2$ (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 \times 10^{-15} \text{ m}^2$. 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 \times 10^{-18} \text{ m}^2$ (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⁻¹). 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
<table><caption>**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₃.</caption><thead><tr><th></th><th>Seawater</th><th>Brine</th></tr></thead><tbody><tr><td>pH</td><td>8.22</td><td>5.370</td></tr><tr><td>Temp (°C)</td><td>25.0</td><td>90</td></tr><tr><td>Error</td><td>0.07</td><td>-4.05</td></tr><tr><td>Alkalinity</td><td>2.406 × 10<sup>-3</sup></td><td>3.185 × 10<sup>-3</sup></td></tr><tr><td>Ba (mol/kgw)</td><td>-</td><td>3.780 × 10<sup>-4</sup></td></tr><tr><td>Ca (mol/kgw)</td><td>1.066 × 10<sup>-2</sup></td><td>1.053</td></tr><tr><td>Cl (mol/kgw)</td><td>5.657 × 10<sup>-1</sup></td><td>5.932</td></tr><tr><td>Fe (mol/kgw)</td><td>-</td><td>6.221 × 10<sup>-3</sup></td></tr><tr><td>K (mol/kgw)</td><td>1.058 × 10<sup>-2</sup></td><td>3.404 × 10<sup>-2</sup></td></tr><tr><td>Mg (mol/kgw)</td><td>5.507 × 10<sup>-2</sup></td><td>-</td></tr><tr><td>Mn (mol/kgw)</td><td>-</td><td>2.023 × 10<sup>-3</sup></td></tr><tr><td>Na (mol/kgw)</td><td>4.854 × 10<sup>-1</sup></td><td>3.312</td></tr><tr><td>Pb (mol/kgw)</td><td>-</td><td>2.223 × 10<sup>-4</sup></td></tr><tr><td>S(6) (mol/kgw)</td><td>2.926 × 10<sup>-2</sup></td><td>9.422 × 10<sup>-4</sup></td></tr><tr><td>Si (mol/kgw)</td><td>7.382 × 10<sup>-5</sup></td><td>6.646 × 10<sup>-4</sup></td></tr><tr><td>Zn (mol/kgw)</td><td>-</td><td>3.502 × 10<sup>-3</sup></td></tr></tbody></table>

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

<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

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.

## 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><img src="image_8.png" /><figcaption>**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).</figcaption></figure>

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 1: 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).

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
**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><img src="image_9.png" /></figure>

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.

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.

### 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.

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.

## 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
<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

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:

$$Pe = \frac{vL}{D} \qquad (8)$$

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
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

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

$$Da_I = \frac{R}{v}, \qquad (9)$$

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

$$Da_{II} = \frac{RL}{D} \qquad (10)$$

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^{-10}$ m s$^{-1}$ for barite growth (Kuwahara et al. 2016) the Damköhler I number is $c \cdot 10^{-4}$ with a fluid velocity of $c \cdot 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 \cdot 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
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

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 geo-physical-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/R00772/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.

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.

Belousov BP (1959) Periodicheski deistvuyushchaya reaktsia i ee mekhanism [Periodically acting reaction and its mechanism]. Sbornik referatov po radio-tsionnoi meditsine, 1958 [Collection of abstracts on radiation medicine, 1958]. Moscow: Medgiz, 145–47.

Bjørlykke K (1993) Fluid flow in sedimentary basins. *Sedimentary Geology* 86, 137–58.

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.

Bons PD, Koehn D and Jessell MW (2008) *Microdynamics Simulation*. New York: Springer. Lecture Series in Earth Sciences, p. 406.

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.

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.

Carman PC (1937) Fluid flow through granular beds. *Transactions of the Institution of Chemical Engineers* 15, 150.

Chapman B, Jost G and Van der Pas R (2007) *Using OpenMP: Portable Shared Memory Parallel Programming*. Cambridge, MA: MIT Press.

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.

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.
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.

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.

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.

Faraday M (1827) *Chemical Manipulation, Being Instructions to Students in Chemistry*. London: W. Phillips, 656 p.

Fredrich JT, Digiovanni AA and Noble DR (2006) Predicting macroscopic transport properties using microscopic image data. *Journal of Geochemical Research-Solid Earth* **111**. doi:10.1029/2005JB003774.

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.

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.

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.

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.

Ghani I, Koehn D, Toussaint R and Passchier CW (2013) Dynamic development of hydrofracture. *Pure and Applied Geophysics* **170**, 1685–703.

Ghani I, Koehn D, Toussaint R and Passchier C (2015) Dynamics of hydrofracturing and permeability evolution in layered reservoirs. *Frontiers in Physics* **3**, 67.

Glansdorff & Prigogine (1971) *Thermodynamic Theory of Structure, Stability and Fluctuations*. Hoboken: Wiley-Interscience, 306p.

Haynes DW, Cross KC, Bills RT and Reed MH (1995) Olympic Dam ore genesis: a fluid-mixing model. *Economic Geology*, **90**, 281–307.

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.

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 K Gessner, TG Blenkinsop and P Sorjonen-Ward), pp. 147–64. Geological Society of London, Special Publication no. 453.

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.

Hrennikoff A (1941) Solution of problems of elasticity by the framework method. *Journal of Applied Mechanics* **8**, A619–715.

Jamtveit B, Putnis CV and Malthe-Sørenssen A (2009) Reaction induced fracturing during replacement processes. *Contributions to Mineralogy and Petrology* **157**, 127–33.

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.

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.

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.

Korrani AKN, Sepehrnoori K and Delshad M (2015) Coupling IPhreeqc with UTCHEDM to model reactive flow and transport. *Computers & Geosciences* **82**, 152–69.

Kuwahara Y, Liu W, Makio M and Otsuka K (2016) In situ AFM study of crystal growth on a barite (001) surface in BaSO₄ solutions at 30°C. *Minerals* **6**, 117.

Lax PD and Wendroff B (1960) Systems of conservation laws. *Communications on Pure and Applied Mathematics* **13**, 217–37.

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.

Lester DR, Ord A and Hobbs BE (2012) The mechanics of hydrothermal systems: II. Fluid mixing and chemical reactions. *Ore Geology Reviews* **49**, 45–71.

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.

Loukola-Ruskeenemiemi K (1999) Origin of black shales and the serpentinite-associated Cu-Zn-Co ores at Outokumpu, Finland. *Economic Geology* **94**, 1007–28.

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.

Moldovanyi EP and Walter LM (1992) Regional trends in water chemistry, Smackover Formation, southwest Arkansas: geochemical and physical controls. *AAPG Bulletin* **76**, 864–94.

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.

Mostaghimi P, Liu M and Arns CH (2016) Numerical simulation of reactive transport on micro-CT images. *Mathematical Geosciences* **48**, 963–83.

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**, 10–21.

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.

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, A Schmidt, J Brugger), pp. 833–35. Geneva: Society for Geology Applied to Mineral Deposits.

Ortoleva P, Chadam J, Merino E and Sen Z (1987) Geochemical self-organization ii: the reactive-infiltration instability. *American Journal of Science* **287**, 1008.

Parkhurst DL and Appelo C (1999) User's guide to PHREEQC (Version 2). *Water-Resources Investigations Report* **99**, 4259.

Phillips OM (1991) *Flow and Reactions in Permeable Rocks*. Cambridge: Cambridge University Press.

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.

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.

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.

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.

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.

**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.

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.

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.

Szymczak P and Ladd AJC (2009) Wormhole formation in dissolving fractures. *Journal of Geophysical Research* **114**. doi: 10.1029/2008JB006122.

Tartakovsky AM and Meakin P (2006) Pore scale modeling of immiscible and miscible smoothed particle hydrodynamics. *Advances in Water Resources* **29**, 1464–78.

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.

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.

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.

Wilkinson J (2001) Fluid inclusions in hydrothermal ore deposits, Lithos 55, 229–72.

Zaikin AN and Zhabotinsky AM (1970) Concentration wave propagation in two-dimensional liquid-phase self-oscillating system. *Nature* **225**, 535–7.

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.