# Effective reaction kinetics of steady mixing fronts in porous media

Satoshi Izumoto, Gauthier Rousseau, Tanguy Le Borgne, Joris Heyman

► To cite this version:

Satoshi Izumoto, Gauthier Rousseau, Tanguy Le Borgne, Joris Heyman. Effective reaction kinetics of steady mixing fronts in porous media. Journal of Fluid Mechanics, 2025, 1013, pp.A4. 10.1017/jfm.2025.10223. insu-05120842

HAL Id: insu-05120842

https://insu.hal.science/insu-05120842v1

Submitted on 19 Jun 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Effective reaction kinetics of steady mixing fronts in porous media

Satoshi Izumoto¹ ID, Gauthier Rousseau¹ ID, Tanguy Le Borgne¹ ID and Joris Heyman¹ ID

¹Univ. Rennes, CNRS, Géosciences Rennes UMR 6118, Rennes 35000, France

**Corresponding author:** Joris Heyman, joris.heyman@univ-rennes.fr

(Received 29 August 2024; revised 30 April 2025; accepted 2 May 2025)

Mixing-induced reactions play an important role in a wide range of porous media processes. Recent advances have shown that fluid flow through porous media leads to chaotic advection at the pore scale. However, how this impacts Darcy-scale reaction rates is unknown. Here, we measure the reaction rates in steady mixing fronts using a chemiluminescence reaction in index-matched three-dimensional porous media. We consider two common mixing scenarios for reacting species, flowing either in parallel in a uniform flow or towards each other in a converging flow. We study the reactive properties of these fronts for a range of Péclet numbers. In both scenarios, we find that the reaction rates significantly depart from the prediction of hydrodynamic dispersion models, which obey different scaling laws. We attribute this departure to incomplete mixing effects at the pore scale, and propose a mechanistic model describing the pore-scale deformations of the front triggered by chaotic advection and their impact on the reaction kinetics. The model shows good agreement with the effective Darcy-scale reaction kinetics observed in both uniform and converging flows, opening new perspectives for upscaling reactive transport in porous media.

**Key words:** porous media, coupled diffusion and flow, laminar reacting flows

## 1. Introduction

Reactive processes in porous media play a central role in the transport, transformation and degradation of chemical and biological substances in a wide range of systems, including soils and rocks (Lichtner, Steefel & Oelkers 2018), engineered porous media (Perfect et al. 2014; Huang et al. 2022) and biological porous media (Khaled & Vafai 2003;

© The Author(s), 2025. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Goirand, Le Borgne & Lorthois (2021). In many applications, reactions are driven by transport processes that bring dissolved reactants into contact (Rolle & Le Borgne 2019; Valocchi, Bolster & Werth 2019). This reactive transport process is generally modelled using the advection dispersion reaction equation (Dentz et al. 2011), which describes mixing processes using an effective dispersion coefficient that lumps pore-scale velocity fluctuations into a single coefficient. During the last decades, there have been increasing observations of the breakdown of this model due to incomplete mixing of reactants at the pore scale (Raje & Kapoor 2000; Gramling, Harvey & Meigs 2002; Klenk & Grathwohl 2002; Battiato et al. 2009; Battiato & Tartakovsky 2011; De Anna et al. 2014).

Although different phenomenological frameworks have been proposed (Olsson & Grathwohl 2007; Benson & Meerschaert 2008; Edery, Scher & Berkowitz 2010; Sanchez-Vila, Fernández-Garcia & Guadagnini 2010; Chiogna & Bellin 2013; Ding et al. 2013; Hochstetler & Kitanidis 2013; Hochstetler et al. 2013; Porta et al. 2013; Paster, Bolster & Benson 2014; Bolster, Paster & Benson 2016; Ginn 2018), it is still unclear how to correctly upscale reactive processes at the Darcy scale. At the pore scale, the interplay between fluid deformations, molecular diffusion and reaction is conveniently described in a Lagrangian frame using the lamellar theory of mixing (Ranz 1979; Le Borgne, Dentz & Villermaux 2015; Villermaux 2019). The latter represents mixing interfaces as elongated lamellar structures whose stretching under the action of flow enhances mixing and reaction rates (Bandopadhyay et al. 2017; Izumoto et al. 2023). However, it is unknown how to upscale such a pore-scale description to three-dimensional (3-D) Darcy-scale porous media. Furthermore, it is unclear how the mixing dynamics at the pore scale interacts with non-uniform flows at larger scale. A case of particular interest is converging flows in which reactive fluids flow towards each other in porous media flows (Hester et al. 2017; Marzadri et al. 2017; Bresciani, Kang & Lee 2019; de Vriendt 2021; Izumoto et al. 2023). This leads to a stagnation point of the flow associated with persistent fluid compression, which can sustain reactions.

Using two-dimensional (2-D) microfluidic devices, several studies have succeeded in quantifying the evolution of concentrations of reacting species at the pore scale in mixing fronts. Reactions producing light or a fluorescent product have been used to image reaction rates in 2-D micromodels (Willingham, Werth & Valocchi 2008; De Anna et al. 2014; Izumoto et al. 2023), yielding high-resolution maps of local reaction intensity. Although 2-D porous micromodels greatly facilitate imaging at the pore scale, their mixing dynamics fundamentally differ from 3-D porous media flows, which sustain chaotic advection (Lester et al. 2013, 2016). Chaotic advection has recently been shown to control dilution rates and scalar gradients on the pore scale (Heyman et al. 2020), affecting incomplete mixing and reactions (Sanquer et al. 2024) and fluid-solid reactions (Aquino, Le Borgne & Heyman 2023). Experimentally, it is particularly challenging to resolve reactive transport processes at the pore scale in 3-D porous media due to the opacity of the grains, and most studies perform measurements of the reaction rate at the outlet of columns (Ham et al. 2007). Recent studies have resolved mixing-induced reactions on the pore scale using magnetic resonance imaging (Markale et al. 2021), or fluid-grain refractive index matching together with chemical reactions producing fluorescent compounds (Sanquer et al. 2024). However, these high-resolution imaging techniques are currently limited to length scales of a few tens of grains, which is not sufficient to bridge the gap to Darcy scales. Thus the impact of chaotic mixing on the large-scale reaction efficiency remains unclear. Colorimetric reactions have been used in refractive index matching media on the Darcy scale to measure the integrated mass of reaction products in an invading front (Gramling et al. 2002; Edery et al. 2015), demonstrating the persistence of the signature of incomplete mixing at large scales. These experiments were obtained in transient conditions
in which a fluid flows and replaces another fluid. Such longitudinal mixing conditions are expected to be dominated by longitudinal shear rather than chaotic advection, which in steady 3-D flows occurs in the transverse direction of the flow (Lester et al. 2018). Hence the impact of pore-scale chaotic mixing on the Darcy-scale effective reaction reaction front properties remains unknown.

Here, we investigate the properties of the Darcy-scale reaction front of an irreversible bimolecular reaction $A+B \rightarrow C$ in steady transverse mixing fronts through porous media. In our flow cells, reactants were continuously injected by two inlets, and mixing occurred in the direction perpendicular to the mean flow direction. We combined the refractive index matching technique with chemiluminescence to map in two dimensions the effective Darcy-scale reaction rates integrated over the gap of the cell. We considered two common scenarios of mixing fronts in porous media: (i) a co-flow configuration where two reactive fluids flow in parallel in a uniform Darcy-scale flow, and (ii) a saddle flow where two reactive fluids flow towards each other in a converging flow around a stagnation point. The paper is organised as follows. We first present the theory for mixing-induced reactions in parallel and converging flows when considering diffusion or dispersion, thus neglecting the effect of incomplete mixing at the pore scale. We then compare these classical predictions with the experimental measurements, showing a significantly different evolution of the effective reaction rate with distance along the front and Péclet number because of incomplete mixing. Finally, we derive a new mechanistic model linking the pore-scale chaotic mixing to the effective reaction front properties at the Darcy scale, capturing experimental data for the range of investigated Péclet numbers, in both parallel and converging flows.

## 2. Reactive front kinetics under diffusion and dispersion

In this section, we present the theory for bimolecular reactions $A+B \rightarrow C$ in mixing fronts when considering diffusion or hydrodynamic dispersion, for both parallel and converging flows. Deviations of experimental data from this reference model, which ignores incomplete mixing at the pore scale, will be discussed in the following sections to quantify the role of micro-scale mixing dynamics in macro-scale reactive front properties.

Given a flow velocity $v$, the dispersion tensor $D$, and $k$ the kinetic constant of a bimolecular reaction, the transport equation for the Darcy-scale concentration of the reactants $C_A$ and $C_B$ follows:

$$ \frac{\partial C_A}{\partial t} + v \cdot \nabla C_A = \nabla \cdot (D \nabla C_A) - k C_A C_B, \\ \frac{\partial C_B}{\partial t} + v \cdot \nabla C_B = \nabla \cdot (D \nabla C_B) - k C_B C_A. \qquad (2.1) $$

The dispersion tensor lumps the effect of pore-scale velocity fluctuations and molecular diffusion into a single tensor, classically modelled by (Bear 1988; Delgado 2007)

$$ D = (D_m + \alpha_T |v|) I + (\alpha_L - \alpha_T) (v \otimes v)/|v|, \qquad (2.2) $$

where $\alpha_L$ and $\alpha_T$ are the longitudinal and transverse dispersivities, respectively, $D_m$ is the molecular diffusivity, and $I$ is the identity matrix. Two dimensionless numbers, the Péclet number $Pe$ and the Damköhler number $Da$, characterise respectively the ratio of diffusion and advection times and the ratio of diffusion and reaction times. Setting $L$ as a characteristic length scale, $U$ as a characteristic velocity scale, $C_A^0$ as the concentration of reactant A far from the mixing front,
<figure><img src="image_3.png" /><figcaption>Figure 1. Flow and boundary conditions corresponding to the co-flow mixing front. Reactants A and B are segregated and co-injected at the left-hand boundary. They flow, mix and react in the domain $x > 0$. The colour scale indicates local reaction rates obtained in numerical simulations of the hydrodynamic dispersion approximation (see § A.4).</figcaption></figure>

$$Pe = \frac{UL}{D_m} \qquad (2.3)$$

and

$$Da = \frac{kC_A^0 L^2}{D_m}. \qquad (2.4)$$

We consider a mixing front in which reactant concentrations vary in the $(x, y)$ plane, with $x$ the coordinate along the front, and $y$ the coordinate in the transverse direction, and are invariant in the $z$ direction. The two reactants are initially fully segregated (i.e. $(C_A, C_B) = (C_A^0, 0)$ for $y > 0$, and $(C_A, C_B) = (0, qC_A^0)$ for $y < 0$). We define $q = C_B^0/C_A^0$ as the ratio of reactant concentrations on both sides of the front. The reaction rate is $R = kC_A C_B$. We characterise the properties of the reaction front as a function of the position $x$ along the flow direction using three metrics: the maximum local reaction rate $R_{max}$ defined as

$$R_{max}(x) = \max_{y} R(x, y), \qquad (2.5)$$

the total reaction intensity $I$ defined as the integrated reaction rate along $y$,

$$I(x) = \int dy R(x, y), \qquad (2.6)$$

and the width of the reactive front $s_R$ defined from the second spatial moment of the reaction rate as

$$s_R^2 = \frac{\int dx y^2 R(x, y)}{\int dx R(x, y)}. \qquad (2.7)$$

From these definitions, $s_R \sim I/R_{max}$ is expected.

## 2.1. Diffusion-reaction front (co-flow)

We first consider a mixing front with two reactive fluids flowing in parallel with constant velocity $v = (U, 0)$ (figure 1) and mixing transversely through diffusion only (low Pe), such that (2.2) reduces to $D \equiv D_m$. As the flow is steady and uniform, the steady advection-diffusion problem simplifies to an unsteady diffusion problem in the Lagrangian coordinate system, where the Lagrangian time is defined by $t = x/U$. From (2.1), the concentration of solute species is governed by the equations
$$ 
\begin{aligned}
\frac{\partial C_A}{\partial t} &= D_m \frac{\partial^2 C_A}{\partial y^2} - k C_A C_B, \\
\frac{\partial C_B}{\partial t} &= D_m \frac{\partial^2 C_B}{\partial y^2} - k C_A C_B.
\end{aligned}
\qquad (2.8) $$

The properties of such a reactive front were studied by Gálfi & Rácz (1988), Taitelbaum et al. (1991) and Larralde et al. (1992), and we recall here their main results for completeness. For an instantaneous reaction (infinite Da), reactants react as soon as they meet. Thus they remain fully segregated, and the reaction intensity is driven solely by the mixing flux across the interface. For finite Da, two effective reaction regimes must be distinguished. At times t smaller than the characteristic reaction time $t_r = 1/(C_A^0 k)$, i.e. a distance smaller than $t_r U$, the reaction is slow compared to diffusion, and the front is in a reaction-limited regime. The concentration profiles of reactants can thus be approximated with the conservative concentration profiles. The characteristics of the reactive front follow typical diffusive scaling. When replacing time by position along front $t = x/U$, this leads to,

$$ s_r \sim L Pe^{-1/2} \left(\frac{x}{L}\right)^{1/2}, \quad R_{max} \sim (C_A^0)^2 k \quad \text{and} \quad I \sim (C_A^0)^2 L k Pe^{-1/2} \left(\frac{x}{L}\right)^{1/2}. \quad (2.9) $$

In contrast, at times larger than the reaction time ($t > t_r$, $x > t_r U$), the front becomes mixing-limited, resulting in different scaling of the reactive front properties with distance, Pe and $D_m$ (Taitelbaum et al. 1991). As shown in Larralde et al. (1992), an approximate solution can be obtained by a perturbation method that uses a spatial linearisation of the conservative profile $F = C_A - C_B$ near the front origin, while neglecting nonlinear terms (see Appendix A). This leads to

$$ 
\begin{aligned}
s_r &\sim L Da^{-1/3} Pe^{-1/6} \left(\frac{x}{L}\right)^{1/6}, \\
R_{max} &\sim \frac{C_A^0 D_m}{L^2} Da^{1/3} Pe^{2/3} \left(\frac{x}{L}\right)^{-2/3}, \\
I &\sim \frac{C_A^0 D_m}{L} Pe^{1/2} \left(\frac{x}{L}\right)^{-1/2}.
\end{aligned}
\qquad (2.10) $$

Note that although the reaction is fast compared to mixing in this regime, the reaction kinetics, encoded in Da, still appears in the scaling laws for $s_r$ and $R_{max}$. As Da increases, the width of the reactive zone becomes narrower because the time in which the two species can co-exist at the same location decreases. As $s_r$ decreases, the maximum reaction rate increases proportionally such that the product $s_r \times R_{max}$ is independent of Da. This is a consequence of the independence of the reaction intensity $I \sim s_r \times R_{max}$ of Da, as it is proportional to the mixing rate in the mixing-limited regime. For infinite Da, $s_r$ goes to zero, $R_{max}$ goes to infinity, and $I$ remains proportional to the mixing rate. For simplicity, we did not report here the scaling with $q$ for the case of an asymmetric reaction in these expressions (see Appendix A for the full scaling laws).

## 2.2. Advection–diffusion–reaction with compression (saddle flow)

We now consider the case of a reactive front in converging flows, considering a saddle flow around a stagnation point located at $(x, y) = (0, 0)$:

$$ v = \gamma \begin{pmatrix} x \\ -y \end{pmatrix}, \qquad (2.11) $$
<figure><img src="image_4.png" /><figcaption>Figure 2. Flow and boundary conditions corresponding to the saddle mixing front. Reactants A and B are segregated and co-injected at the top and bottom boundaries, respectively. They flow, mix and react in the domain. The stagnation point is located at $(x, y) = (0, 0)$. The colour scale indicates local reaction rates obtained in numerical simulations of the hydrodynamic dispersion approximation (see § A.4).</figcaption></figure>

where $\gamma$ is the compression rate (Haller 2015). In this flow field, opposing flows converge and diverge at the stagnation point along directions $y$ and $x$, respectively, leading to the formation of one stable manifold ($y$) and one unstable manifold ($x$) that are perpendicular to each other, and an exponential compression of fluid parcels. For this system, we define the characteristic velocity as $U = \gamma L$ and the Péclet number as $Pe = \gamma L^2/D_m$.

The saddle flow leads to the formation of a reactive front along the $x$ direction, since the two reactants are supplied in opposite directions on the $y$ axis (figure 2). In a Lagrangian coordinate system attached to a fluid particle travelling on the $x$ axis, the concentration of the reactant approximately follows:

$$ \frac{\partial C_A}{\partial t} - \gamma y \frac{\partial C_A}{\partial y} = D_m \frac{\partial^2 C_A}{\partial y^2} - k C_A C_B. \quad (2.12) $$

Equation (2.12) implicitly assumes that longitudinal scalar gradients are much smaller than transverse ones, which is generally the case in steady mixing front. A detailed derivation is provided in Rousseau et al. (2023).

As first shown by Ranz (1979) in the conservative case, this equation can be solved by a change of variables that integrates the deformation of fluid elements in a Lagrangian reference system (see Villermaux (2019) for a review). Under this transformation, the advection-diffusion equation simplifies to a diffusion equation with analytical solutions. Bandopadhyay et al. (2017) and Izumoto et al. (2023) have shown that reactive transport can also be studied in this Lagrangian framework by using the reaction-diffusion theory of Larralde et al. (1992). Two regimes can be distinguished in the saddle-flow configuration: first, a kinetic-limited regime when $\gamma \gg t_r^{-1}$, that is, when reaction at the front is limited by the time for reaction to occur rather than by the supply of reactants; second, a mixing-limited regime for $\gamma \ll t_r^{-1}$, where reaction is limited by the supply of reactant under the action of compression. In the kinetic-limited regime, Izumoto et al. (2023) have shown that the properties of the reaction front under diffusion and compression are

$$ s_R \sim L Pe^{-1/2}, \quad R_{max} \sim (C_A^0)^2 k \quad \text{and} \quad I \sim (C_A^0)^2 kL Pe^{-1/2}, \quad (2.13) $$

while in the mixing-limited regime, the front scales as

$$ \begin{aligned} s_R &\sim L D a^{-1/3} P e^{-1/6}, \\ R_{max} &\sim \frac{C_A^0 D_m}{L^2} D a^{1/3} P e^{2/3}, \\ I &\sim \frac{C_A^0 D_m}{L} P e^{1/2}. \end{aligned} \quad (2.14) $$

More details on the derivation are given in Appendix A and in Izumoto et al. (2023).
## 2.3. Advection-dispersion-reaction in porous media: co-flow

Here, we investigate how the presence of a homogeneous porous media of typical length scale $L=d$, the grain diameter, and the resulting hydrodynamic dispersion (Delgado 2007) affect the properties of the reaction front for both parallel and converging flows. We assume that $\alpha_T U$ is large compared to $D_m$ (high Pe), so molecular diffusivity can be neglected in (2.2). In the co-flow scenario, the velocity is constant, so the dispersion tensor is also constant. The dispersive mixing front is therefore very similar to the diffusive one. Hence the scaling laws corresponding to diffusive transport ((2.9) and (2.10)) still hold for dispersive transport, with the molecular diffusion coefficient $D_m$ replaced by the transverse dispersion coefficient

$$D_T = \alpha_T U = D_m Pe \frac{\alpha_T}{d}, \qquad (2.15)$$

the diffusive Péclet number Pe replaced by a dispersive Péclet number

$$Pe_D = \frac{LU}{D_T} = \frac{d}{\alpha_T}, \qquad (2.16)$$

and the diffusive Damköhler number Da replaced by a dispersive Damköhler number

$$Da_D = \frac{Da}{Pe} \frac{d}{\alpha_T}. \qquad (2.17)$$

We thus obtain in the kinetic-limited regime $x < t_r U$,

$$s_r \sim (\alpha_T x)^{1/2}, \quad R_{max} \sim (C_A^0)^2 k \quad \text{and} \quad I \sim (C_A^0)^2 k (\alpha_T x)^{1/2}, \qquad (2.18)$$

and in the mixing-limited regime $x > t_r U$,

$$s_r \sim d Da^{-1/3} Pe^{1/3} \left(\frac{d}{\alpha_T}\right)^{-1/2} \left(\frac{x}{d}\right)^{1/6}, \\ R_{max} \sim \frac{C_A^0 D_m}{d^2} Da^{1/3} Pe^{2/3} \left(\frac{x}{d}\right)^{-2/3}, \\ I \sim \frac{C_A^0 D_m}{d} Pe \left(\frac{\alpha_T}{x}\right)^{1/2}. \qquad (2.19)$$

## 2.4. Advection-dispersion-reaction in porous media: saddle flow

For the case of saddle flow, Rousseau et al. (2023) have shown that the dispersive front has a shape similar to that of the dispersive co-flow configuration, with a prefactor in the mixing width. The velocity varies spatially along the front as $v(x, 0) = (\gamma x, 0)$. Thus the transverse dispersivity is spatially dependent:

$$D_T = \alpha_T \gamma x = D_m Pe x \alpha_T / d^2. \qquad (2.20)$$

The dispersive Péclet and Damköhler numbers thus read $Pe_D = d^2/(\alpha_T x)$ and $Da_D = Da Pe^{-1} d^2/(\alpha_T x)$, respectively. In the kinetic-limited regime, the diffusive scaling (2.13) thus transforms to the dispersive scaling

$$s_r \sim (\alpha_T x)^{1/2}, \quad R_{max} \sim (C_A^0)^2 k \quad \text{and} \quad I \sim (C_A^0)^2 k (\alpha_T x)^{1/2}. \qquad (2.21)$$
<table><caption>Table 1. Summary of expected scaling laws as a function of Péclet number Pe and distance x at Darcy scale using the hydrodynamic dispersion equation, where x is the distance measured along the downstream flow direction.</caption><thead><tr><th>Flow config.</th><th>Regime</th><th>Equation</th><th>Reactive width s<sub>R</sub></th><th>Maximum reaction rate R<sub>max</sub></th><th>Total intensity I</th></tr></thead><tbody><tr><td>Co-flow</td><td>Kin. lim.</td><td>(2.18)</td><td>x<sup>1/2</sup></td><td>Const.</td><td>x<sup>1/2</sup></td></tr><tr><td>Co-flow</td><td>Mix. lim.</td><td>(2.19)</td><td>Pe<sup>1/3</sup> x<sup>1/6</sup></td><td>Pe<sup>2/3</sup> x<sup>-2/3</sup></td><td>Pe x<sup>-1/2</sup></td></tr><tr><td>Saddle flow</td><td>Kin. lim.</td><td>(2.21)</td><td>x<sup>1/2</sup></td><td>Const.</td><td>x<sup>1/2</sup></td></tr><tr><td>Saddle flow</td><td>Mix. lim.</td><td>(2.22)</td><td>Pe<sup>1/3</sup> x<sup>1/2</sup></td><td>Pe<sup>2/3</sup></td><td>Pe x<sup>1/2</sup></td></tr></tbody></table>

In the mixing-limited regime, the diffusive scaling (2.14) transforms into the dispersive scaling

$$ 
\begin{align} 
s_r &\sim D a^{-1/3} P e^{1/3} (\alpha_T x)^{1/2}, \notag \\ 
R_{max} &\sim \frac{C_A^0 D_m}{d^2} D a^{1/3} P e^{2/3}, \notag \\ 
I &\sim \frac{C_A^0 D_m}{d^2} P e (\alpha_T x)^{1/2}, \tag{2.22} 
\end{align} 
$$

where Pe = γd²/D_m, and Da is defined in (2.4).

All theoretical scaling laws for reactive transport in the presence of hydrodynamic dispersion are summarised in table 1. In § A.4, we numerically validate these scaling laws for conditions close to the experiments. In the next section, we compare these analytical scaling laws with experimental data. Note that these results are derived by assuming that mixing can be modelled by hydrodynamic dispersion, hence ignoring incomplete mixing at the pore scale. Therefore, a discrepancy of these predictions with observations can be used to assess the impact of the pore-scale mixing dynamics on macroscopic reaction rates in experiments.

## 3. Experimental results

We investigated steady reactive fronts that form between two reactive fluids flowing in parallel (co-flow) or towards each other (saddle flow) in 3-D porous media. Our imaging set-up measures two-dimensional maps of the effective reaction rate integrated over several pore diameters across the cell gap. To this end, we used the principle of chemiluminescence and refractive index matching. The chemiluminescence reaction occurs through the mixing of reactants at the pore scale, producing photons that travel freely through a transparent, index-matched porous material. A camera records the total light emitted over the thickness of the cell, which corresponds to a macroscale effective reaction rate mapped across the mixing front. Experimental methods are described in the following.

### 3.1. Methods

#### 3.1.1. Chemiluminescence reaction

We use the chemiluminescence of luminol (Izumoto et al. 2023), which involves a catalytic reaction of H₂O₂ with Co²⁺ followed by an oxidation reaction of the luminol compound with OH· and O₂· radicals (Uchida et al. 2004). These reactions can be written as

$$ \text{luminol} + \text{H}_2\text{O}_2 + 2\text{OH}^- \xrightarrow{\text{Co}^{2+}} \text{3-aminophthalatedianion} + \text{N}_2 + \text{H}_2\text{O} + \text{photon} \quad (3.1) $$
<figure><img src="image_5.png" /><figcaption>Figure 3. (a) Co-flow set-up and (b) saddle-flow set-up. The distance between the inlet/outlet and the stagnation point is 103 mm. For the porous medium, we used a larger cell; the distance between the inlet/outlet and the stagnation point is 208 mm. In saddle flow, a is 303 mm² for the Hele-Shaw cell, and 811 mm² for the cell for porous media. The thick arrows indicate the flow direction.</figcaption></figure>

luminol is oxidised to 3-aminophtalatedianion with an emission of a blue photon of wavelength $\lambda \approx 450$ nm. Such a reaction can be approximated (Matsumoto & Matsuo 2015) as a bimolecular reaction:

$$A + B \rightarrow C + \text{photon}, \qquad (3.2)$$

where A and B represent the H₂O₂ and the luminol species, respectively. The total reaction rate is thus directly proportional to the intensity of the emitted light. In practice, we prepared two solutions. One was a mixture of 1 mM luminol, 7 mM NaOH and 0.01 mM CoCl₂ (termed luminol solution), and the other was a mixture of 0.5 mM H₂O₂ and 3.9 mM NaCl (termed H₂O₂ solution). Izumoto et al. (2023) have estimated the reaction constant $k = 0.08 \text{ s}^{-1} \text{ mM}^{-1}$ by mixing the two solutions in a beaker and measuring the intensity of light over time. This leads to the characteristic time scale $t_r = 1/(kC_0) = 2.5$ s, where $C_0 = 0.5$ mM is the concentration of luminol in the batch solution (i.e. half of the concentration of the solution before mixing). Note that after the fast reaction phase, the luminol reaction continues with a slower reaction. The light emitted from this slow reaction was subtracted from the original images by defining a threshold in the light intensity following the method of Izumoto et al. (2023). Note also that the commercial H₂O₂ solution may contain stabiliser compounds leading to a departure from the expected bimolecular kinetic. This did not occur for our 30% H₂O₂ solution, which contained the Sigma-Aldrich stabiliser.

### 3.1.2. Flow cells

To impose the co-flow and saddle-flow conditions, we designed two flow cells. The first cell includes two inlets in two separated triangular-shaped branches (figure 3a), leading to a co-flow configuration. By injecting two different solutions from each of the inlets, they start mixing at the edge of the separator, and they flow in parallel towards the outlet at the other side of the cell, placed 220 mm from the start of the mixing zone. The second cell has four flow branches (figure 3b), with two inlets and two outlets. The cell boundary is chosen to be a streamline of the saddle flow (2.11), i.e. $y = \pm a/x$, with a constant. The dimensions of the cells are given in figure 3. The flow cells are kept empty or filled by a transparent granular material matched to the water refractive index (fluorinated ethylene propylene, FEP) of size 2 mm. Empty cells, referred to as Hele-Shaw cells, are used to study the advection-diffusion-reaction front. They are thin (2 mm) to maintain viscous
<table><caption>Table 2. Summary of experimental runs and conditions. Constant parameters are the Damköhler number Da = 1600, the molecular diffusion coefficient Dm = 10⁻⁹ m² m⁻² s⁻¹, and the characteristic length L = d = 2 mm. For co-flow, Pe = UL/Dm. For saddle flow, U = γL, leading to Pe = γL²/Dm. For each flow configuration (columns) the different runs are indicated inside brackets.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Co-flow</th><th colspan="2">Saddle flow</th></tr><tr><th>Hele-Shaw cell</th><th>Porous media</th><th>Hele-Shaw cell</th><th>Porous media</th></tr></thead><tbody><tr><td>Injection rate (ml h<sup>-1</sup>)</td><td>{183}</td><td>{864, 594, 409, 281, 193, 133, 91.4, 62.8, 43.2}</td><td>{432, 297, 204, 140, 96.6}</td><td>{1913, 1316, 904.6, 622.1, 427.8, 294.2, 202.3, 139.1, 95.66}</td></tr><tr><td>Velocity (Darcy scale) (mm s<sup>-1</sup>)</td><td>{0.845}</td><td>{1.79, 1.23, 0.845, 0.581, 0.400, 0.275, 0.189, 0.130, 0.0894}</td><td></td><td></td></tr><tr><td>Compression rate × 10<sup>3</sup> (s<sup>-1</sup>)</td><td>0</td><td>0</td><td>{43.4, 29.8, 20.5, 14.1, 9.70}</td><td>{43.4, 29.8, 20.5, 14.1, 9.70, 6.67, 4.59, 3.16, 2.17}</td></tr><tr><td>Pe</td><td>{1690}</td><td>{3575, 2458, 1690, 1162, 799.3, 549.7, 377.0, 259.9, 178.7}</td><td>{8678, 5968, 4104, 2822, 1941}</td><td>{8678, 5968, 4104, 2822, 1941, 1334, 917.6, 631.0, 433.9}</td></tr></tbody></table>

flows (low Reynolds number). Porous media cells are thicker (12 mm) to allow for the packing of a few grain diameters. Since the refractive index of FEP (1.34) is close to that of water (1.33), this allowed us to visualise an integrated image of the reaction rate over the whole cell depth. The porosity of the packed FEP was estimated to be 0.37 by comparing the saturated and unsaturated weights of the cell. Since the scalar gradients developing in co-flows and saddle flows mainly involve mixing in the direction transverse to flow (y), we assume that Taylor-Aris dispersion due to the Poiseuille velocity profile, which operates mainly in the direction of flow, only weakly alters mixing rates. For saddle-flow experiments, all applied compression rates were such that the characteristic compression time was larger than the reaction time, γ⁻¹ > t_r. Hence they are in the mixing-limited regime. A summary of the experimental runs and associated parameters is provided in table 2.

### 3.1.3. Protocol

Before each experiment, we filled the cell with deionised water. For experiments with porous media, we injected CO₂ gas before filling the cell with water. Hence any trapped CO₂ gas dissolved in the water, and there were no remaining bubbles. Then we injected luminol and H₂O₂ solutions from two different inlets at a fixed injection rate. We imposed nine flow rates for each of the experimental configurations. The maximum Reynolds number varied between 0.18 and 3.6 for co-flow cells, and between 0.22 and 4.5 for saddle-flow cells (Re = Ud/v, with v = 10⁻⁶ m² s⁻¹ the kinematic viscosity, d = 2 mm the grain size, and U = Q/Aφ, where Q is the volumetric injection rate, A is the cross-sectional area of the tank at injection, and φ is the porosity). We thus assume that inertial effects are small in our experiments. For each flow rate, we waited for the front to stabilise before acquiring images with a full-frame high-sensitivity 14-bit camera (SONY alpha7s) with a macro lens (MACRO GOSS F2.8/90). The image resolution was 0.046 mm per pixel for all cases. For porous media experiments, we triplicate the experiments by repacking the FEP to obtain results independently of the specific grain packing. The images were rescaled by the bit depth of the camera to obtain the normalised reaction rate.

The luminol reaction can be approximated as a bimolecular reaction after subtracting light emission from the slow reaction kinetics at later times (Izumoto et al. 2023).
<figure><img src="image_6.png" /><figcaption>Figure 4. Experimental reaction rate fields in the Hele-Shaw cell for (a) the co-flow configuration ($Pe = 3575$) and (b) the saddle-flow configuration ($Pe = 8678$). Flow is from left to right. The white bar represents 10 mm. The colour scale indicates local reaction rates $R$ normalised by the bit depth of the camera.</figcaption></figure>

We measured this weak light emission as follows. We first mixed the luminol solution and the $H_2O_2$ solution in a beaker, and then injected this solution into Hele-Shaw and porous media cells. After more than 30 minutes, we started taking images. The intensity of these pictures corresponds to the light emission from the long-lasting reaction. We subtracted these image intensities from the images taken in reactive transport experiments following the protocol of Izumoto et al. (2023).

We estimated the transverse width of the reactive front $s_R$ by computing the spatial standard deviation of a normalised transverse reaction rate profile. The maximum local reaction rate $R_{max}$ was taken as the maximum, and the total reaction intensity $I$ from the image intensity profile. The width was obtained experimentally by

$$s_R(x) = \left( \int y^2 P(x, y) dy - \left( \int y P(x, y) dy \right)^2 \right)^{1/2}, \quad (3.3)$$

where $y$ is the transverse position, and $P(x, y)$ is the normalised image intensity at the position $(x, y)$. The maximum image intensity of the profile was chosen as $R_{max}(x) = \max_y P(x, y)$, and the image intensity was integrated over the measured line to calculate $I_R(x) = \int P(x, y) dy$.

## 3.2. Reactive fronts in the Hele-Shaw cell

We first used Hele-Shaw cells without porous medium to confirm that the chemiluminescence reaction is a good proxy for bimolecular reactive fronts, whose theoretical behaviour is now well understood in the absence of porous media (Gálfi & Rácz 1988). In the co-flow configuration, we observe a steady reactive front that forms at the interface between the flowing reactants reactants (figure 4a). Luminescence, which is proportional to the local reaction rate, is maximal at the entrance of the cell, and rapidly decreases as it moves downstream. The dependence of the measured reaction zone width, maximum reaction rate and total intensity with advective time (figure 5a) is well captured by the theoretical scaling laws expected for a bimolecular reaction in a mixing-limited reactive front (2.10). We could not observe the kinetic-limited regime, which occurs at short times ($t < t_r$), probably because the cell geometry did not allow for a perfectly 2-D co-injection of reactants.

In the saddle-flow configuration (figure 4b), the reaction front is steady and invariant along the $x$ axis. In fact, the constant compression rate induced by the saddle-flow configuration maintains a fixed scale of scalar fluctuations (Rousseau et al. 2023). We plot the characteristics of the reaction zone against the Péclet number in a window of 100 pixels centred on the stagnation point (figure 5b). The data agree well with the scaling laws expected for the mixing-limited regime (2.14), which is expected since $\gamma^{-1} > t_r$ for all compression rates used experimentally. The consistency of these observations with
<figure><img src="image_7.png" /></figure>

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

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

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

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

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

Figure 5. Experimental reaction front properties (red circles) in the Hele-Shaw cells for (a) co-flow at $Pe = 1690$ and (b) saddle-flow configurations. Errors bars are calculated by taking the standard deviation over a window of length $L$. The black dashed lines show the theoretical scaling laws given by (2.10) (co-flow) and (2.14) (saddle flow). The red continuous lines are power-law fits.

the theoretical scaling laws (Gálfi & Rácz 1988; Larralde et al. 1992) confirms that the chemiluminescence reaction is well approximated by a bimolecular reaction.

### 3.3. Steady reactive fronts in porous media

Experimental images of reactive fronts obtained in porous media cells in co-flow and saddle-flow configurations are shown in figure 6 for increasing Péclet numbers. The presence of grains clearly enhances the transverse spreading of the reactive zone compared to the open flow case (figure 4). In contrast to the Hele-Shaw cell, the reactive mixing interface that forms in the saddle flow in porous media broadens with $x$ (figure 6b), since the dispersion coefficient increases with $x$. This broadening was also observed in conservative transport experiments (Rousseau et al. 2023). We present in figure 7 the characteristics of the reactive front and its dependence on the distance $x$ for the co-flow and for the saddle flow. In co-flow, we normalise the distance by the characteristic reaction distance $x_r = U/t_r$ to highlight the two different regimes in space. For saddle flow, we normalise $x$ by grain size since there is only one regime in space. We averaged these characteristics over triplicate experiments with different grain packings.

For co-flow experiments, we observe first a kinetic-limited regime, where the reaction intensity $I$ increases with distance and then, further downstream, a mixing-limited regime, where $I$ decays (figure 7a1). The transition distance between the two regimes is obtained at $x \approx x_r$. In the kinetic-limited regime, $R_{max}$ and $I$ increase faster than predicted by the hydrodynamic dispersion theory (2.19). In the mixing-limited regime ($x > x_r$), the decay of $R_{max}$ and $I$ is closer to the hydrodynamic dispersion model (2.10). For saddle-flow experiments, the evolution of $R_{max}$ and $I$ with distance is also significantly faster than predicted by the hydrodynamic dispersion theory for the mixing-limited regime (2.22).
<figure><img src="image_13.png" /><figcaption>Figure 6. Normalised steady-state reaction rate fields for (a) co-flow and (b) saddle flow, obtained from chemiluminescence experiments in porous media. The reaction rate is normalised by the maximum reaction rate at the highest Pe for each configuration. The white scale bar has length 10 mm. The porosity appears as dark patches in the mixing front. For (a), the left-hand edge corresponds to the start of mixing and for (b), the left-hand edge corresponds to the stagnation point.</figcaption></figure>

Interestingly, the measured reactive mixing width $s_R$ in both co-flow and saddle flows in the presence of porous media approximately follows the hydrodynamic dispersion predictions (2.18), (2.19), (2.22). This suggests that the hydrodynamic dispersion theory provides a robust framework to predict the spatial extent of the mixing fronts, as observed by Rousseau *et al.* (2023) for conservative mixing, although it does not capture actual mixing and reaction rates.

We plot the properties of the experimental reaction front as a function $Pe$ measured at $x = 10d$ in figure 8. At this position, all experiments are in the mixing-limited regime. For the intensity of the reaction front intensity (figure 8a1,b1) and the maximum reaction rate (figure 8a2,b2), the evolution with the Péclet number is stronger than predicted by the hydrodynamic dispersion theory. This discrepancy suggests that incomplete pore-scale mixing plays an important role, as discussed in the following. Note that the evolution of the reactive width with the Péclet number (figure 8a3,b3) is consistent with the scaling predicted by hydrodynamic dispersion, since this spatially integrated measure
<figure><img src="image_14.png" /></figure>

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

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

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

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

<figure><img src="image_19.png" /><figcaption>Figure 7. Reactive mixing properties measured in porous media. (a) The reaction front properties over distance normalised by the characteristic reaction distance $x_r$ in co-flow porous media experiments, showing (a1) reaction intensity, (a2) maximum reaction rate, and (a3) reaction width. The dashed lines show the scaling laws expected from the hydrodynamic dispersion theory (2.18), (2.19). (b) Reaction front properties over distance normalised by the grain diameter in saddle-flow porous media experiments, showing (b1) reaction intensity, (b2) maximum reaction rate, and (b3) reaction width. The dashed lines show the scaling expected from the hydrodynamic dispersion theory (2.22). Continuous lines show the fit to the data. The error margin is estimated from the minimum and maximum envelope of the three experimental replicates obtained with different packing realisations.</figcaption></figure>

is not modified by the presence of pore-scale gradients, as also found in conservative experiments (Rousseau *et al.* 2023).

## 4. Discussion

Experimental observations show that reaction rates in mixing fronts under flow in porous media largely depart from those predicted by the hydrodynamic dispersion framework. In this section, we discuss the origin of this discrepancy and investigate the role played by pore-scale chaotic advection (Lester, Metcalfe & Trefry 2013; Heyman *et al.* 2020), which
<figure><img src="image_20.png" /></figure>

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

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

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

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

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

Figure 8. The reaction front properties as a function of Pe for porous media experiments measured at fixed distance $x = 20d \pm 10d$. (a) Co-flow porous media with (a1) reaction intensity, (a2) maximum reaction rate, and (a3) reaction width. (b) Saddle-flow porous media with (b1) reaction intensity, (b2) maximum reaction rate, and (b3) reaction width. Error bars are estimated from the variability of reaction front properties at various distances between $x = 10d$ and $x = 30d$. The 95% confidence interval is shown. Note that the width at the lower two Pe values, and the maximum reaction rate and intensity at the lowest Pe, could not be obtained because of the small image intensity. The dashed lines show the scaling expected from the hydrodynamic dispersion theory ((2.18), (2.19) and (2.22)). The scaling of the effective reaction intensity $I$ predicted by the micro-scale mixing theory ((4.16) and (4.20)) is shown as thick continuous lines.

is expected to sustain chemical gradients and drives the rate of mixing at the pore scale. We focus on the total intensity of the reaction *I* along the reactive front, because it is a robust integrated measure of effective reaction rates.

## 4.1. *Stretching-enhanced reactive front kinetics under co-flow*

At the pore scale, we define the mixing interface as the isoconcentration surface where $C_A = C_B$. This surface coincides with the conservative isoconcentration surface $F=0$, with the conservative component $F = C_A - C_B$. In a given cross-section perpendicular to the flow (figure 9), the length of this interface results from a balance between stretching, which continuously deforms the interface (figure 9) and aggregation/dilution, which reduces its length by merging and diluting the stretched sections of the front (Le Borgne, Dentz & Villermaux 2013; Villermaux 2018, 2019).
<figure><img src="image_26.png" /><figcaption>Figure 9. Sketch of the micro-scale mixing interface between reactants A and B. (a) Starting from a given mixing interface (black line), fluid stretching leads to an incremental elongation of the interface (red line). The latter is balanced by two main processes that contribute to the reduction of the interface. The first is dilution, which describes the decay of concentration of an isolated lamella of A into the B side (or the opposite). The second is aggregation, which describes the overlap of difference sections of the front. (b) Considering a lamella of A being stretched into the B side, the section of width equal to the Batchelor scale dilutes exponentially, leading to a receding of the interface and thus a flux of A into the B side.</figcaption></figure>

Due to stretching, the mixing front forms an ensemble of lamellae that are elongated in one direction and compressed in the other (figure 9a). The elongation of a lamella is defined as $\rho(t) = \ell(t)/\ell_0$, where $\ell_0$ and $\ell(t)$ are the initial and final lengths of the lamella, respectively. The chaotic nature of pore-scale stretching (Heyman et al. 2020; Souzy et al. 2020) implies that this elongation should be exponential, e.g. $\rho = e^{\lambda t}$, with

$$ \lambda \sim \frac{U}{d} \sim \frac{D_m}{d^2} Pe, \qquad (4.1) $$

the Lyapunov exponent. At pore scale, the thickness of lamella $s$ is reduced by compression and enhanced by diffusion, according to

$$ \frac{ds}{dt} = -\lambda s + \frac{D_m}{s}. \qquad (4.2) $$

When compression and diffusion balance each other, i.e. when the temporal derivative on the right-hand side of (4.2) is zero, the lamellar width converges to the Batchelor scale (Batchelor, Howells & Townsend 1959),

$$ s \to s_B = \sqrt{\frac{D_m}{\lambda}} \sim d Pe^{-1/2}. \qquad (4.3) $$

This equilibrium is achieved at the mixing time (Villermaux 2019)

$$ t_m \sim \frac{\log Pe}{2\lambda}, \qquad (4.4) $$

which is obtained when the lamella width reaches the Batchelor scale by compression from the initial width, $s_0 e^{-\lambda t_m} = s_B$.

After mixing time, further elongation of lamellae of constant width $s_B$ in a finite-size domain implies that they overlap and aggregate, limiting the growth of the total interface length (figure 9a). Furthermore, lamellae of A that are stretched in the B domain (or B into A) start to dilute exponentially, leading to receding of the interface after an incremental deformation (figure 9b). In this regime, the balance between constant stretching and destruction of the interface thus leads to a saturation of the front elongation
to a maximum value $\rho_m$ equal to the elongation at the mixing time (Villermaux 2018):

$$ \rho_m = \rho(t_m) = e^{\lambda t_m} \sim \sqrt{Pe}, \qquad (4.5) $$

where we have used (4.4) and (4.1).

We first consider an element of the interface that is stretched in the B domain at constant stretching rate (figure 9b). After the mixing time, the mixing interface reaches a steady state, and the tip of the lamella, of width equal to the Batchelor scale, dilutes exponentially. This dilution leads to a flux of A to the B side, which is proportional to the stretching rate:

$$ j_A \sim C_A \lambda \sim C_A \frac{D_m}{d^2} Pe, \qquad (4.6) $$

where $C_A$ is the characteristic concentration of A. When the mixing interface reaches a steady length $\rho_m$, it is composed of a characteristic number of elongated elements $n \propto \rho_m/L$, with L the domain size. Using (4.5), we estimate

$$ n = \beta \sqrt{Pe}, \qquad (4.7) $$

with $\beta$ a constant prefactor to be determined. The total flux across the interface due to the dilution of these stretched segments of the interface is thus

$$ J_A = nj_A = \beta \frac{D_m}{d^2} C_A Pe^{3/2}. \qquad (4.8) $$

Initially, reactants are fully segregated, and $C_A$ and $C_B$ remain close to their initial concentrations $C_A^0$ and $qC_A^0$. At a further distance downstream, called $x_\Delta$, dispersion inevitably smooths scalar gradients at the interface such that $C_A < C_A^0$, limiting the flux of A in (4.6), and equivalently for the flux of B. We focus here on the first regime, where $x < x_\Delta$ and $C_A \approx C_A^0$ ($C_B \approx qC_A^0$), which covers a large part of our measurements. Note that in the case of fast reactions (large Damköhler numbers), reactions occur mostly at the mixing interface where concentrations are similar. Here, we consider the general case of arbitrary Damköhler numbers where reactions can be slow compared to mixing rates. In this case, reactions also occur further away from the mixing interface. On the B side, the concentration of A becomes rapidly much smaller than that of B as it is diluted by stretching-enhanced diffusion. The evolution of the mass of A on the B side, denoted $m_{A|B}$, can thus be approximated as a first-order reaction process:

$$ \frac{dm_{A|B}}{dt} = -kC_B m_{A|B} + J_A \approx -\frac{q}{t_r} m_{A|B} + \beta \frac{C_A^0 D_m}{d^2} Pe^{3/2}, \qquad (4.9) $$

where we used $C_B \approx qC_A^0$ for $x < x_\Delta$ and for the general case of an asymmetric reaction front. This equation leads to

$$ m_{A|B} = \beta \frac{C_A^0 D_m t_r}{qd^2} Pe^{3/2} (1 - e^{-qt/t_r}). \qquad (4.10) $$

Similar arguments can be used to estimate $m_{B|A}$, leading to

$$ m_{B|A} = \beta \frac{C_A^0 D_m t_r q}{d^2} Pe^{3/2} (1 - e^{-t/t_r}). \qquad (4.11) $$

Thus the total reaction intensity developing on the two sides of the mixing front is

$$ I = kC_B^0 m_{A|B} + kC_A^0 m_{B|A} = \beta q \frac{C_A^0 D_m}{d^2} Pe^{3/2} (2 - e^{-t/t_r} - e^{-qt/t_r}). \qquad (4.12) $$
S. Izumoto, G. Rousseau, T. Le Borgne and J. Heyman

Replacing $t$ by $x/U$ and using the definitions of $Pe$ and $Da$ (see (2.3) and (2.4)), the last equation is

$$I = \beta q \frac{C_A^0 D_m}{d^2} Pe^{3/2} \left( 2 - \exp \left( -\frac{Da}{Pe} \frac{x}{d} \right) - \exp \left( -q \frac{Da}{Pe} \frac{x}{d} \right) \right). \quad (4.13)$$

This equation predicts a linear growth of reaction intensity followed by a plateau (figure 10a). In the kinetic-limited regime $t \ll t_r$, the consumption rate is much smaller than the mixing flux. Hence the mass of A in the B domain grows linearly as

$$m_{A|B} = J_A t, \quad (4.14)$$

which corresponds to a Taylor expansion of (4.11) at small times $t \ll t_r$. The intensity of the reaction is thus

$$I = \beta q \frac{(1+q) D_m C_A^0}{d^2} Da Pe^{3/2} \frac{t}{t_r} \sim \beta q \frac{(1+q) D_m C_A^0}{d^2} Da Pe^{1/2} \frac{x}{d}. \quad (4.15)$$

In the mixing-limited regime $t \gg t_r$, the reaction is fast compared to mixing, and the reaction intensity reaches a constant:

$$I = J_A = \beta q \frac{D_m C_A^0}{d^2} Pe^{3/2}. \quad (4.16)$$

This constant is proportional to the flux of the solute across the mixing front, and corresponds to the asymptotic limit of (4.13) at large times.

The model of (4.13) is in relatively good agreement with experimental data for the range of distances and Péclet numbers considered (figure 10a). The prefactor $\beta$ in (4.7) is the unique fitting parameter of the model, and is the same for all experiment runs. As a reference, we also plot in figure 10(a) the prediction of the dispersion model, whose parameters are determined by a fit to the lowest Péclet number. Although the dispersion model captures the data relatively well for the lowest Péclet number, it diverges significantly for the largest Péclet numbers, with up to one order of magnitude difference from the data. In contrast, the mixing model captures the initial growth followed by a plateau observed for the evolution of the reaction intensity with distance (figure 10a), and the dependence of the reaction intensity and reaction maximum with the Péclet number (figure 8a1, a2). Thus it provides a mechanism to quantify incomplete mixing at the pore scale and explain the failure of the hydrodynamic dispersion model. The agreement is less good for the lowest Péclet number ($Pe = 179$), for which the model tends to overestimate reaction intensities. At low Péclet number, the fluctuation scale (set by the Batchelor scale) reaches the grain size, and incomplete mixing is expected to be less significant. Using the estimate of the Lyapunov exponent $\lambda \approx 0.2U/d$ measured in bead packs, Heyman et al. (2020) estimated that the Batchelor scale should reach the grain size for $Pe \le 5$. For the reaction fronts considered, this appears to occur at larger Péclet numbers, of the order of $Pe \le 200$, i.e. for a fluctuation scale of the order of 16% of the grain size. The complete mixing limit $Pe = 5$ was obtained assuming stretching in an unbounded domain (Heyman et al. 2020). However, the presence of grain boundaries significantly enhances diffusive homogenisation due to confinement (Hamada, Cueto-Felgueroso & de Anna 2020). Therefore, it is not surprising that the effective limit for complete mixing is larger than 5 in confined environments.

At a distance of approximately 20 grains, the reaction intensity departs from the plateau and begins to decay (figure 10a). We attribute this to the limitation of the supply of fresh reactants to the interface by dispersive mass transfer. To describe the conditions for this to occur, we define $\Delta$ as the characteristic size over which fluid deformations occur in the
transverse direction (Heyman et al. 2020). When the transverse dispersion of the front is less than $\Delta$, the stretching events at the mixing front can bring fresh fluid with $C_A = C_A^0$ from outside the mixing front. In contrast, for $x > x_\Delta$, the concentration of reactant A near the mixing interface $F = 0$ decays due to dispersion. In random bead packs, the mixing width generally follows the classical dispersive behaviour (Delgado 2007; Heyman et al. 2020; Rousseau et al. 2023; Sanquer et al. 2024):

$$s = \sqrt{4\alpha_T x}. \qquad (4.17)$$

We have measured the transverse dispersivity using a conservative tracer in the same set-up (Rousseau et al. 2023) and found $\alpha_T = 0.07d$. We define the dispersion distance as the distance at which $s = \Delta$:

$$x_\Delta = \frac{\Delta^2}{4\alpha_T}. \qquad (4.18)$$

Considering the observation $x_\Delta \approx 20d$ leads to $\Delta \approx 2.4d$, a value compatible with the typical scale at which the filaments are stretched on the pore scale (Heyman et al. 2020).

## 4.2. Stretching-enhanced reactive front kinetics under saddle flow

In the saddle-flow configuration, the experimental conditions are such that reaction is mixing-limited, $\gamma \ll t_r^{-1}$. We assume that the pore-scale mixing mechanism is similar to the one described for the co-flow situation, with a local stretching rate, and thus Péclet number, that increases with distance. Replacing $Pe$ by

$$Pe_x \sim \frac{\gamma x d}{D_m}, \qquad (4.19)$$

in (4.16), we obtain

$$I = \beta'q \frac{C_A^0}{(d D_m)^{1/2}} (\gamma x)^{3/2}. \qquad (4.20)$$

The prefactor $\beta'$ in (4.20) is the sole fitting parameter and is expected to be the same for all experiments performed in the saddle-flow configuration. The pore-scale mixing theory agrees well with the data for the distance range considered and the Péclet numbers (figure 10b). As a reference, we also plot the predictions of the dispersion model in figure 10(b). Although it matches the experimental curves at distances smaller than the pore size, i.e. for low velocities in the saddle flow, the dispersion model largely diverges from the data at larger distances for all Péclet numbers. In contrast, the mixing model captures the observed evolution of the reaction intensity with both distance (figure 10b) and Péclet number (figure 8b1). Note that, contrary to the co-flow configuration, there is no transition to the dispersive scaling at large distances. In saddle flows, the supply of reactants on both sides of the front increases linearly with distance (figure 2), allowing for a continuous supply of fresh reactants. Hence there is no limitation of reactant supply by dispersion, and the effect of pore-scale mixing persists over the whole range of experimental observations.

## 5. Conclusion

We investigated experimentally and theoretically the dynamics of a bimolecular reaction in steady mixing fronts in porous media. Two common flow configurations were considered, namely, a uniform flow and a converging flow. We first derived the effective reaction rates on the basis of the classical hydrodynamic dispersion model, which assumes well-mixed
<figure><img src="image_27.png" /><figcaption>Figure 10. Comparison of theoretical predictions for the mixing model (dashed lines, (4.13) and (4.20)) and the dispersive model (dotted lines, (2.18), (2.19) and (2.22)) with data (solid lines) for the reaction intensity versus distance for (a) the co-flow and (b) the saddle flow. Independent parameters are $Da = 1600$, $D_m = 10^{-9} \text{ m}^2 \text{ s}^{-1}$, $C_A^0 = 1 \text{ mM}$, and $q = 2$. Prefactors used for fitting models are 0.2 (dispersive co-flow), $5 \times 10^{-3}$ (dispersive saddle flow), $\beta = 2 \times 10^{-2}$ (mixing model co-flow) and $\beta' = 3 \times 10^{-5}$ (mixing model saddle flow).</figcaption></figure>

pore-scale conditions. We then experimentally measured the effective reaction rates in a cell filled with 3-D granular material. For this, we combined chemiluminescence and refractive-index matching, allowing us to precisely map the Darcy-scale reaction rates, integrated over the cell gap, across the mixing front.

Experimental reaction rates in the mixing zone show a significant departure from the predictions expected by modelling mixing with hydrodynamic dispersion. In co-flow, differences persist for approximately 20 grain diameters, while in saddle flow, they persist throughout the front. Therefore, these findings highlight the role of incomplete mixing at the pore scale in reactive fronts. Building on the recent demonstration of the chaotic nature of mixing at the pore scale in 3-D porous media, we propose a new description of micro-scale reactive mixing dynamics to explain our Darcy-scale observations. This model leads to new scaling laws for effective reaction rates in porous media, as a function of both space and Péclet number.

Although predictions accurately capture the general experimental trends, additional data are needed to uncover the details of micro-scale reactive mixing and its impact at the Darcy scale. Such data would require resolving both the pore-scale dynamics and the Darcy-scale effective reaction rates, a range of scales that is difficult to achieve experimentally with optical methods. Numerical simulations of large flow domains (Sole-Mari, Bolster & Fernández-Garcia 2022) or porous flow analogues (Heyman, Lester & Le Borgne 2021) may be able to provide such data in the future.

Finally, our findings suggest that micro-scale mixing processes in porous media have an important effect on macroscopic reaction rates, an effect that persists for long time and distances along mixing fronts. The framework proposed to capture these processes, based on the physics of pore-scale flows (Lester et al. 2013), opens new avenues for upscaling reactive transport at Darcy scale, beyond the hydrodynamic dispersion paradigm.

**Acknowledgements.** Views and opinions expressed are those of the author(s) only, and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

**Funding.** Funded by the European Union (ERC, CHORUS, 101042466).

**Declaration of interests.** The authors report no conflict of interest.
**Data availability statement.** The data and codes to reproduce the findings of this study are
openly available in Github at https://github.com/jorishey1234/reactivefronts or upon reasonable request to
joris.heyman@univrennes.fr.

**Appendix A. Theory of reactive mixing at finite Damköhler number**

A.1. Reaction–diffusion front

Gálfi & Rácz (1988) and Larralde et al. (1992) proposed an approximate solution to the equation of advection–diffusion–reaction in the case of one-dimensional fronts. The theory allows for the derivation of the asymptotic scaling with time and Damköhler number of the effective reaction front properties at the Darcy scale. Here, we recall the main steps of the derivation. We take dimensionless variables for space $y = \tilde{y}/L$ and time $t = \tilde{t}/t_D$, with $L$ a length scale, and $t_D = L^2/D_m$ the diffusion time, and dimensionless concentrations $C_i = \tilde{C}_i/C_A^0$. The ratio of initial reactant concentrations is $q = C_B^0/C_A^0$. In our experiments, $q = 2$. The dimensionless reaction–diffusion equation reads

$$
\begin{align}
\frac{\partial C_A}{\partial t} &= \frac{\partial^2 C_A}{\partial y^2} - D a C_A C_B, \nonumber \\
\frac{\partial C_B}{\partial t} &= \frac{\partial^2 C_B}{\partial y^2} - D a C_A C_B, \tag{A1}
\end{align}
$$

where $Da = t_D/t_r$ is the Damköhler number, and $t_r = 1/(kC_A^0)$ is the characteristic reaction time. It is easy to show that $F = C_A - C_B$ is a conserved quantity obeying an advection–diffusion equation. For initially segregated reactants with initial conditions $F = 1$ for $y < 0$ and $F = -q$ for $y > 0$, $F$ follows

$$
F(y, t) = \frac{1-q}{2} - \frac{1+q}{2} \operatorname{erf} \left( \frac{y}{\sqrt{4t}} \right). \quad (A2)
$$

The front position is located at $F = 0$, i.e. at

$$
y_f(t) = -\sqrt{2D_f t}, \tag{A3}
$$

where the constant $D_f$ is a solution of

$$
\mathrm{erf}\left(-\sqrt{D_f/2}\right) = (1-q)/(1+q). \qquad (\mathrm{A}4)
$$

Note that if $q=1$ (equal reactant concentrations), then $D_f=0$ and $y_f=0$ at all times.
For the chemiluminescence reaction considered, $q=2$ and $D_f \approx 0.1855$. Inserting (A2)
into (A1), a nonlinear equation is obtained for $C_A$:

$$
\frac{\partial C_A}{\partial t} = \frac{\partial^2 C_A}{\partial y^2} + D a C_A (F(y, t) - C_A). \quad (A5)
$$

Equation (A5) can be approximated if the penetration of reactants is limited to a small region near $y = y_f$, i.e. for intermediate Damköhler numbers (Gálfi & Rácz 1988; Larralde et al. 1992). Neglecting the $C_A^2$ term and taking an expansion of the error function near the front position $y_f$, we get

$$
\frac{\partial C_A}{\partial t} = \frac{\partial^2 C_A}{\partial y^2} - q' D a C_A \frac{y - y_f}{\sqrt{t}}, \quad (A6)
$$
with $q' = (1 + q) \exp(-D_f/2)/(2\sqrt{\pi})$ (Gálfi & Rácz 1988). The solution to this equation can be found by recognising the Airy differential problem on the right-hand side of (A6), such that a solution for $C_A$ is

$$C_A(y, t) \approx B t^\alpha \operatorname{Ai}\left(\frac{y - y_f}{s_r(t)}\right), \qquad (A7)$$

where $\operatorname{Ai}$ is the Airy function, and $s_r(t) = (q' Da)^{-1/3} t^{1/6}$ is a dimensionless reactive width. The dimensional reactive width reads

$$\tilde{s}_r(t) = Ls_r = L(q' Da)^{-1/3} (D_m \tilde{t}/L^2)^{1/6}. \qquad (A8)$$

Finding $\alpha$ and $B$ is possible by equating the nonlinear term with the spatial derivative in (A5), giving $\alpha = -1/3$ and $B \sim q'^{2/3} Da^{-1/3}$.

The reaction rate is also the source term of (A6), such that

$$R(y) \sim q' Da C_A \frac{y - y_f}{\sqrt{t}} = R_{max} \frac{y - y_f}{s_r(t)} \operatorname{Ai}\left(\frac{y - y_f}{s_r(t)}\right), \qquad (A9)$$

with the maximal reaction rate $R_{max} \sim q'^{4/3} Da^{1/3} t^{-2/3}$. The dimensional max reaction rate is

$$\tilde{R}_{max} = \frac{C_A^0}{t_D} R_{max} = \frac{C_A^0}{t_D^{1/3}} q'^{4/3} Da^{1/3} \tilde{t}^{-2/3} = \frac{C_A^0 D_m^{1/3}}{L^{2/3}} q'^{4/3} Da^{1/3} \tilde{t}^{-2/3}. \qquad (A10)$$

The total reaction intensity is

$$I = 2 \int_{0}^{\infty} R(y) dy \propto R_{max} s_r \sim q't^{-1/2}. \qquad (A11)$$

Its dimensional equivalent is

$$\tilde{I} = C_A^0 \frac{L}{t_D} I = C_A^0 D_m^{1/2} q' \tilde{t}^{-1/2}. \qquad (A12)$$

Note that this result can also be recovered by considering the infinitely fast reaction case, where the intensity of reaction is exactly equal to the diffusive flux across the interface:

$$\tilde{I} \sim D_m \frac{\partial C_A}{\partial y} \sim D_m \frac{C_A^0}{\sqrt{D_m \tilde{t}}}, \qquad (A13)$$

which is equivalent to (A12). In fact, the finite and infinite Damköhler cases yield equal scaling for the total reaction intensity. The difference emerges in the existence of a finite reaction width and a maximum reaction rate.

## A.2. Advection–diffusion–reaction front

In the case of an advection at constant velocity $U$, the preceding results simply hold with $\tilde{x} = \tilde{t}U$, or equivalently, $t = x Pe$. This implies that

$$s_R \sim (q' Da)^{-1/3} Pe^{-1/6} x^{1/6}, \qquad (A14)$$

$$R_{max} \sim q'^{4/3} Da^{1/3} Pe^{2/3} x^{-2/3} \qquad (A15)$$

and

$$I \sim q' Pe^{1/2} x^{-1/2}. \qquad (A16)$$
### A.3. Advection–diffusion–reaction with compression

In the presence of a constant compression rate $\gamma = U/(2L)$, the conservative profile $F$ reaches a steady state:

$$F \rightarrow \frac{1-q}{2} - \frac{1+q}{2} \operatorname{erf}\left(\frac{\tilde{y}}{s_B}\right), \qquad (A17)$$

with the Batchelor scale $s_B = \sqrt{2D_m/\gamma}$ (Villermaux 2019). Since $Pe = L^2\gamma/D_m$, this is also

$$F = \frac{1-q}{2} - \frac{1+q}{2} \operatorname{erf}(Pe^{1/2} y). \qquad (A18)$$

The front location $F = 0$ is now fixed to $y_f$, the solution of

$$\operatorname{erf}(Pe^{1/2} y_f) = \frac{1-q}{1+q}. \qquad (A19)$$

The above derivations still hold, without the time dependence. The equation for $C_A$ reads

$$0 = \frac{\partial^2 C_A}{\partial y^2} + Da C_A (F - C_A), \qquad (A20)$$

whose approximate solution is

$$C_A \sim B \operatorname{Ai}\left(\frac{y - y_f}{s_r}\right), \qquad (A21)$$

with $s_r = (q' Da)^{-1/3} Pe^{-1/6}$ and $B \sim q'^{2/3} Da^{-1/3} Pe^{1/3}$. As before, $B$ is determined by scaling the nonlinear term with the spatial derivative in (A20). Then the reaction rate is

$$R \sim R_{max} \frac{y}{s_r} \operatorname{Ai}\left(\frac{y}{s_r}\right), \qquad (A22)$$

with

$$R_{max} \sim q'^{4/3} Da^{1/3} Pe^{2/3}. \qquad (A23)$$

The total reaction intensity is

$$I \sim s_r R_{max} = q' Pe^{1/2}. \qquad (A24)$$

In dimensional form, these read

$$\begin{gathered} \tilde{s}_R \sim L(q Da)^{-1/3} Pe^{-1/6}, \quad \tilde{R}_{max} \sim C_0 D_m L^{-2} q'^{4/3} Da^{1/3} Pe^{2/3} \\ \text{and} \quad \tilde{I} \sim q' C_0 D_m L^{-1} Pe^{1/2}. \end{gathered} \qquad (A25)$$

### A.4. Numerical simulation of dispersion–reaction fronts

We validated the theoretical scaling of reactive fronts under dispersion in uniform flow (2.18)–(2.22) against numerical simulations of the full 2-D advection–dispersion–reaction problem (2.1). To numerically solve the equations, we use the OpenFOAM finite-volume software (https://openfoam.org). The domain is a 2-D rectangle; for co-flow, $x \in [0, 150]$ mm with 300 grid cells, $y \in [-25, 25]$ mm with 400 cells, and for saddle flow, $x \in [-150, 150]$ mm with 600 cells, $y \in [-25, 25]$ mm with 400 cells. We fixed $D_m = 5 \times 10^{-10} \text{ m}^2 \text{ s}^{-1}$, and varied $U$ or $\gamma$. We used an isotropic dispersivity tensor with $\alpha_T = \alpha_L = 0.06d$ m, where $d$ is the grain size (2 mm). The reaction rate constant
<figure><img src="image_28.png" /></figure>

Figure 11. Normalised reaction rate fields at Pe = 799, 1690 and 3575 for co-flow, and Pe = 1941, 4104 and 8678 for saddle flow (lowest at the top row, highest at the bottom row) from simulations for porous media. The reaction rate is normalised by the maximum reaction rate at the highest Pe for each configuration. For co-flow, the left-hand edge corresponds to the start of mixing, and for saddle flow, the left-hand edge corresponds to the stagnation point. The white dashed lines in the lowest Pe in experimental images show the streamlines.

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

Figure 12. Scaling of reactive front properties with distance to injection in co-flow with hydrodynamic dispersion obtained by numerical simulations: (a) width of the reaction front, (b) maximum reaction rate, and (c) reaction intensity. Black dashed lines stand for the theoretical prediction (2.18) and (2.19).

is fixed at $k = 0.08 \text{ mM}^{-1} \text{ s}^{-1}$. We used a kinetic rate that is 10 times smaller than the experimental value such that the width of the simulated fronts is approximately the same as the one observed experimentally. This net difference highlights the fact that the porous medium efficiently spreads the reactants but without completely mixing them (retarding reaction). To match the observed reactive width with a hydrodynamic dispersion model, it is necessary to slow down the reaction kinetics so that the reactive width approaches the dispersive width. The solute concentration at the inlet boundary of the co-flow (at $x=0$) was $(C_A, C_B) = (1, 0)$ mM for $y>0$, and $(C_A, C_B) = (0, 1)$ mM for $y<0$. The inlet boundary condition of the saddle flow was $(C_A, C_B) = (1, 0)$ mM at $y=50$, and $(C_A, C_B) = (0, 1)$ mM at $y=-50$ mm. For the outlet boundaries ($x=150$ mm for co-flow, and $y=\pm50$ mm for saddle flow), we imposed a zero gradient for all species. We used the Euler method as a temporal discretisation scheme, and a linear interpolation scheme for interpolating face-centred values from cell-centred values. We varied the
<figure><img src="image_30.png" /><figcaption>Figure 13. Scaling of reactive front properties with distance to injection in saddle-flow with hydrodynamic dispersion obtained by numerical simulations: (a) width of the reaction front, (b) maximum reaction rate, and (c) reaction intensity. Black dashed lines stand for the theoretical prediction in the kinetic-limited regime (2.21).</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 14. Simulated reactive front over Pe at 150 mm in co-flow under hydrodynamic dispersion: (a) width of the reaction, (b) maximum reaction rate, and (c) reaction intensity. Black dashed lines are the hydrodynamic dispersion predictions for the mixing-limited regime at fixed time (2.19).</figcaption></figure>

velocity $U$ (co-flow) and compression rate $\gamma$ (saddle flow) in the same range as in the experiments. The characteristic length scale to define the Péclet number in the co-flow is the diameter of the grain, $L = 2$ mm. In saddle flow, the characteristic velocity is $U = \gamma L$, with the length $L = 100$ mm corresponding to our experimental set-up (see § 3.1.2). Such definitions resulted in simulations with Pe ranging from 179 to 3575 for co-flow, and from 434 to 8678 for saddle flow. The spatial behaviour of the reactive fronts in the presence of dispersion is plotted in figure 11.

### A.5. Simulation results

The dependences of the maximum reaction rate $R_{max}$, total reaction intensity $I$ and reaction width $s_R$ with position $x$ along the front in co-flow and saddle flows are plotted in figures 12 and 13 respectively, as well as the theoretical predictions (2.18)–(2.22). The predictions show very good agreement with the simulations. Next, the front properties are plotted against Pe at a given location $x$ in figures 14 and 15 for the co-flow and saddle
<figure><img src="image_32.png" /></figure>

Figure 15. Simulated reactive front over Pe at 150 mm in saddle flow under hydrodynamic dispersion: (a) width of the reaction, (b) maximum reaction rate, and (c) reaction intensity. The green dashed lines are hydrodynamic dispersion predictions of reaction limited regime (2.21), whereas the black dashed lines are those of mixing-limited regime (2.22). The blue dashed lines show fitted lines with fitted exponent.

flows, respectively. In co-flow, the theoretical scaling laws compare well with simulations. On the other hand, in saddle flow, the scaling laws of the mixing-limited regime are verified only for low Pe. This is because the longitudinal derivatives that were dropped in (2.12) $(\gamma x \partial C_A/\partial x$ and $\alpha_L \gamma |y| \partial^2 C_A/\partial x^2$) become large compared to $kC_A C_B$, and the approximation becomes not valid at high Pe. Thus the approximation is valid for low to moderate Pe. To compare with the experiments, the scaling exponents were derived by fitting the simulation results.

REFERENCES

AQUINO, T., LE BORGNE, T. & HEYMAN, J. 2023 Fluid-solid reaction in porous media as a chaotic restart process. *Phys. Rev. Lett.* **130** (26), 264001.

BANDOPADHYAY, A., LE BORGNE, T., MÉHEUST, Y. & DENTZ, M. 2017 Enhanced reaction kinetics and reactive mixing scale dynamics in mixing fronts under shear flow for arbitrary Damköhler numbers. *Adv. Water Resour.* **100**, 78–95.

BATCHELOR, G.K., HOWELLS, I.D. & TOWNSEND, A.A. 1959 Small-scale variation of convected quantities like temperature in turbulent fluid. Part 2. The case of large conductivity. *J. Fluid Mech.* **5** (1), 134.

BATTIATO, I. & TARTAKOVSKY, D.M. 2011 Applicability regimes for macroscopic models of reactive transport in porous media. *J. Contam. Hydrol.* **120–121** (C), 18–26.

BATTIATO, I., TARTAKOVSKY, D.M., TARTAKOVSKY, A.M. & SCHEIBE, T. 2009 On breakdown of macroscopic models of mixing-controlled heterogeneous reactions in porous media. *Adv. Water Resour.* **32** (11), 1664–1673.

BEAR, J. 1988 *Dynamics of Fluids in Porous Media*. Dover Publications.

BENSON, D.A. & MEERSCHAERT, M.M. 2008 Simulation of chemical reaction via particle tracking: diffusion-limited versus thermodynamic rate-limited regimes. *Water Resour. Res.* **44** (12), W12201.

BOLSTER, D., PASTER, A. & BENSON, D.A. 2016 A particle number conserving Lagrangian method for mixing-driven reactive transport. *Water Resour. Res.* **52** (2), 1518–1527.

BRESCIANI, E., KANG, P.K. & LEE, S. 2019 Theoretical analysis of groundwater flow patterns near stagnation points. *Water Resour. Res.* **55** (2), 1624–1650.

CHIOGNA, G. & BELLIN, A. 2013 Analytical solution for reactive solute transport considering incomplete mixing within a reference elementary volume. *Water Resour. Res.* **49** (5), 2589–2600.

DE ANNA, P., JIMENEZ-MARTINEZ, J., TABUTEAU, H., TURUBAN, R., LE BORGNE, T., DERRIEN, M. & MÉHEUST, Y. 2014 Mixing and reaction kinetics in porous media: an experimental pore scale quantification. *Environ. Sci. Technol.* **48** (1), 508–516.

DELGADO, J.M.P.Q. 2007 Longitudinal and transverse dispersion in porous media. *Chem. Engng Res. Des.* **85** (9), 1245–1252.
DENTZ, M., LE BORGNE, T., ENGLERT, A. & BIJELJIC, B. 2011 Mixing, spreading and reaction in heterogeneous media: a brief review. *J. Contam. Hydrol.* **120–121** (C), 1–17.

DING, D., BENSON, D.A., PASTER, A. & BOLSTER, D. 2013 Modeling bimolecular reactions and transport in porous media via particle tracking. *Adv. Water Resour.* **53**, 56–65.

EDERY, Y., DROR, I., SCHER, H. & BERKOWITZ, B. 2015 Anomalous reactive transport in porous media: experiments and modeling. *Phys. Rev. E* **91** (5), 052130.

EDERY, Y., SCHER, H. & BERKOWITZ, B. 2010 Particle tracking model of bimolecular reactive transport in porous media. *Water Resour. Res.* **46** (7), W12201.

GÁLFI, L. & RÁCZ, Z. 1988 Properties of the reaction front in an $A + B \rightarrow C$ type reaction-diffusion process. *Phys. Rev. A* **38** (6), 3151-3154.

GINN, T.R. 2018 Modeling bimolecular reactive transport with mixing-limitation: theory and application to column experiments. *Water Resour. Res.* **54** (1), 256–270.

GOIRAND, F., LE BORGNE, T. & LORTHOIS, S. 2021 Network-driven anomalous transport is a fundamental component of brain microvascular dysfunction. *Nat. Commun.* **12** (1), 7295.

GRAMLING, C.M., HARVEY, C.F. & MEIGS, L.C. 2002 Reactive transport in porous media: a comparison of model prediction with laboratory visualization. *Environ. Sci. Technol.* **36** (11), 2508–2514.

HALLER, G. 2015 Lagrangian coherent structures. *Annu. Rev. Fluid Mech.* **47** (1), 137–162.

HAM, P.A.S., PROMMER, H., OLSSON, ÅH., SCHOTTING, R.J. & GRATHWOHL, P. 2007 Predictive modelling of dispersion controlled reactive plumes at the laboratory-scale. *J. Contam. Hydrol.* **93** (1–4), 304–315.

HAMADA, M., CUETO-FELGUEROSO, L. & DE ANNA, P. 2020 Diffusion limited mixing in confined media. *Phys. Rev. Fluids* **5** (12), 124502.

HESTER, E.T., CARDENAS, M.B., HAGGERTY, R. & APTE, S.V. 2017 The importance and challenge of hyporheic mixing. *Water Resour. Res.* **53** (5), 3565–3575.

HEYMAN, J., LESTER, D.R. & LE BORGNE, T. 2021 Scalar signatures of chaotic mixing in porous media. *Phys. Rev. Lett.* **126** (3), 034505.

HEYMAN, J., LESTER, D.R., TURUBAN, R., MÉHEUST, Y. & LE BORGNE, T. 2020 Stretching and folding sustain microscale chemical gradients in porous media. *Proc. Natl Acad. Sci. USA* **117** (24), 13359–13365.

HOCHSTETLER, D.L. & KITANIDIS, P.K. 2013 The behavior of effective rate constants for bimolecular reactions in an asymptotic transport regime. *J. Contam. Hydrol.* **144** (1), 88–98.

HOCHSTETLER, D.L., ROLLE, M., CHIOGNA, G., HABERER, C.M., GRATHWOHL, P. & KITANIDIS, P.K. 2013 Effects of compound-specific transverse mixing on steady-state reactive plumes: insights from pore-scale simulations and Darcy-scale experiments. *Adv. Water Resour.* **54**, 1–10.

HUANG, G., WEI, X., GU, Y., KANG, Z., LAO, L., LI, L., FAN, J. & SHOU, D. 2022 Heterogeneously engineered porous media for directional and asymmetric liquid transport. *Cell Rep. Phys. Sci.* **3** (1), 1–25.

IZUMOTO, S., HEYMAN, J., HUISMAN, J.A., DE VRIENDT, K., SOULAINE, C., GOMEZ, F., TABUTEAU, H., MÉHEUST, Y. & LE BORGNE, T. 2023 Enhanced mixing and reaction in converging flows: theory and pore-scale imaging. *Water Resour. Res.* **59** (8), e2023WR034749.

KHALED, A.-R.A. & VAFAI, K. 2003 The role of porous media in modeling flow and heat transfer in biological tissues. *Intl J. Heat Mass Transfer* **46** (26), 4989–5003.

KLENK, I.D. & GRATHWOHL, P. 2002 Transverse vertical dispersion in groundwater and the capillary fringe. *J. Contam. Hydrol.* **58** (1–2), 111–128.

LARRALDE, H., ARAUJO, M., HAVLIN, S. & STANLEY, H.E. 1992 Reaction front for $A + B \rightarrow C$ diffusion-reaction systems with initially separated reactants. *Phys. Rev. A* **46** (2), 855–859.

LE BORGNE, T., DENTZ, M. & VILLERMAUX, E. 2013 Stretching, coalescence, and mixing in porous media. *Phys. Rev. Lett.* **110** (20), 1–5.

LE BORGNE, T., DENTZ, M. & VILLERMAUX, E. 2015 The lamellar description of mixing in porous media. *J. Fluid Mech.* **770**, 458–498.

LESTER, D.R., DENTZ, M., LE BORGNE, T. & DE BARROS, F.P.J. 2018 Fluid deformation in random steady three-dimensional flow. *J. Fluid Mech.* **855**, 770–803.

LESTER, D.R., METCALFE, G. & TREFRY, M.G. 2013 Is chaotic advection inherent to porous media flow? *Phys. Rev. Lett.* **111** (17), 174101.

LESTER, D.R., TREFRY, M.G. & METCALFE, G. 2016 Chaotic advection at the pore scale: mechanisms, upscaling and implications for macroscopic transport. *Adv. Water Resour.* **97**, 175–192.

LICHTNER, P.C., STEEFEL, C.I. & OELKERS, E.H. 2018 *Reactive Transport in Porous Media*. vol. 34 Walter de Gruyter.

MARKALE, I., CIMMARUSTI, G.M., BRITTON, M.M. & JIMÉNEZ-MARTÍNEZ, J. 2021 Phase saturation control on mixing-driven reactions in 3D porous media. *Environ. Sci. Technol.* **55** (13), 8742–8752.
MARZADRI, A., DEE, M.M., TONINA, D., BELLIN, A. & TANK, J.L. 2017 Role of surface and subsurface processes in scaling N₂O emissions along riverine networks. *Proc. Natl Acad. Sci. USA* **114** (17), 4330–4335.

MATSUMOTO, R. & MATSUO, R. 2015 Measurement of local chemical reaction rate in a microchannel by using luminol chemiluminescence. *Expl Therm. Fluid Sci.* **67**, 88–95.

OLSSON, A. & GRATHWOHL, P. 2007 Transverse dispersion of non-reactive tracers in porous media: a new nonlinear relationship to predict dispersion coefficients. *J. Contam. Hydrol.* **92** (3–4), 149–161.

PASTER, A., BOLSTER, D. & BENSON, D.A. 2014 Connecting the dots: semi-analytical and random walk numerical solutions of the diffusion–reaction equation with stochastic initial conditions. *J. Comput. Phys.* **263**, 91–112.

PERFECT, E., CHENG, C.L., KANG, M., BILHEUX, H.Z., LAMANNA, J.M., GRAGG, M.J. & WRIGHT, D.M. 2014 Neutron imaging of hydrogen-rich fluids in geomaterials and engineered porous media: a review. *Earth-Sci. Rev.* **129**, 120–135.

PORTA, G.M., CHAYNIKOV, S., THOVERT, J., RIVA, M., GUADAGNINI, A. & ADLER, P.M. 2013 Numerical investigation of pore and continuum scale formulations of bimolecular reactive transport in porous media. *Adv. Water Resour.* **62**, 243–253.

RAJE, D.S. & KAPOOR, V. 2000 Experimental study of bimolecular reaction kinetics in porous media. *Environ. Sci. Technol.* **34** (7), 1234–1239.

RANZ, W.E. 1979 Applications of a stretch model to mixing, diffusion, and reaction in laminar and turbulent flows. *AIChE J.* **25** (1), 41–47.

ROLLE, M. & LE BORGNE, T. 2019 Mixing and reactive fronts in the subsurface. *Rev. Mineral. Geochem.* **85** (1), 111–142.

ROUSSEAU, G., IZUMOTO, S., LE BORGNE, T. & HEYMAN, J. 2023 Dispersion versus diffusion in mixing fronts. *Water Resour. Res.* **59** (11), e2023WR035848.

SANCHEZ-VILA, X., FERNÁNDEZ-GARCIA, D. & GUADAGNINI, A. 2010 Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation. *Water Resour. Res.* **46** (12), W12510.

SANQUER, H., HEYMAN, J., HANNA, K. & LE BORGNE, T. 2024 Microscale chaotic mixing as a driver for chemical reactions in porous media. *Environ. Sci. Technol.* **58** (20), 8899–8908.

SOLE-MARI, G., BOLSTER, D. & FERNÁNDEZ-GARCIA, D. 2022 A closer look: high-resolution pore-scale simulations of solute transport and mixing through porous media columns. *Transp. Porous Med.* **146**, 1–27.

SOUZY, M., LHUISSIER, H., MÉHEUST, Y., LE BORGNE, T. & METZGER, B. 2020 Velocity distributions, dispersion and stretching in three-dimensional porous media. *J. Fluid Mech.* **891**, A16.

TAITELBAUM, H., HAVLIN, S., KIEFER, J.E., TRUS, B. & WEISS, G.H. 1991 Some properties of the A + B → C reaction–diffusion system with initially separated components. *J. Stat. Phys.* **65** (5–6), 873–891.

UCHIDA, S., SATOH, Y., YAMASHIRO, N. & SATOH, T. 2004 Determination of hydrogen peroxide in water by chemiluminescence detection, (II). *J. Nucl. Sci. Technol.* **41** (9), 898–906.

VALOCCHI, A.J., BOLSTER, D. & WERTH, C.J. 2019 Mixing-limited reactions in porous media. *Transp. Porous Med.* **130** (1), 157–182.

VILLERMAUX, E. 2018 On shapes and forms: population balance dynamics of corrugated stirred fronts. *C. R. Phys.* **19** (5), 306–315.

VILLERMAUX, E. 2019 Mixing versus stirring. *Annu. Rev. Fluid Mech.* **51** (1), 245–273.

DE VRIENDT, K. 2021 Mixing and chemical reaction hotspots in saline–freshwater mixing zones. PhD thesis. Universitat Politècnica de Catalunya.

WILLINGHAM, T.W., WERTH, C.J. & VALOCCHI, A.J. 2008 Evaluation of the effects of porous media structure on mixing-controlled reactions using pore-scale modeling and micromodel experiments. *Environ. Sci. Technol.* **42** (9), 3185–3193.