# Hydrological and Geological Controls for the Depth Distribution of Dissolved Oxygen and Iron in Silicate Catchments

Ivan Osorio-leon, Camille Bouchez, Eliot Chatton, Nicolas Lavenant, Laurent

Longuevergne, Tanguy Le Borgne

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

## To cite this version:

Ivan Osorio-leon, Camille Bouchez, Eliot Chatton, Nicolas Lavenant, Laurent Longuevergne, et al. Hydrological and Geological Controls for the Depth Distribution of Dissolved Oxygen and Iron in Silicate Catchments. Water Resources Research, 2023, 59 (8), pp.e2023WR034986. 10.1029/2023wr034986. insu-04182271

HAL Id: insu-04182271

https://insu.hal.science/insu-04182271

Submitted on 17 Aug 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>
# Water Resources Research®

## RESEARCH ARTICLE

10.1029/2023WR034986

### Special Section:

The Future of Critical Zone
Science: Towards Shared
Goals, Tools, Approaches and
Philosophy

### Key Points:

* We derive a reduced analytical model to predict the depth distribution of dissolved O₂ and Fe²⁺ in silicate catchments

* The theory, validated from a reactive transport model, quantifies the hydrological and geological controls through two dimensionless numbers

* We use this framework to analyze and explain the O₂/Fe²⁺ distributions in two sites with contrasted chemical properties

### Supporting Information:

Supporting Information may be found in the online version of this article.

### Correspondence to:

I. Osorio-Leon and C. Bouchez,
idosoriole@unal.edu.co;
camille.bouchez@univ-rennes.fr

### Citation:

Osorio-Leon, I., Bouchez, C., Chatton, E., Lavenant, N., Longuevergne, L., & Le Borgne, T. (2023). Hydrological and geological controls for the depth distribution of dissolved oxygen and iron in silicate catchments. *Water Resources Research*, 59, e2023WR034986. https://doi.org/10.1029/2023WR034986

Received 28 MAR 2023

Accepted 29 JUL 2023

### Author Contributions:

**Conceptualization:** Ivan Osorio-Leon, Camille Bouchez, Tanguy Le Borgne

**Data curation:** Ivan Osorio-Leon, Camille Bouchez, Eliot Chatton, Nicolas Lavenant, Tanguy Le Borgne

© 2023. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

## Hydrological and Geological Controls for the Depth Distribution of Dissolved Oxygen and Iron in Silicate Catchments

Ivan Osorio-Leon¹ ID, Camille Bouchez¹ ID, Eliot Chatton¹ ID, Nicolas Lavenant¹, Laurent Longuevergne¹ ID, and Tanguy Le Borgne¹ ID

¹University Rennes—CNRS, Géosciences Rennes—UMR 6118, Rennes, France

**Abstract** Dissolved Oxygen (DO) plays a key role in reactive processes and microbial dynamics in the critical zone. Recent observations showed that fractures can provide rapid pathways for oxygen penetration in aquifers, triggering unexpected biogeochemical processes. In the shallow subsurface, DO reacts with electron donors, such as Fe²⁺ coming from mineral dissolution. Yet, little is known about the factors controlling the spatial heterogeneity and distribution of oxygen with depth. Here we present a reduced analytical model describing the coupled evolution of DO and Fe²⁺ as a function of fluid travel time in silicate catchments. Our model, validated from fully resolved reactive transport simulations, predicts a linear decay of DO with time, followed by a rapid non-linear increase of Fe²⁺ concentrations up to a far-from-equilibrium steady-state. The relative effects of geological and hydrological forcings are quantified through a Damköhler number (Da) and a lithological number (Λ). We use this framework to investigate the depth distribution of DO and Fe²⁺ in two catchments with similar environmental contexts but contrasted hydrochemical properties. We show that hydrochemical differences are explained by small variations in Da but orders of magnitude variations in Λ. Therefore, we demonstrate that the hydrological and geological drivers controlling hydrochemistry in silicate catchments can be discriminated by analyzing jointly the O₂ and Fe²⁺ evolution with depth. These findings provide a new conceptual framework to understand and predict the evolution of DO in modern groundwater, which plays an important role in critical zone processes.

**Plain Language Summary** Dissolved Oxygen (DO) is involved in important biogeochemical reactions, such as the aerobic respiration of microbes, rock-weathering or contaminant degradation. The general view is that the deeper subsurface of continents is mostly anoxic. However, recent observations have shown that cracks in rocks can allow oxygen to reach deeper down, thus extending deeper the influence of oxygen in reactions. Here, we present a modeling framework validated with field data to understand and predict the hydrological and geological controls on DO evolution in silicate catchments, shedding new light on its influence on rock-weathering and microbial life in the continents.

## 1. Introduction

Oxygen is central in redox reactions because it is the most abundant and readily available electron acceptor in the environment (Korom, 1992; Stumm & Morgan, 1996) and offers a strong redox potential. In aquatic environments, dissolved oxygen (DO) is mostly produced by photosynthesis and consumed by aerobic respiration of organic matter (Mader et al., 2017). In the subsurface, the transport of DO and CO₂ by fluid flow triggers the weathering of silicate rocks (Fletcher et al., 2006; Kim et al., 2017; Li et al., 2017; Singha & Navarre-Sitchler, 2022), which represent a quarter of the outcropping Earth rocks (Hartmann & Moosdorf, 2012). The cycle of DO is thus closely related to the geochemical cycles of carbon (Bar-on et al., 2018; Petsch et al., 2004), iron (Kappler et al., 2021; Melton et al., 2014; Napieralski et al., 2019) and sulfur (Canfield, 2004; Gu et al., 2020). Besides abiotic processes, recent studies have pointed out that weathering reactions, particularly involving DO, are often mediated by microbes (Erable et al., 2012; Kappler et al., 2021; Napieralski et al., 2019) because the redox potential of the reaction offers an energy source for microbes to thrive (Emerson et al., 2010). DO thus also exerts a key ecological role by impacting the biodiversity (Hancock et al., 2005; Humphreys, 2009; Malard & Hervant, 1999) and activity of aerobic microbial metabolisms (Druschel et al., 2008; Mader et al., 2017; Maisch et al., 2019).

When DO is not depleted in soils, either because of a limited soil-thickness or a low organic matter availability, oxic water can enter the saturated zone in the bedrock. It then reacts with the electron donors available
**Formal analysis:** Ivan Osorio-Leon, Camille Bouchez, Eliot Chatton, Nicolas Lavenant

**Funding acquisition:** Camille Bouchez, Tanguy Le Borgne

**Investigation:** Ivan Osorio-Leon, Camille Bouchez, Tanguy Le Borgne

**Methodology:** Ivan Osorio-Leon, Camille Bouchez, Eliot Chatton, Tanguy Le Borgne

**Project Administration:** Camille Bouchez, Laurent Longuevergne, Tanguy Le Borgne

**Resources:** Camille Bouchez, Nicolas Lavenant

**Software:** Ivan Osorio-Leon

**Supervision:** Camille Bouchez, Tanguy Le Borgne

**Validation:** Ivan Osorio-Leon, Camille Bouchez, Tanguy Le Borgne

**Visualization:** Ivan Osorio-Leon

**Writing – original draft:** Ivan Osorio-Leon, Camille Bouchez, Tanguy Le Borgne

**Writing – review & editing:** Ivan Osorio-Leon, Camille Bouchez, Laurent Longuevergne, Tanguy Le Borgne

in the subsurface, such as Fe(II), Mn(II), CH₄, H₂, HS⁻ or petrogenic organic carbon (Canfield, 2004; Kartsen Pedersen, 1997; Tebo et al., 2005). Iron is the most abundant redox-sensitive element in the Earth crust (Frey & Reed, 2012), it is linked to biogeochemical cycles of carbon, sulfur and nitrogen (Casar et al., 2021; Kendall et al., 2012). The reduction of DO by iron can occur either by Fe(II)-sites on mineral surfaces (White et al., 1985) or by Fe²⁺ dissolved in water subsequently to the release of structural Fe(II) by mineral dissolution (White & Yee, 1985). The most common Fe(II)-bearing primary minerals in hard-rocks are silicates such as biotite (Aquilina et al., 2018; Hampl et al., 2021; Holbrook et al., 2019; Malmström et al., 1996), pyroxene (Behrens et al., 2015) and hornblende (Fletcher et al., 2006) and in less proportion, sulfides such as pyrite (Gu et al., 2020).

Several field works documenting weathering profiles in hard-rocks, either from outcrops (Antoniellini et al., 2017) or borehole cores (Bazilevskaya et al., 2013; Dideriksen et al., 2010; Hampl et al., 2021; Holbrook et al., 2019), have suggested that DO transport by subsurface flow could explain the presence of secondary minerals and weathering induced fracturing (WIF) in ferrous silicates (Bazilevskaya et al., 2013; Kim et al., 2017) from deep regolith. Likewise, evidence from field measurements in fractured-rock aquifers has shown that DO can effectively persist in deep aquifers (Bucher et al., 2009; DeSimone et al., 2014; Edmunds et al., 1984; Ruff et al., 2023; Sullivan et al., 2016; Winograd & Robertson, 1982). Based on field observations, the presence or absence of DO in the subsurface has been attributed to lithological differences of the bedrock (Malard & Hervant, 1999; Winograd & Robertson, 1982) or to the spatial distribution of oxidation fronts (Liao et al., 2022). However, it is still challenging to predict the expected depth of DO in silicate catchments.

Fractures can provide fast transport pathways in the subsurface and therefore influence reactive transport processes (Deng & Spycher, 2019). The transport and fate of DO in fractured rocks have been studied in the framework of risk assessments for DO penetration to nuclear waste repository sites. Numerical studies have simulated the advance of the redox front in the matrix of granitic rocks (Macquarrie et al., 2010; Spiessl et al., 2008; Trinchero et al., 2019) and in fracture networks (Trinchero et al., 2017) providing insights into the mechanisms driving oxygen transport in fractured rocks. Three main approaches have been used to model the reactive transport of oxygen in fractured media (Macquarrie et al., 2010): (a) Equivalent porous medium (EPM) models for fracture zones, (b) discrete fracture models with reactive in fill and no interaction with the adjacent rock matrix, and (c) fracture–matrix models in which the fracture and/or matrix may contain reactive minerals. The effect of matrix diffusion is particularly important in driving the long term (thousands of years) transport of oxygen, relevant for applications such as nuclear waste storage. In this context, approximate analytical solutions have been obtained (Sidborn & Neretnieks, 2007, 2008) by assuming that the dissolution of Fe(II) bearing minerals is rapid compared to transport, which is relevant for large time scales (thousands of years).

In the context of catchment dynamics, modern groundwater residence times are on the order of tens to hundreds of years (Sprenger et al., 2019). At these time scales, DO diffusion in the matrix may be assumed to be relatively small compared to advective fluxes while dissolution rates cannot assumed to be fast (Sidborn & Neretnieks, 2008; Trinchero et al., 2018). In this regime, how the interplay between dissolution rates and advective transport dynamics drives DO transport has not been quantified analytically. The coupling between transport and reaction rates at catchment scale has been conceptualized using the Damköhler number, defined as the ratio of the characteristic transport and reaction time scales (Maher, 2010; Maher & Chamberlain, 2014). In this approach, structures such as fractures and rock matrix are not represented explicitly but instead integrated into an effective fluid travel time, representing the time during which fluid has been exposed to reactive minerals (Seeboonruang & Ginn, 2006). In this study, we use a travel time formulation to develop reduced analytical model for coupled transport and reaction of DO and Fe²⁺ in silicate catchments.

We thus develop approximate analytical solutions for coupled oxygen and iron transport and reaction in a simple conceptual water-rock interaction model. We use a travel time formulation and a simple setting of two minerals accounting for dissolution and precipitation processes. We then use a full numerical model with a more realistic set of minerals to validate the approximations of the analytical approach. The reduced analytical model provides a mechanistic understanding of geological and hydrological controls of oxygen transport and reaction through two dimensionless numbers: Damköhler (Da) and a Lithological number (Λ). The former quantifies the relative effect of Fe(II)-bearing minerals abundance and transport velocity, while the latter quantifies the balance between dissolution and precipitation of primary and secondary minerals, respectively. We use this framework to interpret field data in two critical zone observatories with contrasted chemical properties.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Conceptual model for the transport and reactivity of Dissolved Oxygen (DO) and Fe²⁺ along a flowpath (a). (b) In the unweathered rock, DO is consumed owing to the following reaction network (i) Fe(II)-bearing minerals dissolution, promoted by groundwater acidity, releases Fe²⁺ (ii) aqueous oxidation of dissolved iron by DO and precipitation of Fe(III) oxyhydroxydes (iii) incongruent mineral dissolution releases Fe²⁺ and forms Fe(II)-bearing clay. OAT indicates the Oxic-Anoxic Transition along the pathway.</figcaption></figure>

## 2. Reactive Transport of Dissolved Oxygen and Iron

### 2.1. Model Conceptualization

The different conceptual models used for the reactive transport of oxygen in fractured media may consider advection in the fracture, reaction in the fracture, diffusion and reaction in the matrix (Macquarrie et al., 2010). It is not known how to solve the problem analytically with all these processes. However, analytical solutions can be derived by considering the dominant processes at different time scales. At long time scales (thousands of years), matrix diffusion is a key process, while reaction can be considered to be fast, which facilitates the derivation of analytical expressions for the oxic front propagation (Sidborn & Neretnieks, 2007, 2008). Here we focus on the short term regime (10¹⁻10² years), consistent with catchment residence times. We thus neglect the effect of matrix diffusion but include the impact of reaction kinetics, leading to both transport-limited and reaction-limited dynamics. Therefore, our model considers a fracture zone as an equivalent porous media with reactive infill and neglects matrix diffusion, following a common assumption (Trinchero et al., 2018).

To unravel the respective roles of hydrological and geological controls on DO and Fe²⁺, we consider a travel time formulation (Maher, 2010). For a given flowpath, we model a reaction network in a parcel of transported fluid (i.e., we use a Lagrangian framework) as a kinetic system controlled by travel time (Figure 1). We thus assume that transport and mineral dissolution reactions are at steady state and identify travel time with reaction time in the model. We quantify interactions between oxic recharge water and Fe(II)-bearing minerals in a fracture zone. We assume advection dominated transport in the fracture zone and thus neglect diffusive DO inputs from the surface. We derive approximate analytical expressions of both DO and Fe²⁺ as a function of fluid travel time. We assume that the flowpath crosses a shallow regolith zone followed by an unweathered zone in which DO reacts with Fe²⁺ produced from minerals dissolution (Figure 1). In a second step, we relate travel time to depth using an average transport velocity and derive equations for the evolution with depth of DO and Fe²⁺.

### 2.2. Geochemical System

We study the kinetic coupling between two geochemical reactions: the dissolution of iron-bearing minerals to produce dissolved iron (Fe²⁺) followed by its oxidation with DO transported from recharge fluids (Figure 1). Even
if the oxidation of structural Fe(II) has been evidenced for Fe(II)-bearing minerals like Biotite (Buss et al., 2008) and Chlorite (Liao et al., 2022), here we neglect this process since its kinetics are still poorly constrained. Therefore, and similarly to previous works (Macquarrie et al., 2010; Sidborn & Neretnieks, 2008), we assume that Fe(II)-bearing minerals dissolve through a non-oxidative mechanism that liberates Fe²⁺ while consuming acidity, and that Fe(II)-bearing clays (Fe-clay) precipitate.

We model the dissolution and precipitation kinetics for mineral reactions using the Transition State Theory, expressing the reaction rate for a mineral j ($R_j$) as (Lasaga, 1984):

$$R_j = S_{V,j} \omega^{-1} (k_{j,H} a_H^{n_H} + k_j + k_{j,OH} a_{OH}^{n_{OH}}) e^{-\frac{E_{a,j}}{RT}} \left(1 - \frac{Q_j}{K_{sp,j}}\right), \\ \text{with } S_{V,j} = \frac{S_{M,j} \Phi_j M M_j}{V_{m,j}} \quad (1)$$

With ($k_{j,H}$, $k_j$, $k_{j,OH}$) [mol m⁻² s⁻¹] the intrinsic reaction constants at 25°C at acid, neutral and basic pH, respectively; a the activity and n the affinity factor of the indicated ion accounting for reaction catalysis by pH, $E_{a,j}$ the activation energy, $Q_j$ the activity product and $K_{sp,j}$ the solubility product for mineral species j, $\omega$ is the fracture porosity, R is the ideal gas constant and T the absolute temperature. $S_{V,j}$ [$m²·m_{pm}⁻³$] is the bulk surface area of the mineral species j, a function of the mineral specific surface area $S_{M,j}$, the mineral volume fraction $\Phi_j$, the mineral molecular mass $MM_j$ and the mineral molar volume $V_{m,j}$. A glossary containing the main variables in this manuscript is available in Table 1.

The aqueous reaction between DO and dissolved iron, and the corresponding kinetic law, are respectively:

$$4Fe^{2+} + O_2 + 4H^+ \rightarrow 4Fe^{3+} + 2H_2O \quad (2)$$

$$R_{ox} = k_{ox}C_{Fe^{2+}}P_{O_2}C_{OH-}^2 \quad \text{for pH > 4.5,} \quad (3)$$

with $R_{ox}$ the rate of Fe²⁺ oxidation, $k_{ox}$ the intrinsic reaction constant of oxidation (1.3 × 10¹²M⁻² atm⁻¹ s⁻¹), $C_i$ the concentration of species i in [mol L⁻¹] and $P_{O_2}$ the partial pressure of oxygen in [atm] (Singer & Stumm, 1970).

Since groundwater in silicate catchments is commonly slightly acid to near-neutral (DeSimone et al., 2014), pH is here buffered at 7, Equations 1 and 2 reduce to:

$$R_j = S_{V,j} k'_j \omega^{-1} \left(1 - \frac{Q_j}{K_{sp,j}}\right), \quad \text{with } k'_j = k_j |^{25^\circ C} e^{-\frac{E_{a,j}}{RT}} \quad (4)$$

and

$$R_{ox} = k_{ox}^* C_{Fe^{2+}} C_{O_2}, \quad \text{with } k_{ox}^* = \frac{k_{ox}}{K_H} C_{OH-}^2 \quad (5)$$

## 2.3. Reduced Analytical Model

At any moment, the change on DO and dissolved iron concentrations with respect to travel time (τ) is the result of the iron release from Fe(II)-bearing minerals dissolution $R_d$, and iron retention processes, resulting from Fe²⁺ oxidation $R_{ox}$ and precipitation of clay, $R_p$. For the purpose of simplicity, we assume that there is only one species of Fe(II)-bearing mineral that dissolves and one species that precipitates:

$$\begin{cases} \frac{dC_{Fe^{2+}}}{d\tau} = v_d R_d - v_p R_p - R_{ox} \\ \frac{dC_{DO}}{d\tau} = -\frac{1}{4} R_{ox} \end{cases} \quad (6)$$

where $v_d$ and $v_p$ correspond to stoichiometric coefficients accounting for the number of Fe²⁺ per mineral formula in the dissolving or precipitating minerals respectively. $R_d$ and $R_p$ are described by Equation 4 and $R_{ox}$ by Equation 5. In order to solve the system of Equation 6, we consider two regimes: (a) initially oxic conditions, (b) anoxic conditions once oxygen has been depleted.
<table><caption>**Table 1**

Glossary of Main Variables in the Text</caption><thead><tr><th>Nomenclature</th><th>Units</th><th>Variable</th></tr></thead><tbody><tr><td>a</td><td>[mol.kg<sup>-1</sup>]</td><td>Activity</td></tr><tr><td>C</td><td>[mol.kg<sup>-1</sup>]</td><td>Concentration</td></tr><tr><td>Da</td><td></td><td>Damköhler number</td></tr><tr><td>DO</td><td></td><td>Dissolved Oxygen</td></tr><tr><td>E<sub>a</sub></td><td>[kcal mol<sup>-1</sup>]</td><td>Activation energy</td></tr><tr><td>Fe<sup>2+</sup></td><td></td><td>Dissolved Fe(II)</td></tr><tr><td>Fe - clay</td><td></td><td>Fe-rich clay (secondary mineral)</td></tr><tr><td>Φ</td><td>[-]</td><td>Volume fraction</td></tr><tr><td>γ</td><td>[mol Fe.kg<sup>-1</sup>.s]</td><td>Reducing or precipitation capacity of mineral</td></tr><tr><td>Q</td><td>[-]</td><td>Ionic Activity Product</td></tr><tr><td>Γ</td><td>[-]</td><td>Q for all ionic species in mineral excepting iron</td></tr><tr><td>j</td><td></td><td>Index standing for b: biotite, c: Fe-clay, d: dissolving mineral, p: precipitating mineral</td></tr><tr><td>K<sub>H</sub></td><td>[mol.atm<sup>-1</sup>.kg<sup>-1</sup>]</td><td>Henry's constant</td></tr><tr><td>K<sub>sp</sub></td><td>[-]</td><td>Solubility product for mineral</td></tr><tr><td>K<sub>w</sub></td><td>[mol<sup>2</sup> L<sup>-2</sup>]</td><td>Auto-dissociation constant of water</td></tr><tr><td>k</td><td>[mol m<sup>-2</sup> s<sup>-1</sup>]</td><td>Kinetic constant for mineral dissolution</td></tr><tr><td>k<sub>ox</sub></td><td>[kg<sub>w</sub> mol<sup>-1</sup> s<sup>-1</sup>]</td><td>Kinetic constant for iron oxidation by DO</td></tr><tr><td>Λ</td><td>[-]</td><td>Lithological parameter</td></tr><tr><td>ν</td><td>[-]</td><td>number of Fe<sup>2+</sup> atoms per mineral formula</td></tr><tr><td>R</td><td>[mol.kg<sup>-1</sup>.s<sup>-1</sup>]</td><td>Reaction rate</td></tr><tr><td>s</td><td>[mol.kg<sup>-1</sup>]</td><td>Mineral solubility</td></tr><tr><td>S<sub>M</sub></td><td>[m<sup>2</sup> g<sup>-1</sup>]</td><td>Specific surface area</td></tr><tr><td>S<sub>V</sub></td><td>[m<sup>2</sup>.m<sup>-3</sup>]</td><td>Bulk surface area</td></tr><tr><td>τ</td><td>[y]</td><td>Mean fluid travel time</td></tr><tr><td>τ<sub>c</sub></td><td>[y]</td><td>Characteristic time for DO consumption</td></tr><tr><td>τ<sub>f</sub></td><td>[y]</td><td>Characteristic time for DO transport</td></tr><tr><td>u&#773;<sub>a</sub></td><td>[m y<sup>-1</sup>]</td><td>Apparent vertical velocity</td></tr><tr><td>ω</td><td>[-]</td><td>Porosity</td></tr><tr><td>z<sub>c</sub></td><td>[m]</td><td>Reference depth for DO transport</td></tr></tbody></table>

### 2.3.1. Oxic Regime

At short travel times, DO is in excess with respect to Fe<sup>2+</sup> and pH is close to neutrality. Under this condition, very little iron can persist in solution as it is rapidly oxidized according to Equation 2, which has a large kinetic constant. Due to this solute limitation and considering that primary minerals are more abundant than secondary minerals in unweathered rocks, we assume that $R_p \ll R_d$. Therefore, Equation 6 simplifies to:

$$ \left\{ \begin{aligned} \frac{dC_{Fe^{2+}}}{d\tau} + R_{ox} &= v_d R_d \\ \frac{dC_{DO}}{d\tau} &= -\frac{1}{4} R_{ox} \end{aligned} \right. \qquad (7) $$
As Fe(II)-bearing minerals are highly under-saturated, the saturation state of the mineral ($\frac{Q_j}{K_{sp,j}}$) tends to zero in Equation 4. Rearranging Equation 7 with this assumption yields:

$$ \frac{dC_{Fe^{2+}}}{d\tau} - 4\frac{dC_{DO}}{d\tau} = v_d S_{V,d} k'_d \omega^{-1}. \quad (8) $$

Since the oxidation reaction (consuming DO and Fe²⁺) is much faster than the dissolution reaction (releasing Fe²⁺), dissolved iron concentrations are low and vary little because of this kinetic limitation. Thus, $\frac{dC_{Fe^{2+}}}{d\tau} \ll \frac{dC_{DO}}{d\tau}$, which leads to

$$ \frac{dC_{DO}}{d\tau} = -\frac{1}{4}v_d S_{V,d} k'_d \omega^{-1}. \quad (9) $$

Hence, the decay of oxygen is controlled by the amount of dissolved iron produced from mineral dissolution. We define the reducing capacity of a mineral j ($\gamma_j$) as the flux of iron produced by mineral reaction per unit of volume of fluid:

$$ \gamma_j = v_j S_{V,j} k'_j \omega^{-1}, \quad (10) $$

By integration, considering the reducing capacity ($\gamma_d$) constant, the evolution of DO concentration with the fluid travel time $\tau$ is expressed by:

$$ C_{DO} = C_{DO}(0) \left(1 - \frac{\tau}{\tau_c}\right), \quad \forall \tau < \tau_c, \quad (11) $$

where $\tau_c$ is the characteristic time required to consume DO and reach anoxic conditions,

$$ \tau_c = \frac{4C_{DO}(0)}{\gamma_d} \quad (12) $$

The concentration of dissolved iron $C_{Fe^{2+}}$ respects the kinetic law of Equation 5:

$$ \frac{dC_{DO}}{d\tau} = -\frac{1}{4}k_{ox}^* C_{Fe^{2+}} C_{DO}. \quad (13) $$

Hence,

$$ C_{Fe^{2+}} = -\frac{4}{k_{ox}^* C_{DO}} \frac{dC_{DO}}{d\tau}. \quad (14) $$

Inserting Equations 11 and 9 into the above expression leads to the evolution of iron concentration as a function of travel time in the oxic regime,

$$ C_{Fe^{2+}} = \frac{4}{k_{ox}^*} \left( \frac{1}{\tau_c - \tau} \right), \quad \forall \tau < \tau_c \quad (15) $$

This analytical model thus yields the following solutions for the coupled evolution of DO and Fe²⁺ for the oxic regime.

$$ \begin{cases} C_{DO} = C_{DO}(0) \left(1 - \frac{\tau}{\tau_c}\right), \forall \tau < \tau_c \\ C_{Fe^{2+}} = \frac{4}{k_{ox}^*} \left(\frac{1}{\tau_c - \tau}\right), \forall \tau < \tau_c \end{cases} \quad (16) $$

We introduce the non-dimensional Damköhler number Da (Maher, 2010) as the ratio of the characteristic timescale for DO transport $\tau_f$ over the timescale for oxygen consumption $\tau_c$ (Equation 12):

$$ Da = \frac{\tau_f}{\tau_c}. \quad (17) $$

with $\tau_f$ the characteristic time for DO transport up to a reference depth $z_c$ while flowing at an apparent vertical velocity $\overline{v}_a$:

$$ \tau_f = \frac{z_c}{\overline{v}_a}. \quad (18) $$
Da is thus proportional to the ratio of the reducing capacity ($\gamma_d$) to the apparent vertical velocity ($\bar{v}_a$):

$$Da = \frac{\gamma_d}{\bar{v}_a} \times \frac{z_c}{4C_{DO}(0)}. \qquad (19)$$

The system of equations for DO and Fe²⁺ can be expressed in terms of the Damköhler regime, such as:

$$\begin{cases} C_{DO} = C_{DO}(0)(1 - D a \frac{\tau}{\tau_t}), & \forall \tau < \tau_c \\ C_{Fe^{2+}} = \frac{4D a}{k_{ox}^s} \left( \frac{1}{\tau_t - D a \tau} \right), & \forall \tau < \tau_c \end{cases} \qquad (20)$$

For $Da > 1$, the timescale of DO transport $\tau_t$ is longer than the timescale of DO consumption $\tau_c$. Thus, DO supply is transport-limited and the conditions transition from oxic to anoxic along flow paths. Conversely, for $Da < 1$ the timescale of DO consumption is longer than transport. In this case, DO is not depleted because transport overcomes DO consumption and the system remains oxic.

### 2.3.2. Anoxic Regime

Once the available DO has been depleted (i.e., $\tau \ge \tau_c$), the concentration of Fe²⁺ is no longer limited by oxidation ($R_{ox}$ becomes negligible) and then Equation 6 reduces to:

$$\begin{cases} C_{DO} = 0 \\ \frac{dC_{Fe^{2+}}}{d\tau} = v_d R_d - v_p R_p \end{cases} \qquad (21)$$

Since the dissolution rate of primary silicates is slow (Helgeson et al., 1969) and thermodynamics at low P-T impose a high $K_{ps}$, we assume that under reducing conditions, the dissolving Fe(II)-bearing mineral is still highly under-saturated ($\frac{Q_d}{K_{sp,d}} \ll 1$) and releases elements in solution that increase the saturation of Fe(II)-bearing clay (Fe-clay) which then precipitates ($\frac{Q_p}{K_{sp,p}} \gg 1$). Thus, Fe²⁺ is released from the dissolution of primary mineral and part of it precipitates in secondary minerals. According to the definition of $R_j$ (Equation 4) and considering the above approximations for $\frac{Q_j}{K_{sp,j}}$, Equation 21 simplifies to:

$$\frac{dC_{Fe^{2+}}}{d\tau} = \gamma_d - \gamma_p \frac{Q_p}{K_{sp,p}}, \qquad (22)$$

with both $\gamma_d$ and $\gamma_p$, the reducing and precipitation capacities of respectively, the dissolving Fe(II)-bearing primary mineral and the precipitating Fe-clay, defined by Equation 10.

For a general composition of the mineral Fe-clay, the corresponding solubility reaction may be written as:

$$Fe_{v_p} M g_a Al_b Si_c O_{10}(OH)_{m_1} + nH^+ \xrightarrow{K_{sp,p}} v_p Fe^{2+} + aMg^{2+} + bAl^{3+} + cSiO_2^c + m_2 H_2O \qquad (23)$$

Assuming that activity coefficients are close to 1, the ionic activity product ($Q$) for the above reaction can be defined as follows:

$$Q_p = \frac{C_{Mg}^a C_{Al}^b C_{SiO_2}^c}{10^{-n \cdot pH}} C_{Fe^{2+}}^{v_p} = \Gamma_p C_{Fe^{2+}}^{v_p}, \qquad (24)$$

where $C_i$ is the concentration of species i and $\Gamma_p$ the concentration product $\forall i \ne Fe^{2+}$ among the elements present in the Fe-clay mineral. The solubility product ($K_{sp,p}$) may be reformulated in an analogous way:

$$K_{sp,p} = \frac{(a \cdot s_p)^a (b \cdot s_p)^b (c \cdot s_p)^c}{10^{-n \cdot pH}} (v_p \cdot s_p)^{v_p} = \Gamma_p^s s_p^{v_p}, \qquad (25)$$

where $s_p$ is the solubility of Fe-clay and $\Gamma_p^s$ is the solubility product $\forall i \ne Fe^{2+}$. Replacing the expressions for $Q_p$ and $K_{sp,p}$ in Equation 22 and rearranging gives:

$$-\frac{\Lambda}{\gamma_d} \frac{dC_{Fe^{2+}}}{dt} = \left(\frac{C_{Fe^{2+}}}{s_p}\right)^{v_p} - \Lambda, \qquad (26)$$
with the non-dimensional lithological parameter

$$ \Lambda = \frac{\gamma_d \Gamma_p^s}{\gamma_p \Gamma_p} \qquad (27) $$

which quantifies the relative reducing capacity between primary to secondary Fe(II)-bearing minerals ($\gamma_d/\gamma_p$) and the deviation from chemical equilibrium ($\Gamma_p^s/\Gamma_p$). Owing to the definition of the reducing capacity (Equation 10), the lithological parameter introduced herein synthesizes the competition between rates of precipitation and dissolution of Fe(II)-bearing minerals.

The solution of the differential Equation 26 depends on the stoichiometric coefficient for iron per mineral formula of Fe(II)-bearing clay ($\nu_p$). Here, we assume $\nu_p = 2$, which characterizes a typical Fe(II)-rich clay composition according to Sugimori et al. (2008). Solving Equation 26 leads to the following solutions under anoxic conditions:

$$ \begin{cases} C_{DO} = 0 \\ C_{Fe^{2+}} = s_p \sqrt{\Lambda} \tanh\left(\frac{\gamma_d}{s_p \sqrt{\Lambda}} \tau\right), \forall \tau \ge \tau_c. \end{cases} \qquad (28) $$

For travel times larger than $\tau > 2 \frac{s_p \sqrt{\Lambda}}{\gamma_d}$, the transient term in Equation 28 rapidly tends to 1 and the concentration of dissolved iron approaches a steady-state expressed by:

$$ \lim_{\tau \to \infty} C_{Fe^{2+}} = s_p \sqrt{\Lambda} \qquad (29) $$

Under anoxic conditions, the concentration of dissolved iron in Equation 29 represents a far-from-equilibrium steady-state that is only affected by geological factors (the balance between dissolution and precipitation rates of Fe(II)-bearing minerals), here synthesized through the lithological parameter $\Lambda$ and the solubility of Fe-clay $s_p$.

## 2.4. Fully Resolved Reactive Transport Simulations

### 2.4.1. Base Case Simulation

To test the validity of the approximate analytical expressions derived above, we perform reactive transport simulations using the code CrunchFlow (Steefel et al., 2015). As in the analytical model, here we use a Lagrangian approach to simulate a parcel of water traveling along a flow path, following the kinetics of a batch system. We assume an homogeneous bedrock composition and that mineral dissolution reactions are at steady state. We use an effective travel time formulation in which the simulation time in Crunchflow represents the contact time between groundwater and the minerals in the host rock. We simulate the evolution of DO and Fe$^{2+}$, as a result of reactions between oxic water with a composition typical of recharge groundwater (see Table S1 in Supporting Information S1) and a more realistic set of minerals with respect to the simplified analytical case. We consider a silicate rock lithology containing quartz, feldspars, biotite and muscovite, with an initial porosity ($\omega$) of 1% (see Table S2 in Supporting Information S1 for details on the mineralogical composition). In unweathered silicate rocks, biotite is an ubiquitous Fe(II)-bearing mineral present in high proportions (Bazilevskaya et al., 2013; Hampl et al., 2021; Holbrook et al., 2019; Kim et al., 2017) and it has been pointed out as a key mineral for DO consumption (Spiessl et al., 2008). Unlike other reactive-transport models in which each primary mineral is modeled with an independent specific surface area ($S_M$), for example, Heidari et al. (2017), in this work we used the same $S_M$ for all primary minerals. We assume that this $S_M$ is representative of the texture and the exposed rock-surface in fractures that is then affected by the mineral volume fraction ($\Phi$) to obtain the bulk surface area ($S_V$) in Equation 1.

Precipitation and dissolution of secondary minerals are allowed to simulate incongruent dissolution of silicate phases. The incongruent dissolution of feldspars is simulated with the formation of kaolinite according to mineral stability diagrams (see Figure S1 in Supporting Information S1). Biotite weathering is commonly described as an incongruent dissolution process in which the hydration and progressive replacement of inter-layer cations (i.e., $K^+$) forms a wide range of Fe-bearing clays (Fe-clay) and Fe(III) oxyhydroxides depending on the leaching and redox conditions (Acker & Bricker, 1992; Dideriksen et al., 2010; Hampl et al., 2021; Murakami et al., 2003; Scott & Amonette, 1985; Sequeira Braga et al., 2002; Tullborg et al., 2008). In the present simulations, we consider
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Simulation results for the Base Case scenario, obtained with Crunchflow for a lithological composition with $\gamma_d = 2 \times 10^{-3} \text{mol} \cdot \text{m}_{\text{pm}}^{-3} \cdot \text{y}^{-1}$, $\Lambda = 4 \times 10^{4}$. (a) Evolution of Dissolved Oxygen (DO) and iron concentration as a function of travel time ($\tau$). $C_{\text{DO}}^{\text{ox}}$ and $C_{\text{Fe}}^{\text{ox}}$ correspond to concentrations under oxic conditions, whereas $C_{\text{DO}}^{\text{anox}}$ and $C_{\text{Fe}}^{\text{anox}}$ represent concentrations at anoxic conditions. $\tau_c$ is the characteristic time for DO depletion. (b and c) correspond to actual reaction rates for minerals ($R_{\text{mineral}}[\text{mol} \cdot \text{m}_{\text{pm}}^{3} \cdot \text{s}^{-1}]$) as a function of fluid travel time for primary and secondary minerals, respectively.</figcaption></figure>

the precipitation of goethite since it is a typical Fe-oxyhydroxyde found in oxic weathering fronts (Dideriksen et al., 2010). For anoxic conditions, Sugimori et al. (2008) have documented that Fe-rich corrensite forms as the dominant secondary phase. However, since thermodynamic and kinetic parameters for corrensite are poorly known, we allow the formation of chlorite in our simulation because this mineral has a similar stability field and paragenetic relationships with corrensite (Beaufort et al., 1997).

Figure 2 presents the evolution of DO and Fe²⁺ aqueous concentrations as well as actual reaction rates ($R_{\text{mineral}}$) for minerals with mean fluid travel time ($\tau$), as simulated with the numerical model over a period of 1,500 years. We use this extended simulation time to illustrate the stability of the numerical model in the anoxic regime. Although, we prevent the reader that the purpose of analytical model is to focus on travel times that are relevant for modern groundwater residence times in aquifers, typically <1,000 years (Sprenger et al., 2019), where diffusion can be neglected (Sidborn & Neretnieks, 2008; Trinchero et al., 2018).

DO decreases linearly until it gets depleted and a non-linear transition from oxic to anoxic conditions is observed. The evolution of the reaction rates of mineral species with the mean fluid travel time (Figures 2b and 2c) indicates that iron is being released from biotite during both oxic and anoxic conditions ($R_{\text{biotite}} < 0$ indicates dissolution). Under oxic conditions, the Fe²⁺ concentration is driven by the presence of DO and its oxidation into Goethite ($R_{\text{goethite}} > 0$ indicates precipitation). On the other hand, when $\tau \ge \tau_c$, the Fe²⁺ concentration is driven by the precipitation of Fe(II)-bearing clay as indicated by the change on $R_{\text{chlorite}}$ from dissolution (oxic regime) to precipitation (anoxic regime) conditions.

### 2.4.2. Validation of Analytical Expressions

We compare fully resolved numerical simulations with the reduced analytical solutions derived for the evolution of both DO and Fe²⁺ concentrations under oxic (Equation 16) and anoxic (Equation 28) conditions (Figure 2a). Thermodynamic parameters from Table 2 and mineral compositions are used to calculate the parameters of the analytical solutions. The analytical expressions for both DO and Fe²⁺ accurately capture the fully resolved numerical simulations and the predicted characteristic time $\tau_c$ (Equation 12) coincides with the non-linear oxic-anoxic transition from the full numerical model, which validates the assumptions made to derive the reduced analytical model.

To further evaluate the validity of the analytical approximations, a sensitivity analysis is carried out by varying the main parameters identified from the reduced analytical model. With the full numerical model and the reduced
<caption>**Table 2**</caption>

Geochemical Parameters Used in the Numerical Modeling and to Calculate Non-Dimensional Numbers in the Analytical Solutions

<table><thead><tr><th rowspan="2">Parameter</th><th rowspan="2"></th><th colspan="2">Value</th><th rowspan="2">Source</th></tr><tr><th>biotite<sup>a</sup></th><th>Fe-clay<sup>b</sup></th></tr></thead><tbody><tr><td>Molar volume [cm<sup>3</sup> mol<sup>-1</sup>]</td><td>V<sub>m</sub></td><td>149.65</td><td>215.88</td><td>(1)</td></tr><tr><td>Molecular mass [g mol<sup>-1</sup>]</td><td>MM</td><td>417.3</td><td>713.5</td><td></td></tr><tr><td>Solubility product [-]</td><td>K<sub>sp</sub></td><td>10<sup>41.1</sup></td><td>10<sup>47.6</sup></td><td>(2)</td></tr><tr><td>Solubility @ pH = 7 [mol.kg<sup>-1</sup>]</td><td>s</td><td>10<sup>-3.49</sup></td><td>10<sup>-6.85</sup></td><td></td></tr><tr><td>Activation energy [kcal mol<sup>-1</sup>]</td><td>E<sub>a</sub></td><td>5.26</td><td>21.03</td><td>(3)</td></tr><tr><td>kinetic dissolution constant [mol m<sup>-2</sup> s<sup>-1</sup>]</td><td>k<sub>d</sub></td><td>10<sup>-12.55</sup></td><td>10<sup>-12.52</sup></td><td>(3)</td></tr><tr><td>Aqueous parameters</td><td></td><td></td><td></td><td></td></tr><tr><td>Kinetic oxidation constant [L<sup>2</sup> mol<sup>-2</sup> atm<sup>-1</sup> s<sup>-1</sup>]</td><td>k<sub>ox</sub></td><td></td><td>10<sup>12.12</sup></td><td>(4)</td></tr><tr><td>Modified k<sub>ox</sub> from Equation 5 [L mol<sup>-1</sup> s<sup>-1</sup>]</td><td>k<sup>*</sup><sub>ox</sub></td><td></td><td>10.2</td><td></td></tr><tr><td>Initial DO concentration [mol.kg<sup>-1</sup>]</td><td>C<sub>DO</sub>(0)</td><td></td><td>2.62 &times; 10<sup>-4</sup></td><td></td></tr><tr><td>Henry's constant for O<sub>2</sub> [mol atm<sup>-1</sup> L<sup>-1</sup>]</td><td>K<sub>H</sub></td><td></td><td>10<sup>-2.89</sup></td><td></td></tr><tr><td>Water auto-dissociation constant [mol<sup>2</sup>.L<sup>-2</sup>]</td><td>K<sub>w</sub></td><td></td><td>10<sup>-14</sup></td><td></td></tr><tr><td>pH</td><td></td><td></td><td>7.0</td><td></td></tr></tbody></table>

Note. References in column *Source* correspond to (1) (Robie & Philip, 1962), (2) (Malmstrom et al., 1995), (3) (Palandri & Kharaka, 2004), (4) (Singer & Stumm, 1970).

<sup>a</sup>Values for phlogopite. <sup>b</sup>Based on properties for chlorite.

analytical solutions in Equation 16, we evaluate (a) the influence of the reducing capacity (Equation 10) and of the mean fluid travel time on DO and Fe<sup>2+</sup> concentrations (Figures 3a and 3b) in the oxic regime; (b) the influence of the reducing capacity of biotite on the characteristic time for oxygen consumption τ<sub>c</sub> (Figure 3c); (c) the influence of the clay precipitation capacity on the steady-state concentration of Fe<sup>2+</sup> (Figure 3d). The sensitivity on the reducing capacity of biotite is tested by varying the specific surface area (S<sub>M,d</sub>) and the volume fraction of biotite (Φ<sub>d</sub>), both included in the bulk surface area definition (see Equation 1). The specific surface area S<sub>M</sub> varies over four orders of magnitude, and the volume fraction of biotite in the rock Φ<sub>d</sub> by a factor 20. The sensitivity on the clay precipitation capacity is tested by varying the specific surface area S<sub>M</sub> over 3 orders of magnitude and the initial fraction of clay in the rock Φ<sub>p</sub> over 4 orders of magnitude.

The reduced analytical model is in good agreement with the simulations for the full range of tested parameters and simulation times. As expected, less oxic conditions occur in systems with higher biotite contents (Figures 3a and 3b). Insets in Figures 3a and 3b show that longer mean fluid travel times can counterbalance low reducing capacity (low biotite content, Φ<sub>d</sub> and/or low Specif Surface Area, S<sub>M</sub>), leading to more persistent oxic conditions. This interplay between the reducing capacity and travel time is captured by the characteristic time for oxygen consumption τ<sub>c</sub> (Equation 12) in agreement with simulations (Figure 3c). The characteristic time for oxygen consumption decreases with the reducing capacity (both with Φ<sub>d</sub> and S<sub>M</sub>), leading to a longer persistence of oxic conditions in low reducing capacity rocks. The steady-state Fe<sup>2+</sup> concentration in Figure 3d decreases as the initial clay precipitation capacity (γ<sub>p</sub>) approaches the reducing capacity of biotite (γ<sub>p</sub>) (either due to higher initial Fe(II)-clay content Φ<sub>p</sub>(0) or a higher specific surface area S<sub>M,p</sub>), showing that the precipitation of Fe(II)-clays works as a regulator of dissolved iron at anoxic conditions. This effect explains the various C<sub>Fe2+</sub> levels that are reached under anoxic conditions in Figure 3a, as they relate to different γ<sub>d</sub>/γ<sub>p</sub> ratios.

## 2.5. Hydrological and Geological Drivers for DO and Fe<sup>2+</sup> Concentrations

While the full numerical model allows for a complex and realistic reactive modeling framework, the reduced analytical model provide a framework to understand and quantify the hydrological and geological drivers for the evolution of DO and Fe<sup>2+</sup> in silicate catchments (Figure 4).

The DO concentration predicted by Equation 20 is represented as a function of the fluid travel time (τ) and the reducing capacity (γ<sub>Bt</sub>), the two main parameters expressing the hydrological and geological controls in the
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Sensitivity analysis for the numerical and analytical models using the lithological composition in Table S2 in Supporting Information S1. Variations in $\Phi_d$ or $\Phi_p$ are compensated by changing $\Phi_{quartz}$ to keep porosity in 1%. Subfigures (a) and (b) present sensitivities of, respectively, Dissolved Oxygen and Fe$^{2+}$ concentrations to the Specific Surface Area ($S_{M,d}$). Main plots correspond to concentrations computed at a transit time $\tau = 100$ years under different scenarios of Fe(II)-bearing mineral content ($\Phi_d$), while inserts show different travel time scenarios for the $\Phi_d = 1\%$ case. Analytical solutions are represented by Equations 16 and 28. Subfigure (c) presents the sensitivity of $\tau_c$ to $S_{M,d}$ under the same $\Phi_d$ scenarios (i.e., to the reducing capacity in Equation 12). Subfigure (d) presents the sensitivity of the steady-state Fe$^{2+}$ concentration ($C_{Fe}^*$) to the initial clay precipitation capacity (i.e., to both clay's initial volume fraction $\Phi_p(0)$ and specific surface area $S_{M,p}$). Subscript *d* stands for dissolving mineral (e.g., biotite) and subscript *p* stands for precipitating mineral (e.g., Fe-clay).</figcaption></figure>

Damköhler number (Figure 4). The counteracting effects of fluid travel time and reducing capacity imply that there is a fundamental indeterminacy on the hydrological and geological drivers when considering oxygen alone. A given concentration of oxygen can be obtained by a range of different combinations of hydrological and

<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) Dependence of Dissolved Oxygen (DO) concentrations to geological (x-axis) and hydrological (y-axis) forcings. Multiple combinations of reducing capacity ($\gamma_d$) and characteristic transport time ($\tau_t$) can result in a same prediction of DO concentration (e.g., iso-concentration lines for $C_{DO}$). Stronger geological forcings (higher reducing capacity of a certain lithology as in Equation 10) limit the persistence of DO to very short $\tau_t$ (favor anoxic conditions) and vice versa. Higher hydrological forcings in y-axis (faster groundwater velocities) favor deeper transport of DO for a same characteristic transport time, $\tau_t$ (Equation 30). (b) Dependence of the steady-state iron concentration ($C_{Fe}^*$) to the lithological parameter ($\Lambda$) at anoxic conditions. The lithological parameter is defined in Equation 26.</figcaption></figure>
<table><caption>**Table 3**

*General Characteristics of the Study Site*</caption><thead><tr><th></th><th>Guidel</th><th>Kermadoye</th></tr></thead><tbody><tr><td><strong>Generalities</strong></td><td></td><td></td></tr><tr><td>Latitude</td><td>47°45'16"N</td><td>47°44'50"N</td></tr><tr><td>Longitude</td><td>3°28'51"W</td><td>3°25'38"W</td></tr><tr><td>Altitude (m.a.s.l.)</td><td>15.1</td><td>26.7</td></tr><tr><td>Groundwater flow-regime</td><td>Natural circulation</td><td>Pumped since 1991</td></tr><tr><td><strong>Geology</strong></td><td></td><td></td></tr><tr><td>Dominant bedrock lithology</td><td>Mica-schist, Paragneiss</td><td>Granite</td></tr><tr><td>Age</td><td>Ordovician inf</td><td>Carboniferous</td></tr><tr><td>Mean [min, max] depth to fresh bedrock (m)</td><td>17 [4, 34]</td><td>30 [10, 44]</td></tr><tr><td><strong>Climate</strong></td><td></td><td></td></tr><tr><td>Climate type</td><td>Oceanic</td><td></td></tr><tr><td>Mean annual rainfall (mm)</td><td>924<sup>a</sup></td><td></td></tr><tr><td>Mean annual temperature (°C)</td><td>12.1<sup>a</sup></td><td></td></tr></tbody></table>

geological parameters. In rocks with low reducing capacities, DO concentrations are poorly sensitive to travel time since DO consumption is very slow (Da < 1). For high reducing capacity values, DO concentrations evolve rapidly with travel time due to fast reaction kinetics (Da > 1).

When oxygen is depleted, Fe²⁺ concentrations tend toward a far-from-equilibrium steady-state concentration (Equation 29), which is controlled by the relative abundance of primary and secondary Fe(II)-bearing minerals expressed through the non-dimensional number Λ (Equation 27). As a consequence, Fe²⁺ concentrations under anoxic conditions are mainly driven by the geological context of the subsurface. Therefore, the joint analysis of both DO and Fe²⁺ concentrations gives independent constraints on the potential roles of hydrological and geological processes on the distribution of oxygen in the subsurface, which we discuss using field data in the following section.

## 3. Field Study at the Ploemeur CZO (France)

To evaluate the modeling framework presented above, we test it against extensive field observations of DO and Fe²⁺ concentrations in the subsurface available at the Critical Zone Observatory (CZO) of Ploemeur (France). The Ploemeur CZO includes two subcatchments located at a distance of about 4 km: the Kermadoye site and the Guidel site. Both catchments are characterized by a fractured bedrock with similar lithologies and an oceanic climate but differ significantly in their DO and Fe²⁺ concentrations. It is thus well suited to test the concepts presented above.

### 3.1. Site Presentation

The Ploemeur CZO belongs to the H⁺ hydrogeological network (http://hplus.ore.fr/en/), the French network of critical zone observatories OZCAR (https://www.ozcar-ri.org/fr/) and the e-LTER european infrastructure (https://deims.org/731f3ced-148d-4eb5-aa46-870fa22be713). It is located in the southern part of the Armorican massif in Brittany, France. The region is characterized by the intersection of two main tectonic features: (a) a gently dipping (around 30° to the north) contact zone between a Late Hercynian granite and the surrounding micaschist rock; and (b) a dextral-slip normal fault zone which strike north 20° and dip East 70° (Ruelleu et al., 2010). Both are the main transmissive structures of the fractured-bedrock aquifer, characterized by a relatively large average transmissivity on the order of 10⁻³ m²/s sustained by a well connected fracture network (Jiménez-Martínez et al., 2013; Le Borgne et al., 2007).

The Ploemeur CZO hosts two main catchments, Guidel and Kermadoye (see Table 3). The Kermadoye aquifer has been exploited for drinkable water supply since 1991 at an average pumping rate of 1 Mm³ yr⁻¹. This particularly high production rate is attributed to the presence of the regional contact zone which flows toward the vertical faults (Jiménez-Martínez et al., 2013; Leray et al., 2013; Roques et al., 2016). The Kermadoye catchment is monitored with 22 boreholes, with depths ranging from 50 to 150 m. The majority of them cross the contact zone between

<footnote><sup>a</sup>Data from Lann-Bihoué weather station, averaged from the period 2006–2014.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Geographical location, geological and hydrogeological maps of the Ploemeur Critical Zone Observatory and its two catchments, Guidel and Kermadoye. The geological map is reproduced from Béchennec et al. (2012).</figcaption></figure>

micaschist and granite or are entirely in the granite. The Guidel catchment is not pumped, although it has similar hydraulic properties as the Kermadoye catchment. Natural groundwater flows converge to supply a stream and a wetland. The Guidel catchment is located to the north of the contact zone (see Figure 5). It has been monitored since 2009 with 25 boreholes of depths ranging from 50 to 150 m. These boreholes intersect mostly micaschists. Recent studies have shown that the mixing of oxygen rich and iron rich fluids at fracture intersections or in the wetland, sustains microbial hotspots, dominated by iron-oxidizing bacteria (FeOB) (Bethencourt et al., 2020; Bochet et al., 2020). Understanding and modeling DO and Fe²⁺ concentrations at this site is thus of particular interest.

The comparison of the two sites is particularly interesting as the hydrological forcing is stronger at the Kermadoye site due to pumping. Furthermore, while the two sites have comparable geology, the contact zone is deeper at Guidel site, which may be more influenced by the micaschist lithology. Since the two sites have contrasted DO and Fe²⁺ signatures, they are particularly relevant test grounds to resolve the hydrological and geological controls on these chemical properties.

## 3.2. DO and Fe²⁺ Depth Profiles

The DO and iron concentrations measured in the two catchments are represented as a function of depth in Figure 6. Each point represents the concentration measured in front of the main transmissive fracture in each borehole (see Texts S1 and S2 in Supporting Information S1 for Fe²⁺ and DO measurement methods, respectively). Following the methodology proposed by Chatelier et al. (2011), we converted the fracture depth to an effective depth (depth proxy) representative of larger scale flow path. To do so, we used temperature anomalies measured in boreholes and the depth of origin of water (shallow vs. deep) was derived by comparison with the geothermal gradient (details on the depth conversion methodology are presented in Text S3 in Supporting Information S1). DO concentrations decrease with depth, from values close to saturation down to values close to the detection limit (0.1 mg L⁻¹). Both catchments reach the oxic-anoxic transition, but at significantly different depths. The oxic-anoxic transition takes place between 100 and 150 m depth at the Guidel site. At the Kermadoye site it takes place between 200 and 400 m depth. The decrease in DO is associated with an opposite gradient of Fe²⁺ (see Figure 6b). Fe²⁺ concentrations in Guidel increase up to an average value of 2.8 ± 1.7 ppm, while in Kermadoye, Fe²⁺ concentrations increase up to 6.3 × 10⁻² ± 2.7 × 10⁻² ppm, that is, about two orders of magnitude lower than in Guidel.

Both catchments show apparent differences between recharge and discharge zones: high DO and low Fe²⁺ concentrations relate to recharge groundwater while discharge zones are characterized by anoxic groundwater and high Fe²⁺ concentrations.

## 3.3. Model Application to Field Data

To compare field data to our model, we transform travel time to depth based on an apparent vertical velocity ($\overline{v_a}$) that can be constrained from age tracer data (Ayraud et al., 2008):

$$ \tau = \frac{z}{\overline{v_a}} \qquad (30) $$
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Distribution of Dissolved Oxygen (DO) and total iron ($Fe_T$) concentrations as a function of depth for the two catchments. Measured field data is represented with dots and triangles for sampling points from discharge or recharge zones in the catchments. In the oxic part, analytical solutions are drawn as curves with different values of Da for DO and $Fe^{2+}$ (Equation 31). The blue box in the right-hand side plot highlights iron data from oxic boreholes that is not explained by the model (see Section 3.4 for details). In the anoxic part, analytical solutions for $Fe^{2+}$ are drawn with different values of $\Lambda$ (Equation 32). Sample depth has been inferred using temperature as a proxy (see Section 3.2 for details). OAT corresponds to the Oxic-Anoxic Transition ranges.</figcaption></figure>

Thus, the system of equations describing DO and $Fe^{2+}$ concentration profiles can be expressed as a function of depth:

$$ \text{oxic :} \quad \begin{cases} C_{\text{DO}}(z) = C_{\text{DO}}(0) \left(1 - D a \frac{z}{z_c}\right), & \forall z < \frac{z_c}{D a} \\ C_{\text{Fe}^{2+}}(z) = \frac{4 v_a}{k_{\text{ox}}^*} \left(\frac{1}{z_c - D a} z\right), & \forall z < \frac{z_c}{D a} \end{cases} \qquad (31) $$

$$ \text{anoxic :} \quad \begin{cases} C_{\text{DO}}(z) = 0 \\ C_{\text{Fe}}(z) = s_p \sqrt{\Lambda} \tanh \left( \frac{4 C_{\text{DO}}(0) D a}{s_p \sqrt{\Lambda} \frac{z}{z_c}} \right), & \forall z \ge \frac{z_c}{D a} \end{cases} \qquad (32) $$

with $z_c$ the reference depth in Equation 18. Here, we define $z_c = 100$ m, which corresponds to the average depth of our measurements (Figure 6). The dissolved iron concentration $C_{Fe}$ only depends on both Da and $\Lambda$ for $z < \frac{z_c s_p \sqrt{\Lambda}}{4 C_{\text{DO}}(0) D a}$. At deeper depths, $C_{Fe}$ only depends on the lithological parameter and Equation 32 approaches a steady-state concentration (Equation 29).

### 3.4. Comparison of Field Data and Model Simulations of DO and $Fe^{2+}$

The reduced analytical model for DO and $Fe^{2+}$ for the oxic and anoxic regimes (Equations 20 and 29) are fitted to the measured DO and total soluble iron ($Fe_T$) depth-profiles at the two sites. Based on previous works in the study site (Bochet et al., 2020), discrepancies between $Fe_T$ and $Fe^{2+}$ are negligible for the purpose of our study ($Fe^{2+}$/ $Fe_T$ ratio ranges from 93% to 98% in groundwater).

To fit the model to the data, we manually adjusted both dimensionless numbers Da and $\Lambda$ while the rest of the parameters were fixed to values given in Table 2. The apparent vertical velocity ($\overline{v}_a$), that influences the $Fe^{2+}$ concentration under oxic conditions (Equation 31), was fixed at 3 m·y⁻¹ as assessed by Ayraud et al. (2008) for the study site. Under oxic conditions, the solution for DO has one parameter, the Damköhler number Da, and the solution for $Fe^{2+}$ depends on both Da and on the apparent vertical velocity ($\overline{v}_a$) (Equation 20). Theoretically, the
combined resolution of the two equations for DO and Fe²⁺ should allow to fit Da and $\bar{v}_a$, and thus to estimate the reducing capacity of the bedrock (γ).

The reduced analytical model captures the approximate linear decrease in DO with depth in both catchments (Figure 6). The reduced model is also compatible with the evolution of Fe<sub>T</sub> on both sites. For the Guidel site this evolution is well defined over about three orders of magnitude in concentration. For the Kermadoye site, the model describes well the data at anoxic conditions, however, some discrepancies appear between the model and the field data from shallow oxic samples (see blue box in Figure 6). Indeed, the model predicts iron concentrations that significantly lower than those measured from field samples. We hypothesize that these discrepancies are due to a large Fe(III)/Fe<sub>T</sub> ratio in the field measurements. This can be either due to pH effects (shallow waters are normally slightly acid which favorise the solubility of Fe(III) or to the formation of Fe(III)-colloids, particularly if organic matter is available (Nordstrom, 2011; Serikov et al., 2009; Wolthoorn et al., 2004).

The range of Damköhler numbers can be relatively well constrained as consistent values of Da explain jointly the linear decrease of DO and the non-linear increase of Fe²⁺. The Da values range between 0.7 and 1.7 at the Guidel site and between 0.3 and 0.5 at the Kermadoye site (Figure 6). On the other hand, the value of $\bar{v}_a$ influences the very low concentrations of Fe²⁺ at shallow depths (see Figure S4 in Supporting Information S1). However, the discrepancy between measured and modeled Fe²⁺ data at oxic conditions (blue box in Figure 6) prevent us from directly constraining $\bar{v}_a$ from these data. Under anoxic conditions, Fe²⁺ concentrations show a transient phase that depends on both Da and Λ (Equation 32) and that rapidly converges to the steady-state Fe²⁺ concentration from Equation 29. The later depends only on the lithological properties (i.e., Λ). The parameter Λ is relatively well constrained since the maximum Fe²⁺ is very different between the two sites: Λ ranges between 10¹ and 10² for the Kermadoye site, and between 10⁵ and 10⁶ for the Guidel site (Figure 6). The large difference in the two values of Λ suggests that the contrast of the main lithological properties of the two sites is likely more important than possible differences in hydrological properties.

Therefore, the analysis of DO and Fe²⁺ concentration data under oxic conditions provides a good constrain on the Damköhler number, suggesting that the Guidel catchment is more reaction-limited than the Kermadoye catchment. However, this still leaves an uncertainty on whether this difference may be explained by a hydrological or geological contrast. An additional constrain is given when analyzing Fe²⁺ concentrations in the anoxic regime, suggesting an important lithological difference between the two catchments.

# 4. Discussion

## 4.1. Discriminating Hydrological and Geological Controls From Field Data

Despite geographical, climatic and lithological similarities between the two sites, the DO and Fe²⁺ depth-profiles are significantly different. According to our model, this difference may be explained by a factor 3 in the Damköhler numbers Da, the ratio of the transport to reaction time scales, and a factor of 10⁴ in the Λ number, which represents geochemical properties such as the reducing capacity of Fe(II)-bearing minerals and the deviation from chemical equilibrium for clay precipitation. The difference in Damköhler numbers can be either due to a contrast in transport velocity or dissolution rates; or by both effects acting at the same time. We discuss two end-member hypotheses in the following.

We first consider the hypothesis that would attribute DO and Fe²⁺ differences to a contrast mostly caused by hydrological properties. This may be plausible since the Kermadoye site is pumped while the Guidel site is not. Under this hypothesis, the smaller Da for the Kermadoye site would indicate a factor three increase in average groundwater velocity compared to the Guidel site. In such case, travel times should be shorter in the pumped catchment and the apparent age distribution would be expected to have a stronger contribution of young groundwater. However, a geochemical analysis at the Kermadoye site has shown that pumping tends to increase travel times due to the increase of the contribution of deep groundwater flow paths (Roques et al., 2018). This was confirmed by a hydrogeological model of the site (Leray et al., 2014). Therefore, it is unlikely that flow acceleration due to pumping would explain the factor three contrast in Damköhler number.

We now consider the second hypothesis stating that the difference in Da would be mostly caused by a 3-fold contrast in the reducing capacity of the host rock (γ<sub>d</sub>). The Guidel site is composed mostly of micaschists, while the Kermadoye site is composed of both micaschists and granites (Figure 5). Mineralogical analyses carried on
the two rocks indicate that micaschists contain about 7 times more biotite than granite. Therefore, assuming that about 50% of groundwater flowpaths were interacting with granites and 50% with micaschists, the average reducing capacity of the Kermadoye aquifer would be around 3 times lower than the reducing capacity of the Guidel site. This hypothesis would therefore explain a difference of 3 in the Da of the two catchments. The difference in the lithology of both catchments is also supported by the large difference in the parameter Λ, that is entirely related to geochemical parameters (Equation 27). Higher inherited clay contents in granite may account for the 4 orders of magnitude difference in Λ between the two sites. According to Goldich (1938), granites have higher proportions of plagioclase, which easily weathers into clay. This is further reinforced by the occurrence of kaolins near the granite's edge in the location (Béchennec et al., 2012). However, the large difference in Λ could also be related to mineral surface areas which are hard to constrain and vary by orders of magnitude in rocks (Ackerer et al., 2020; Wild et al., 2019).

In summary, the key element that motivates our interpretation is that the difference in Damköhler numbers is of a factor 3, while the difference in Λ are in the order of 10⁴. An assumption of similar γp for the two catchments, would imply a 10⁴ contrast in γd in order to explain the Λ differences. However, since the Damköhler number also depends on γd, we should also observe a difference of the same order for Da, unless it is compensated by a 10⁴ contrast in the velocities, which seems unlikely. Following this reasoning, our hypothesis that the order of magnitude difference in Λ is mostly related to the clay content, that is, to γp rather than to γd.

The analysis of DO and Fe²⁺ data with the presented modeling framework provides key constraints on the hydrological and geological drivers of reactive transport processes. Although the two considered sites are a priori very similar in terms of hydrological and geological contexts, our findings suggest that differences in the proportion of granite and micaschist lead to a strong contrast in DO penetration with depth and Fe²⁺ concentrations. Moreover, we provide new constraints on the hydrogeological functioning of the sites. While previous studies suggested that groundwater was recharged and transported through the micaschist and deeper collected at the contact zone between granite and micaschist (Leray et al., 2013) at the Kermadoye site, this analysis suggests here that there is also a significant proportion of groundwater flowing through the granite.

## 4.2. Controls of DO and Iron Concentrations in Silicate Catchments

The depth of the oxic-anoxic transition in silicate catchments is controlled by the relative importance of the reducing capacity of the rocks and the effective transport of oxygen from the surface, quantified here by the Damköhler number (Figure 6a). Consistently with our model, linear trends in DO with depth have been seen in various hard-rock systems, including the Clara mine in Germany (Bucher et al., 2009) or the Western Canadian Sedimentary Basin (Ruff et al., 2023).

The reducing capacity can be considered as a value inherited from the geology, that varies slowly in time. For instance, Macquarrie et al. (2010) showed a decrease of 0.2% of biotite content in 2 ky under oxic conditions. On the other hand, groundwater table fluctuations or perturbations of the flow regime (e.g., pumping) can be very rapid (timescales ranging from days to the season (Guillaumot et al., 2022; Jiménez-Martínez et al., 2013; Molenat et al., 1999)). Thus, temporal changes in the oxic-anoxic depth in a particular system are likely mostly due to hydrological fluctuations.

The steady-state Fe²⁺ concentration at anoxic conditions is determined by the Lithological Parameter (Λ) that controls the balance between dissolution and precipitation of primary and secondary minerals, respectively. Here, the rationale to introduce the Λ parameter considers a silicate as secondary mineral but the conceptual framework still holds with other types of secondary minerals such as carbonates (see Text S6 in Supporting Information S1 for further details). As described by Equation 32, the steady-state Fe²⁺ concentration depends on the simplified assumption that after DO depletion, no other oxidants can oxidize Fe²⁺. While this model works for our case study, it is possible that in hydrogeological systems the presence of alternative oxidizing agents such as Nitrates or Mn(IV) can promote iron oxidation under anoxic conditions (Kappler et al., 2021). In such case, the predicted rise of Fe²⁺ concentrations up to the steady-state iron concentration would be shifted deeper until the depletion of the alternative oxidizing species.

## 4.3. Model Limitations and Uncertainty in Parameter Estimations

The numerical and analytical models developed in this study are intended to capture the first-order trends in the catchment-scale distributions of DO and Fe²⁺. We therefore introduced several simplifications that could be lifted
in future studies. The dominance of advective transport over matrix diffusion is justified here by the considered relatively short time scales (10 to hundreds of years), consistent with catchment travel times. However, at longer time scales, matrix diffusion would act as an additional buffer for oxygen (Macquarrie et al., 2010). The assumption of negligible mixing among streamlines is based on the discrete nature of flow pathways in the considered crystalline rocks, where flow paths of different travel times only meet in few localized mixing zones (Bochet et al., 2020). Yet, the effect of mixing on DO and Fe dynamics is an interesting question that should be addressed in future studies.

We did not take into account here the spatial heterogeneity in the lithological (constant reducing capacities, $\gamma$) or hydraulic properties (spatially uniform vertical velocity, $\overline{v}_a$). Since our model omits the water-rock interactions in the soil compartment, we used as initial condition for the water chemistry the composition of a typical recharge water instead of rainwater. This choice allows us to focus on the first order impact of other important controlling hydrological and geological factors in the saturated zone without introducing significant bias in the model (see Text S5 in Supporting Information S1 for further details). Similarly, we did not take into account the contribution of microbiological communities in fluid-mineral reactions or the oxidation of structural iron in minerals, since those processes are still poorly constrained. Although, both processes are expected to increase the reducing capacity of the media and have large scale impacts (Buss et al., 2008; Liao et al., 2022; Napieralski et al., 2019). Thus, they should be addressed in future works.

Despite of the above simplifications, the reduced model proposed here provides a simple yet informative tool to identify the processes controlling the evolution of DO and $Fe^{2+}$ with depth. To this end, our model relies on the representation of spatial heterogeneities through effective parameters. This is a common approach in large-scale hydrological and geological models that allows to overcome the still open question of scaling of physical and biogeochemical properties in earth sciences (Li et al., 2017).

As a consequence of this effective modeling approach, the model parameters cannot be fully determined beforehand, and need to be adjusted to field data. This uncertainty in parameter estimation is inherent to reactive transport modeling at the field scale (White & Peterson, 1990), where effective parameters such as the specific surface need to be fitted to the data (Bao et al., 2017; Heidari et al., 2017; Moore et al., 2012). Beyond the precise estimation of these parameters, the interest of the reduced analytical model is to demonstrate that the two sites have different Damköhler and Lambda numbers, which provides insights into the dominant reactive transport processes occurring at the catchment scale.

### 4.4. Persistence of Dissolved Oxygen in the Subsurface

Subsurface environments for which the timescales of DO transport are shorter than timescales of DO consumption (i.e., transport-limited regime, $Da > 1$) favor the deep transport of DO in silicate catchments, which has a major impact on biogeochemical processes in the critical zone. The occurrence of DO in bedrock is responsible for deep WIF (Buss et al., 2008; Kim et al., 2017; Liao et al., 2022), which has been observed in 100 m-deep rock cores (Antoniellini et al., 2017; Bazilevskaya et al., 2013; Dideriksen et al., 2010; Hampl et al., 2021; Holbrook et al., 2019; Krone et al., 2021). By quantifying the dynamics of DO and the parameters controlling the distribution of DO, we provide controls on the conditions favorable for the active oxidative weathering of the deep subsurface. For instance, Bazilevskaya et al. (2013) attributed the thicker regolith in felsic rocks (e.g., granites) compared to mafic rocks (e.g., micaschists) to higher degree of fracturing and higher advective transport of DO with groundwater. We argue that, besides differences in advective transport, the lower reducing capacity of felsic rocks could also explain deeper penetration of DO in the subsurface.

Furthermore, DO exerts a first control on subsurface microbial processes, such as aerobic respiration or denitrification (Dalsgaard et al., 2014; Kolbe et al., 2019) and therefore structures the habitability of subsurface ecosystems. The respiration of Fe-oxidizing bacteria (FeOB) is of particular interest because of the ubiquity of both iron and FeOB (Kappler et al., 2021; Melton et al., 2014) as well as the coupling of the iron cycle with biogeochemical cycles of carbon, sulfur and nitrogen (Casar et al., 2021). The activity of FeOB at near-neutral pH environments requires the simultaneous presence of $Fe^{2+}$ and microaerobic DO concentrations (Druschel et al., 2008; Maisch et al., 2019) while $Fe^{2+}$ and DO often do not coexist because of the rapid oxidation of $Fe^{2+}$. From observations at the Guidel site, Bochet et al. (2020) suggested that the formation of subsurface FeOB hot-spots is favored by the mixing of oxic recharge water with deep anoxic iron-rich water at fracture intersections. As discussed above,
the steady-state Fe²⁺ concentration (Equation 29) is much lower for granite than for micaschist (Figure 6-b). Therefore, even if oxic and anoxic groundwater mix in granite system, the low dissolved iron concentrations do not favor the formation of FeOB hot-spots because Fe²⁺ is limiting. The depth of formation of subsurface FeOB hot-spots therefore not only depends on the transport of DO but also on the availability of Fe²⁺, which is a function of the geological context.

## 5. Conclusions

In this study, we developed a reduced model to describe the depth-distributions of dissolved O₂ and Fe²⁺ concentrations in silicate catchments, which play a central role in a large range of biogeochemical processes. We derived a reduced analytical model, validated with reactive transport simulations, that quantify the parameters controlling jointly DO and Fe²⁺ evolution in fracture zones. Under oxic conditions, DO concentrations decrease linearly with fluid travel time following a slope that is function of the reducing capacity of the bedrock. In this regime, dissolved Fe²⁺ remains low because its aqueous oxidation by DO is faster than its release by mineral dissolution. At the Oxic-Anoxic Transition, DO is depleted and Fe²⁺ concentrations show a rapid non-linear increase up to a far-from-equilibrium steady-state concentration. This concentration is controlled by the relative abundance of primary to secondary Fe(II)-bearing minerals. These reactive transport dynamics can be understood with two non-dimensional parameters: the Damköhler number Da and the lithological parameter Λ.

We use this framework to interpret DO and Fe²⁺ concentrations measured extensively over two neighboring silicate catchments with similar hydrogeological properties but contrasted chemical properties. The differences in the depth of the oxic-anoxic transition and in the Fe²⁺ steady-state concentration are successfully modeled and explained by differences in Da and Λ. The interpretation of DO alone leads to a fundamental indeterminacy in the respective role of geological and hydrological properties that may explain the difference in Damköhler numbers. However, the joint investigation of DO and Fe²⁺ provides additional constraints and points to the role of a geological contrast, here likely due to a difference in the relative proportion of granite and micaschist in the two sites.

The methodology presented here may be implemented on other sites and contexts, to understand and model the depth of the oxic-anoxic transition. The two non-dimensional numbers can be estimated from field data as a guide for DO transport in the subsurface. Here, we investigated the oxidation of Fe²⁺ by DO and assumed that no other oxidants could oxidize Fe²⁺ after DO depletion. However, the presence of alternative oxidizing agents such as Nitrates or Mn(IV) could promote further iron oxidation under anoxic conditions. In this case, the predicted rise of Fe²⁺ concentrations up to the steady-state concentration (C*<sub>Fe²⁺</sub>) would be shifted deeper until the depletion of all oxidizing species. The derivation of the analytical solutions can also be adapted to other reactions involving dissolved reactants and minerals, as long as a silicate mineral is the dominant control on DO consumption and flow occurs through fractures. While here we considered a simplified approach based on an effective travel time, it would be interesting to investigate the form of the analytical solutions and the corresponding dimensionless parameters when representing explicitly structural heterogeneities and fracture-matrix exchanges.

## Data Availability Statement

Numerical simulation were run with the Crunchflow code (Steefel et al., 2015) and figures were made using Matplotlib version 3.5.1 (Hunter, 2007), available under the Matplotlib license at https:// matplotlib.org/. Data from the Ploemeur fractured rock observatory, which is part of the French network of hydrogeological observatories H⁺ (http://hplus.ore.fr/en/), was used to create this manuscript. All field data as well as Crunchflow input files and databases to build upon this research are available in the H⁺ database by following the permanent DOI: https://doi.org/10.26169/hplus.ploemeur_field_data_dissolved_oxygen_and_iron_landscape s.

## References

* Acker, J. G., & Bricker, O. P. (1992). The influence of pH on biotite dissolution and alteration kinetics at low temperature. *Geochimica et Cosmochimica Acta*, 56(8), 3073–3092. https://doi.org/10.1016/0016-7037(92)90290-Y

* Ackerer, J., Jeannot, B., Delay, F., Weill, S., Lucas, Y., Fritz, B., et al. (2020). Crossing hydrological and geochemical modeling to understand the spatiotemporal variability of water chemistry in a headwater catchment (Strengbach, France). *Hydrology and Earth System Sciences*, 24(6), 3111–3133. https://doi.org/10.5194/hess-24-3111-2020

* Antoniellini, M., Mollema, P., & Del Sole, L. (2017). Application of analytical diffusion models to outcrop observations: Implications for mass transport by fluid flow through fractures. *Water Resources Research*, 53(7), 5545–5566. https://doi.org/10.1111/j.1752-1688.1969.tb04897.x

## Acknowledgments

This research was funded, in whole, by ANR EQUIPEX CRITEX project (ANR-11-EQPX-0011) and the ERC project ReactiveFronts (648377). We thank the French Ministry of Higher Education, Research and Innovation and the ANR IRONSTONE (ANR-21-CE01-0008) for funding. A CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grant's open access conditions. We are very grateful to the french networks of hydrogeological sites H⁺ (https://hplus.ore.fr/en/) and of critical zone observatories OZCAR (https://www.ozcar-ri.org/) for providing access to the Ploemeur CZO. We also thank the technical services of Geosciences Rennes CONDATEAU and GeOHeLiS for their support with dissolved gases and groundwater chemistry analysis. The authors would also like to extend their appreciation to the anonymous reviewers for their insightful comments and suggestions, which significantly improved the quality of this manuscript.
Aquilina, L., Roques, C., Boisson, A., Vergnaud-Ayraud, V., Labasque, T., Pauwels, H., et al. (2018). Autotrophic denitrification supported by biotite dissolution in crystalline aquifers (1): New insights from short-term batch experiments. *Science of the Total Environment*, **619–620**, 842–853. https://doi.org/10.1016/j.scitotenv.2017.11.079

Ayraud, V., Aquilina, L., Labasque, T., Pauwels, H., Molenat, J., Pierson-wickmann, A.-c., et al. (2008). Compartmentalization of physical and chemical properties in hard-rock aquifers deduced from chemical and groundwater age analyses. *Applied Geochemistry*, **23**(June), 2686–2707. https://doi.org/10.1016/j.apgeochem.2008.06.001

Bao, C., Li, L., Shi, Y., & Duffy, C. (2017). Understanding watershed hydrogeochemistry: 1. Development of RT-flux-PIHM. *Journal of the American Water Resources Association*, **53**(3), 2328–2345. https://doi.org/10.1111/j.1752-1688.1969.tb04897.x

Bar-on, Y. M., Phillips, R., & Milo, R. (2018). The biomass distribution on Earth. *Proceedings of the National Academy of Sciences of the United States of America*, **115**(25), 6506–6511. https://doi.org/10.1073/pnas.1711842115

Bazilevskaya, E., Lebedeva, M., Pavich, M., Rother, G., Parkinson, D. Y., Cole, D., & Brantley, S. L. (2013). Where fast weathering creates thin regolith and slow weathering creates thick regolith. *Earth Surface Processes and Landforms*, **38**(8), 847–858. https://doi.org/10.1002/esp.3369

Beaufort, D., Baronnet, A., Lanson, B., & Meunier, A. (1997). Corresidence; A single phase or a mixed-layer phyllosilicate in saponite-to-chlorite conversion series? A case study of Sancerre-Couy deep drill hole (France). *American Mineralogist*, **82**(1–2), 109–124. https://doi.org/10.2138/am-1997-1-213

Béchennec, F., Hallégouët, B., Thiéblemont, D., & Thinon, I. (2012). Notice explicative, Carte géol France (1/50 000), feuille Lorient (383). In *BRGM (Ed.).*, *Orléans: BRGM*.

Behrens, R., Bouchez, J., Schuessler, J. A., Dultz, S., Hewawasam, T., & Von Blanckenburg, F. (2015). Mineralogical transformations set slow weathering rates in low-porosity metamorphic bedrock on mountain slopes in a tropical climate. *Chemical Geology*, **411**, 283–298. https://doi.org/10.1016/j.chemgeo.2015.07.008

Bethencourt, L., Bochet, O., Farasin, J., Aquilina, L., Borgne, T. L., Quaiser, A., et al. (2020). Genome reconstruction reveals distinct assemblages of Gallionellaceae in surface and subsurface redox transition zones. *FEMS Microbiology Ecology*, **96**(5). https://doi.org/10.1093/femsec/fiaa036

Bochet, O., Bethencourt, L., Dufresne, A., Farasin, J., Pédrot, M., Labasque, T., et al. (2020). Iron-oxidizer hotspots formed by intermittent oxic–anoxic fluid mixing in fractured rocks. *Nature Geoscience*, **13**(2), 149–155. https://doi.org/10.1038/s41561-019-0509-1

Bucher, K., Zhu, Y., & Stober, I. (2009). Groundwater in fractured crystalline rocks, the Clara mine, Black Forest (Germany). *International Journal of Earth Sciences*, **98**(7), 1727–1739. https://doi.org/10.1007/s00531-008-0328-x

Buss, H. L., Sak, P. B., Webb, S. M., & Brantley, S. L. (2008). Weathering of the Rio Blanco quartz diorite, Luquillo Mountains, Puerto Rico: Coupling oxidation, dissolution, and fracturing. *Geochimica et Cosmochimica Acta*, **72**(18), 4488–4507. https://doi.org/10.1016/j.gca.2008.06.020

Canfield, D. (2004). Canfield, D. E. The early history of atmospheric oxygen: Homage to Robert M. Garrels. *Annual Review of Earth and Planetary Sciences*, **33**, 1–36. https://doi.org/10.1146/annurev.earth.33.092203.122711

Casar, C. P., Momper, L. M., Kruger, B. R., & Osburn, M. R. (2021). Iron-fueled life in the continental subsurface: Deep mine. *Applied and Environmental Microbiology*, **87**(20). https://doi.org/10.1128/aem.00832-21

Chatelier, M., Ruelleu, S., Bour, O., Porel, G., & Delay, F. (2011). Combined fluid temperature and flow logging for the characterization of hydraulic structure in a fractured karst aquifer. *Journal of Hydrology*, **400**(3–4), 377–386. https://doi.org/10.1016/j.jhydrol.2011.01.051

Dalsgaard, T., Stewart, F. J., Thamdrup, B., De Brabandere, L., Revsbech, N. P., Ulloa, O., et al. (2014). Oxygen at nanomolar levels reversibly suppresses process rates and gene expression in anammox and denitrification in the oxygen minimum zone off northern Chile. *mBio*, **5**(6), 019666–e2014. https://doi.org/10.1128/mBio.01966-14

Deng, H., & Spycher, N. (2019). Modeling reactive transport processes in fractures. *Reviews in Mineralogy and Geochemistry*, **85**(1), 49–74. https://doi.org/10.2138/rmg.2019.85.3

DeSimone, L., Mcmahon, P. B., & Rosen, M. (2014). *The quality of our Nation's waters—Water quality in Principal Aquifers of the United States, 1991-2010 (Tech. Rep.)* (Vol. 1360). U.S. Geological Survey Circular. https://doi.org/10.3133/cir1360

Dideriksen, K., Christiansen, B. C., Frandsen, C., Balic-Zunic, T., Mørup, S., & Stipp, S. L. S. (2010). Paleo-redox boundaries in fractured granite. *Geochimica et Cosmochimica Acta*, **74**(10), 2866–2880. https://doi.org/10.1016/j.gca.2010.02.022

Druschel, G. K., Emerson, D., Sutka, R., Suchecki, P., & Luther, G. W. (2008). Low-oxygen and chemical kinetic constraints on the geochemical niche of neutrophilic iron(II) oxidizing microorganisms. *Geochimica et Cosmochimica Acta*, **72**(14), 3358–3370. https://doi.org/10.1016/j.gca.2008.04.035

Edmunds, W. M., Miles, D., & Cook, J. (1984). A comparative study of sequential redox processes in three British aquifers. *Hydrochemical balances of freshwater systems*, **5**, 55–70. https://doi.org/10.2138/rmg.2019.85.3

Emerson, D., Fleming, E. J., & McBeth, J. M. (2010). Iron-oxidizing bacteria: An environmental and genomic perspective. *Annual Review of Microbiology*, **64**(1), 561–583. https://doi.org/10.1146/annurev.micro.112408.134208

Erable, B., Féron, D., & Bergel, A. (2012). Microbial catalysis of the oxygen reduction reaction for microbial fuel cells: A review. *ChemSusChem*, **5**(6), 975–987. https://doi.org/10.1002/cssc.201100836

Fletcher, R. C., Buss, H. L., & Brantley, S. L. (2006). A spheroidal weathering model coupling porewater chemistry to soil thicknesses during steady-state denudation. *Earth and Planetary Science Letters*, **244**(1–2), 444–457. https://doi.org/10.1016/j.epsl.2006.01.055

Frey, P. A., & Reed, G. H. (2012). The ubiquity of iron. *ACS Chemical Biology*, **7**(9), 1477–1481. https://doi.org/10.1021/cb300323q

Goldich, S. S. (1938). A study in rock-weathering. *The Journal of Geology*, **46**(1), 17–58. https://doi.org/10.1086/624619

Gu, X., Heaney, P. J., Reis, F. D., & Brantley, S. L. (2020). Deep abiotic weathering of pyrite. *Science*, **370**(6515), 370. https://doi.org/10.1126/science.abb8092

Guillaumot, L., Longuevergne, L., Marçais, J., Lavenant, N., & Bour, O. (2022). Frequency domain water table fluctuations reveal recharge in fractured aquifers depends on both intense and seasonal rainfall and unsaturated zone thickness. *Hydrology and Earth System Sciences Discussions*, **1–38**. https://doi.org/10.5194/hess-2022-201

Hampl, F. J., Schiperski, F., Byrne, J. M., Schwerdhelm, C., Kappler, A., Bryce, C., et al. (2021). The role of iron-bearing minerals for the deep weathering of a hydrothermally altered plutonic rock in semi-arid climate (Chilean coastal cordillera). *Chemical Geology*, **604**, 120922. https://doi.org/10.1016/j.chemgeo.2022.120922

Hancock, P. J., Boulton, A. J., & Humphreys, W. F. (2005). Aquifers and hyporheic zones: Towards an ecological understanding of groundwater. *Hydrogeology Journal*, **13**(1), 98–111. https://doi.org/10.1007/s10040-004-0421-6

Hartmann, J., & Moosdorf, N. (2012). The new global lithological map database GLiM: A representation of rock properties at the Earth surface. *Geochemistry, Geophysics, Geosystems*, **13**(12), 1–37. https://doi.org/10.1029/2012GC004370

Heidari, P., Li, L., Jin, L., Williams, J. Z., & Brantley, S. L. (2017). A reactive transport model for Marcellus shale weathering. *Geochimica et Cosmochimica Acta*, **217**, 421–440. https://doi.org/10.1016/j.gca.2017.08.011
Helgeson, H. C., Garrels, R. M., & MacKenzie, F. T. (1969). Evaluation of irreversible reactions in geochemical processes involving minerals and aqueous solutions—II. Applications. *Geochimica et Cosmochimica Acta*, **33**(4), 455-481. https://doi.org/10.1016/0016-7037(69)90127-6

Holbrook, W. S., Marcon, V., Bacon, A. R., Brantley, S. L., Carr, B. J., Flinch, B. A., et al. (2019). Links between physical and chemical weathering inferred from a 65-m-deep borehole through Earth's critical zone. *Scientific Reports*, **9**, 1–11. https://doi.org/10.1038/s41598-019-40819-9

Humphreys, W. F. (2009). Hydrogeology and groundwater ecology: Does each inform the other? *Hydrogeology Journal*, **17**(1), 5–21. https://doi.org/10.1007/s10040-008-0349-3

Hunter, J. D. (2007). Matplotlib: A 2d graphics environment [Software]. Computing in Science & Engineering, **9**, 90–95. https://doi.org/10.1109/MCSE.2007.55

Jiménez-Martínez, J., Longuevergne, L., Le Borgne, T., Davy, P., Russian, A., & Bour, O. (2013). Temporal and spatial scaling of hydraulic response to recharge in fractured aquifers: Insights from a frequency domain analysis. *Water Resources Research*, **49**(5), 3007–3023. https://doi.org/10.1002/wrcr.20260

Kappler, A., Bryce, C., Mansor, M., Lueder, U., Byrne, J. M., & Swanner, E. D. (2021). An evolving view on biogeochemical cycling of iron. *Nature Reviews Microbiology*, **19**(6), 360–374. https://doi.org/10.1038/s41579-020-00502-7

Kendall, B., Anbar, A. D., Kappler, A., & Konhauser, K. O. (2012). The global iron cycle. *Fundamentals of Geobiology*, 65–92. https://doi.org/10.1002/9781118280874.ch6

Kim, H., Stinchcomb, G., & Brantley, S. L. (2017). Feedbacks among O2 and CO2 in deep soil gas, oxidation of ferrous minerals, and fractures: A hypothesis for steady-state regolith thickness. *Earth and Planetary Science Letters*, **460**, 29–40. https://doi.org/10.1016/j.epsl.2016.12.003

Kolbe, T., De Dreuzy, J. R., Abbott, B. W., Aquilina, L., Babey, T., Green, C. T., et al. (2019). Stratification of reactivity determines nitrate removal in groundwater. *Proceedings of the National Academy of Sciences of the United States of America*, **116**(7), 2494–2499. https://doi.org/10.1073/pnas.1816892116

Korom, S. F. (1992). Natural denitrification in the saturated zone: A review. *Water Resources Research*, **28**(6), 1657–1668. https://doi.org/10.1029/92WR00252

Krone, L. V., Hampl, F. J., Schwerdhelm, C., Bryce, C., Ganzert, L., Kitte, A., et al. (2021). Deep weathering in the semi-arid coastal cordillera, Chile. *Scientific Reports*, **11**(1), 1–15. https://doi.org/10.1038/s41598-021-90267-7

Lasaga, A. C. (1984). Chemical kinetics of water-rock interactions. *Journal of Geophysical Research*, **89**(B6), 4009–4025. https://doi.org/10.1029/JB0891B06p04009

Le Borgne, T., Bour, O., Riley, M. S., Gouze, P., Pezard, P. A., Belghoul, A., et al. (2007). Comparison of alternative methodologies for identifying and characterizing preferential flow paths in heterogeneous aquifers. *Journal of Hydrology*, **345**(3–4), 134–148. https://doi.org/10.1016/j.jhydrol.2007.07.007

Leray, S., de Dreuzy, J.-R., Aquilina, L., Vergnaud-Ayraud, V., Labasque, T., Bour, O., & Le Borgne, T. (2014). Temporal evolution of age data under transient pumping conditions. *Journal of Hydrology*, **511**, 555–566. https://doi.org/10.1016/j.jhydrol.2014.01.064

Leray, S., de Dreuzy, J. R., Bour, O., & Bresciani, E. (2013). Numerical modeling of the productivity of vertical to shallowly dipping fractured zones in crystalline rocks. *Journal of Hydrology*, **481**, 64–75. https://doi.org/10.1016/j.jhydrol.2012.12.014

Li, L., Maher, K., Navarre-sitchler, A., Druhan, J., Meile, C., Lawrence, C., et al. (2017). Expanding the role of reactive transport models in critical zone processes. *Earth-Science Reviews*, **165**, 280–301. https://doi.org/10.1016/j.earscirev.2016.09.001

Liao, R., Gu, X., & Brantley, S. L. (2022). Weathering of chlorite from grain to watershed: The role and distribution of oxidation reactions in the subsurface. *Geochimica et Cosmochimica Acta*, **333**, 284–307. https://doi.org/10.1016/j.gca.2022.07.019

Macquarrie, K. T. B., Mayer, K. U., Jin, B., & Spiessl, S. M. (2010). The importance of conceptual models in the reactive transport simulation of oxygen ingress in sparsely fractured crystalline rock. *Journal of Contaminant Hydrology*, **112**(1–4), 64–76. https://doi.org/10.1016/j.jconhyd.2009.10.007

Mader, M., Schmidt, C., van Geldern, R., & Barth, J. A. (2017). Dissolved oxygen in water and its stable isotope effects: A review. *Chemical Geology*, **473**(September), 10–21. https://doi.org/10.1016/j.chemgeo.2017.10.003

Maher, K. (2010). The dependence of chemical weathering rates on fluid residence time. *Earth and Planetary Science Letters*, **294**(1–2), 101–110. https://doi.org/10.1016/j.epsl.2010.03.010

Maher, K., & Chamberlain, C. P. (2014). Hydrologic regulation of chemical weathering and the geologic carbon cycle. *Science*, **343**(6178), 1502–1504. https://doi.org/10.1126/science.1250770

Maisch, M., Lueder, U., Laufer, K., Scholze, C., Kappler, A., & Schmidt, C. (2019). Contribution of Microaerophilic iron(II)-Oxidizers to iron(III) mineral formation. *Environmental Science and Technology*, **53**(14), 8197–8204. https://doi.org/10.1021/acs.est.9b01531

Malard, F., & Hervant, F. (1999). Oxygen supply and the adaptations of animals in groundwater. *Freshwater Biology*, **41**, 1–30. https://doi.org/10.1046/j.1365-2427.1999.00379.x

Malmstrom, M., Banwart, S., Duro, L., Wersin, P., & Bruno, J. (1995). Biotite and chlorite weathering at 25°C (Tech. Rep. No. January).

Malmström, M., Banwart, S., Lewenhagen, J., Duro, L., & Bruno, J. (1996). The dissolution of biotite and chlorite at 25°C in the near-neutral pH region. *Journal of Contaminant Hydrology*, **21**(1–4), 201–213. https://doi.org/10.1016/0169-7722(95)00047-x

Melton, E. D., Swanner, E. D., Behrens, S., Schmidt, C., & Kappler, A. (2014). The interplay of microbially mediated and abiotic reactions in the biogeochemical Fe cycle. *Nature Reviews Microbiology*, **12**(12), 797–808. https://doi.org/10.1038/nrmicro3347

Molenat, J., Davy, P., Gascuel-Odoux, C., & Durand, P. (1999). Study of three subsurface hydrologic systems based on spectral and cross-spectral analysis of time series. *Journal of Hydrology*, **222**(1–4), 152–164. https://doi.org/10.1016/s0022-1694(99)00107-9

Moore, J., Lichtner, P. C., White, A. F., & Brantley, S. L. (2012). Using a reactive transport model to elucidate differences between laboratory and field dissolution rates in regolith. *Geochimica et Cosmochimica Acta*, **93**, 235–261. https://doi.org/10.1016/j.gca.2012.03.021

Murakami, T., Utsunomiya, S., Yokoyama, T., & Kasama, T. (2003). Biotite dissolution processes and mechanisms in the laboratory and in nature: Early stage weathering environment and vermiculitization. *American Mineralogist*, **88**(2–3), 377–386. https://doi.org/10.2138/am-2003-2-314

Napieralski, S. A., Buss, H. L., Brantley, S. L., Lee, S., Xu, H., & Roden, E. E. (2019). Microbial chemolithotrophy mediates oxidative weathering of granitic bedrock. *Proceedings of the National Academy of Sciences of the United States of America*, **116**(52), 26394–26401. https://doi.org/10.1073/pnas.1909970117

Nordstrom, D. K. (2011). Hydrogeochemical processes governing the origin, transport and fate of major and trace elements from mine wastes and mineralized rock to surface waters. *Applied Geochemistry*, **26**(11), 1777–1791. https://doi.org/10.1016/j.apgeochem.2011.06.002

Palandri, J., & Kharaka, Y. (2004). A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modeling.

Pedersen, K. (1997). Microbial life in deep granitic rock. *FEMS Microbiology Reviews*, **20**(3–4), 399–414. https://doi.org/10.1016/S0168-6445(97)00022-3

Petsch, S. T., Bolton, E., Mok, U., & Evans, B. (2004). The weathering of sedimentary organic matter as a control on atmospheric O₂: I. Analysis of a black shale. *American Journal of Science*, **304**(3), 234–249. https://doi.org/10.2475/ajs.304.3.234

Robie, A. R., & Philip, M. B. (1962). *Molar volumes and densities of minerals—TEI-822 (Tech. Rep.).* USGS.
Roques, C., Aquilina, L., Boisson, A., Vergnaud-Ayraud, V., Labasque, T., Longuevergne, L., et al. (2018). Autotrophic denitrification supported by biotite dissolution in crystalline aquifers: (2) transient mixing and denitrification dynamic during long-term pumping. *Science of the Total Environment*, 619–620, 491–503. https://doi.org/10.1016/j.scitotenv.2017.11.104

Roques, C., Bour, O., Aquilina, L., & Dewandel, B. (2016). High-yielding aquifers in crystalline basement: Insights about the role of fault zones, exemplified by Armorican massif, France. *Hydrogeology Journal*, 24(8), 2157–2170. https://doi.org/10.1007/s10040-016-1451-6

Ruelleu, S., Moreau, F., Bour, O., Gapais, D., & Martelet, G. (2010). Impact of gently dipping discontinuities on basement aquifer recharge: An example from Ploemeur (Brittany, France). *Journal of Applied Geophysics*, 70(2), 161–168. https://doi.org/10.1016/j.jappgeo.2009.12.007

Ruff, S. E., Humez, P., de Angelis, I. H., Diao, M., Nightingale, M., Cho, S., et al. (2023). Hydrogen and dark oxygen drive microbial productivity in diverse groundwater ecosystems. *Nature Communications*, 14(1), 3194. https://doi.org/10.1038/s41467-023-38523-4

Scott, A. D., & Amonette, J. (1985). Role of iron in mica weathering. In *Iron in soils and clay minerals* (pp. 537–605). Bad Windsheim: Reidel Publishing Company. https://doi.org/10.1007/978-94-009-4007-9

Seeboonruang, U., & Ginn, T. R. (2006). Upscaling heterogeneity in aquifer reactivity via exposure-time concept: Forward model. *Journal of Contaminant Hydrology*, 84(3–4), 127–154. https://doi.org/10.1016/j.jconhyd.2005.12.011

Sequeira Braga, M. A., Paquet, H., & Begonha, A. (2002). Weathering of granites in a temperate climate (NW Portugal): Granitic saprolites and arenization. *Catena*, 49(1), 41–56. https://doi.org/10.1016/S0341-8162(02)00017-6

Serikov, L. V., Tropina, E. A., Shiyan, L. N., Frimmel, F. H., Metreveli, G., & Delay, M. (2009). Iron oxidation in different types of groundwater of Western Siberia. *Journal of Soils and Sediments*, 9(2), 103–110. https://doi.org/10.1007/s11368-009-0069-x

Sidborn, M., & Neretnieks, I. (2007). Long term redox evolution in granitic rocks: Modelling the redox front propagation in the rock matrix. *Applied Geochemistry*, 22(11), 2381–2396. https://doi.org/10.1016/j.apgeochem.2007.05.007

Sidborn, M., & Neretnieks, I. (2008). Long-term oxygen depletion from infiltrating groundwaters: Model development and application to intra-glaciation and glaciation conditions. *Journal of Contaminant Hydrology*, 100(1–2), 72–89. https://doi.org/10.1016/j.jconhyd.2008.05.010

Singer, P. C., & Stumm, W. (1970). Acidic mine drainage: The rate-determining step. *Science*, 167(3921), 1121–1123. https://doi.org/10.1126/science.167.3921.1121

Singha, K., & Navarre-Sitchler, A. (2022). The importance of groundwater in critical zone science. *Ground Water*, 60(1), 27–34. https://doi.org/10.1111/gwat.13143

Spiessl, S. M., Macquarrie, K. T. B., & Mayer, K. U. (2008). Identification of key parameters controlling dissolved oxygen migration and attenuation in fractured crystalline rocks. *Journal of Contaminant Hydrology*, 95(3–4), 141–153. https://doi.org/10.1016/j.jconhyd.2007.09.002

Sprenger, M., Stummpp, C., Weiler, M., Aeschbach, W., Allen, S. T., Benettin, P., et al. (2019). The demographics of water: A review of water ages in the critical zone. *Reviews of Geophysics*, 57(3), 800–834. https://doi.org/10.1029/2018RG000633

Steefel, C. I., Appelo, C. A. J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., et al. (2015). Reactive transport codes for subsurface environmental simulation [Software]. Computational Geosciences, 19, 445–478. https://doi.org/10.1007/s10596-014-9443-x

Stumm, W., & Morgan, J. J. (1996). Oxidation and reduction: Equilibria and microbial mediation. In J. Schnoor & A. Zehnder (Eds.), *Aquatic chemistry: Chemical equilibria and rates in natural waters* (3rd ed., pp. 425–515). Wiley.

Sugimori, H., Iwatsuki, T., & Murakami, T. (2008). Chlorite and biotite weathering, Fe²⁺-rich corrodine formation, and Fe behavior under low Po₂ conditions and their implication for Precambrian weathering. *American Mineralogist*, 93(7), 1080–1089. https://doi.org/10.2138/am.2008.2663

Sullivan, P. L., Hynek, S. A., Gu, X., Singha, K., White, T., West, N., et al. (2016). Oxidative dissolution under the channel leads geomorphological evolution at the shale hills catchment. *American Journal of Science*, 316(10), 981–1026. https://doi.org/10.2475/10.2016.02

Tebo, B. M., Johnson, H. A., McCarthy, J. K., & Templeton, A. S. (2005). Geomicrobiology of manganese(II) oxidation. *Trends in Microbiology*, 13(9), 421–428. https://doi.org/10.1016/j.tim.2005.07.009

Trinchero, P., Molinero, J., Ebrahimi, H., Puigdomenech, I., Gylling, B., Svensson, U., et al. (2018). Simulating oxygen intrusion into highly performance computing. *Mathematical Geosciences*, 50(5), 49–51. https://doi.org/10.1007/s11004-017-9718-6

Trinchero, P., Puigdomenech, I., Molinero, J., Ebrahimi, H., Gylling, B., Svensson, U., et al. (2017). Continuum-based DFN-consistent numerical framework for the simulation of oxygen infiltration into fractured crystalline rocks. *Journal of Contaminant Hydrology*, 200, 60–69. https://doi.org/10.1016/j.jconhyd.2017.04.001

Trinchero, P., Sidborn, M., Puigdomenech, I., Svensson, U., Ebrahimi, H., Molinero, J., et al. (2019). Transport of oxygen into granitic rocks: Role of physical and mineralogical heterogeneity. *Journal of Contaminant Hydrology*, 220, 108–118. https://doi.org/10.1016/j.jconhyd.2018.12.001

Tullborg, E. L., Drake, H., & Sandström, B. (2008). Palaeohydrogeology: A methodology based on fracture mineral studies. *Applied Geochemistry*, 23(7), 1881–1897. https://doi.org/10.1016/j.apgeochem.2008.02.009

White, A. F., & Peterson, M. L. (1990). *Role of reactive-surface-area characterization in geochemical kinetic models*. ACS Publications.

White, A. F., & Yee, A. (1985). Aqueous oxidation-reduction kinetics associated with coupled electron-cation transfer from iron-containing silicates at 25°C. *Geochimica et Cosmochimica Acta*, 49(5), 1263–1275. https://doi.org/10.1016/0016-7037(85)90015-8

White, A. F., Yee, A., & Flexser, S. (1985). Surface oxidation-reduction kinetics associated with experimental basalt-water reaction at 25°C. *Chemical Geology*, 49(1–3), 73–86. https://doi.org/10.1016/0009-2541(85)90148-2

Wild, B., Daval, D., Beaulieu, E., Pierret, M. C., Viville, D., & Imfeld, G. (2019). In-situ dissolution rates of silicate minerals and associated bacterial communities in the critical zone (Strengbach catchment, France). *Geochimica et Cosmochimica Acta*, 249, 95–120. https://doi.org/10.1016/j.gca.2019.01.003

Winograd, I., & Robertson, F. (1982). Deep oxygenated ground water: Anomaly or common occurrence? *Science*, 216(4551), 1227–1230. https://doi.org/10.1126/science.216.4551.1227

Wolthoorn, A., Temminghoff, E. J., Weng, L., & Van Riemsdijk, W. H. (2004). Colloid formation in groundwater: Effect of phosphate, manganese, silicate and dissolved organic matter on the dynamic heterogeneous oxidation of ferrous iron. *Applied Geochemistry*, 19(4), 611–622. https://doi.org/10.1016/j.apgeochem.2003.08.003

## References From the Supporting Information

Blanc, P., Lassin, A., Piantone, P., Azaroual, M., Jacquemet, N., Fabbri, A., & Gaucher, E. C. (2012). Thermoderm: A geochemical database focused on low temperature water/rock interactions and waste materials [Dataset]. Applied Geochemistry, 27, 2107–2116. https://doi.org/10.1016/j.apgeochem.2012.06.002

Klepikova, M. V., Le Borgne, T., Bour, O., & Davy, P. (2011). A methodology for using borehole temperature-depth profiles under ambient, single and cross-borehole pumping conditions to estimate fracture hydraulic properties. *Journal of Hydrology*, 407(1–4), 145–152. https://doi.org/10.1016/j.jhydrol.2011.07.018
Le Borgne, T., Paillet, F., Bour, O., & Caudal, J. P. (2006). Cross-borehole flowmeter tests for transient heads in heterogeneous aquifers. *Ground Water, 44*(3), 444–452. https://doi.org/10.1111/j.1745-6584.2005.00150.x

Pouladi, B., Bour, O., Longuevergne, L., de La Bernardie, J., & Simon, N. (2021). Modelling borehole flows from distributed temperature sensing data to monitor groundwater dynamics in fractured media. *Journal of Hydrology, 598*, 126450. https://doi.org/10.1016/j.jhydrol.2021.126450

Sugisaki, R., & Taki, K. (1987). Simplified analyses in of He, natural Ne, waters and Ar dissolved. *Geochemical Journal*, 21(1), 23–27. https://doi.org/10.2343/geochemj.21.23