# Enhanced Mixing and Reaction in Converging Flows: Theory and Pore-Scale Imaging

Satoshi Izumoto, Joris Heyman, Johan Alexander Huisman, Kevin de Vriendt, Cyprien Soulaine, Francesco Gomez, Hervé Tabuteau, Yves Méheust, Tanguy Le Borgne

► To cite this version:

Satoshi Izumoto, Joris Heyman, Johan Alexander Huisman, Kevin de Vriendt, Cyprien Soulaine, et al. Enhanced Mixing and Reaction in Converging Flows: Theory and Pore-Scale Imaging. Water Resources Research, 2023, 59 (8), pp.e2023WR034749. 10.1029/2023wr034749. insu-04182342

HAL Id: insu-04182342

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

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_1.png" /></figure>
# Water Resources Research®

## RESEARCH ARTICLE

10.1029/2023WR034749

**Key Points:**

* We developed analytical solutions for mixing-induced reaction rates in converging flows for any Péclet and Damköhler numbers

* We measured pore scale reaction rates in mixing fronts compressed at constant rate by converging flows in millifluidic experiments

* The experiments validated the theory, showing a four times larger reaction rate in porous media compared to Hele-Shaw cells

**Supporting Information:**

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

**Correspondence to:**

S. Izumoto,
satoshi.izumoto@hotmail.co.jp

**Citation:**

Izumoto, S., Heyman, J., Huisman, J. A., De Vriendt, K., Soulaine, C., Gomez, F., et al. (2023). Enhanced mixing and reaction in converging flows: Theory and pore-scale imaging. *Water Resources Research*, *59*, e2023WR034749. https://doi.org/10.1029/2023WR034749

Received 24 FEB 2023

Accepted 18 JUL 2023

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

## Enhanced Mixing and Reaction in Converging Flows: Theory and Pore-Scale Imaging

Satoshi Izumoto¹,² ID, Joris Heyman³ ID, Johan Alexander Huisman² ID, Kevin De Vriendt⁴ ID, Cyprien Soulaine⁵ ID, Francesco Gomez³, Hervé Tabuteau¹, Yves Méheust³ ID, and Tanguy Le Borgne³ ID

¹Institut de Physique de Rennes, University Rennes 1, CNRS, Unité Mixte de Recherche 6118, Rennes, France, ²Agrosphere (IBG-3), Institute of Bio- and Geosciences, Forschungszentrum Jülich, Jülich, Germany, ³University Rennes 1, CNRS, Géosciences Rennes, UMR 6118, Rennes, France, ⁴Institute of Environmental Assessment and Water Research, IDAEA-CSIC, Barcelona, Spain, ⁵Earth Sciences Institute of Orléans, CNRS, Université d'Orléans, BRGM, Orléans, France

**Abstract** Mixing fronts at the interface of opposing flows are compressed at a constant rate. The resulting exponential stretching of fluid elements leads to enhanced chemical gradients and biogeochemical processes. This process is similar as what occurs in the pore space of 3D chaotic flows. However, it is so far not known how such fluid compression controls the amplitude of mixing and reaction rates in porous media. Here we derive analytical predictions for the mixing width, the maximum reaction rate and the reaction intensity in compressed mixing fronts as a function of the Péclet and Damköhler numbers. We developed an experimental setup providing pore scale measurements of mixing and reaction rates in mixing fronts at the interface of converging flows. The theory accurately predicts the scaling of mixing and reaction with the Péclet number both in porous micromodels and simple Hele-Shaw cells. Additionally, we found that the presence of pore scale heterogeneities in the porous micromodels enhances reaction rates by a factor of 4 compared to the Hele-Shaw cells. Using numerical simulations of pore scale velocity fields, we attributed this phenomenon to the enhancement of pore-scale compression due to the presence of grains in accelerating flows. These findings provide new insights into the dynamics of mixing-induced reactions in porous media.

### 1. Introduction

Mixing-induced reactions in porous media affect the transport, transformation and degradation of chemical and biological substances in subsurface environments (Dentz et al., 2011, 2023; Rolle & Le Borgne, 2019; Valocchi et al., 2019). Hence, they play a central role in a large range of applications such as hyporheic zone processes (Bandopadhyay et al., 2018; Hester et al., 2017; Ziliotto et al., 2021), freshwater-saltwater mixing (De Vriendt et al., 2022), soil and groundwater remediation (Fu et al., 2014; Kitanidis & McCarty, 2012; Wang et al., 2022), geological carbon sequestration (Szulczewski et al., 2012; Zoback & Gorelick, 2012), hydrogen storage (Lysyy et al., 2022; Tarkowski, 2019) and geothermal systems (Burté et al., 2019). Reactions in porous media often occur in incompletely mixed conditions, where reactants are partly segregated and reaction rates are limited by mixing (de Anna, Dentz, et al., 2014; Gramling et al., 2002; Sole-Mari et al., 2020). Therefore, many biogeochemical reactions develop in localized mixing hot spots (Bochet et al., 2020; Briggs et al., 2014; Gu et al., 2012; McClain et al., 2003).

Under non-uniform flows, velocity gradients deform fluid elements through shear and compression. This deformation enhances chemical gradients and mixing-induced reactions (Benson et al., 2017; Bolster et al., 2011; de Anna, Dentz, et al., 2014; de Anna et al., 2013; De Barros et al., 2012; Engdahl et al., 2014; Hochstetler & Kitanidis, 2013; Le Borgne et al., 2014, 2015; Lee & Kang, 2020; Soler-Sagarra et al., 2023; Willingham et al., 2008; Ye et al., 2015). The effect of shear on enhanced reaction has been studied by Bandopadhyay et al. (2017), where fluid elements deform linearly in time. Mixing fronts at opposing flows are subject to constant compression rates. As a fluid element travels toward/away from the resulting stagnation points, it deforms exponentially over time (De Vriendt et al., 2022; Hester et al., 2017; Hidalgo et al., 2015). Thus, reactions are strongly enhanced both at the stagnation point and downstream where fluid elements experience constant acceleration. These stagnation points and associate reaction hot spots develop in non-uniform flows (Bresciani et al., 2019), as observed in the hyporheic zone (Hester et al., 2017; Marzadri et al., 2017), in fresh-salt water interfaces (De Vriendt et al., 2022) and in density-driven flows (Hidalgo et al., 2015). In addition, the exponential elongation/compression can be created by engineered injection and pumping under chaotic stretching and folding protocols
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic illustration showing (a) flow around the recharging well in regional flow field. The black lines are the streamlines, blue arrows indicate the flow directions and the orange circles indicate the positions of the stagnation point and the well. (b) The deformation of small fluid elements near the stagnation point where blue and orange represent different solutes. (c) The concentration profile of species A (orange) and B (blue) near the stagnation point. The black line represents the conservative component F. The green line represents the linear approximation of the error function (F) near the stagnation point.</figcaption></figure>

(Mays & Neupauer, 2012; Piscopo et al., 2013; Wang et al., 2022). At the pore-scale, in three dimensional porous media, similar chaotic mixing dynamics are naturally generated by branching and merging at pore junctions as well as by twisting and braiding of streamlines (Heyman et al., 2020; Lester et al., 2013, 2014, 2016; Souzy et al., 2020; Turuban et al., 2018, 2019; Ye et al., 2015). Exponential elongation of mixing interfaces hence lead to a large enhancement of mixing-induced reaction either locally in converging flows or globally in chaotic flows. Yet, it is not known how exponential elongation controls effective mixing and reaction rates over the range of flow rates and reaction kinetics relevant to natural systems. Here, we address this question by theoretical, experimental, and numerical analysis of a reactive mixing front at the interface of converging flows. We focus on bimolecular reactions, which represent the prototype of mixing-induced reaction and are often studied to help understanding the control of fluid mixing on more complex reactions such as redox reactions, precipitation and dissolution (Rolle & Le Borgne, 2019; Valocchi et al., 2019). We derive scaling laws for the width, maximum reaction rate and reaction intensity as a function of Péclet (Pe) and Damköhler (Da) numbers. We compare these predictions to conservative and reactive transport experiments in mixing fronts sustained by opposing flows in a 2D Hele-Shaw cell apparatus, in the absence and in the presence of porous media. Such 2D millifluidic setup allows a precise quantification of mixing and reaction at the pore scale. We therefore investigate the coupling between Darcy scale flow heterogeneity induced by converging flows and pore scale mixing and reaction. We finally present pore scale simulations to explain these results in terms of pore scale fluid stretching and compression.

## 2. Theory

### 2.1. Conservative Mixing Near Stagnation Points in Pore Scale

To create a mixing front at the interface of opposing flows, we create a point injection in a cell where we impose a background flow (Figure 1a). This results in the formation of a stagnation point where the opposing flows cancel each other. In fluid mechanics such superposition of uniform flow and source flow is known as the Rankine half body flow. In hydrogeology, this corresponds to a common situation where fluid injected into a well under regional flow (Figure 1a). Near the stagnation point, the flow field ($v_x, v_y$) can be linearized as:

$$v_x = -\gamma x, v_y = \gamma y \qquad (1)$$

where $\gamma$ is the local compression rate (Tennekes & Lumley, 1972). Owing to fluid incompressibility, this compression rate is equal to the acceleration of fluid elements as they move away from the stagnation point in the mixing front. As discussed in the following, this acceleration is maintained on average over a large part of the mixing front.

The interplay between fluid compression and diffusion can be quantified in the Lagrangian coordinate system formed by elongated solute structures, called lamellae (Ranz, 1979; Villermaux, 2019). The so-called *lamellar*
theory of mixing provides accurate predictions of mixing rates in a range of heterogeneous flows (Heyman et al., 2020; Le Borgne et al., 2015). Here we consider the advection-diffusion problem under exponential compression. The size $\tilde{\delta}$ of an element of fluid (Figure 1b), is therefore compressed along the direction $x$ as (Figure S5 in Supporting Information S1):

$$ \gamma = -\frac{1}{\tilde{\delta}} \frac{d\tilde{\delta}}{dt} \qquad (2) $$

leading to

$$ \tilde{\delta} = \tilde{\delta}_0 e^{-\gamma t} \qquad (3) $$

Such exponential compression results from a constant compression rate $\gamma$, which develops both at stagnation points and in chaotic flows (Heyman et al., 2020). Fluid compression increases concentration gradients, thus competing with molecular diffusion, which tends to broaden solute interfaces (Le Borgne et al., 2015; Villermaux, 2019). The advective compression and diffusive expansion rates balance when $\gamma \approx D/\tilde{\delta}^2$, where $D$ is the diffusion coefficient. This leads to an asymptotically constant mixing width, called the Batchelor scale $s_B$ (Batchelor et al., 1959; Villermaux, 2019):

$$ s_B = \sqrt{\frac{2D}{\gamma}} \qquad (4) $$

where $D$ is the diffusion coefficient. Here we included a prefactor $\sqrt{2}$ to ensure that $s_B$ is equal to one standard deviation of the concentration distribution. For a characteristic length scale $L$, the Péclet number can be defined as the ratio of the diffusion time $L^2/2D$ over the stretching time $\gamma^{-1}$ (Villermaux, 2019), leading to,

$$ Pe = \gamma L^2 / (2D) = (L / s_B)^2. \qquad (5) $$

Considering a mixing front with one conservative species on one side of the front, the concentration gradients along the elongation direction can be assumed to be negligible compared to those developing in the compression direction (Villermaux, 2019). Therefore, the transport problem can be formulated as a one dimensional transport equation along the compression direction $x$ (Figure 1):

$$ \frac{\partial C}{\partial t} - \gamma x \frac{\partial C}{\partial x} = D \frac{\partial^2 C}{\partial x^2}. \qquad (6) $$

This equation can be transformed to a diffusion equation using the following change of variables (Ranz, 1979):

$$ \theta(t) = \int_{0}^{t} d\tau D/\tilde{\delta}(\tau)^2 = \frac{D}{2\gamma\tilde{\delta}_0^2}(e^{2\gamma t} - 1), \qquad (7) \\ \tilde{x} = x/L\delta $$

where $\theta$ is called warped time. In these non-dimensional variables, Equation 6 reduces to:

$$ \frac{L^2}{\tilde{\delta}^2} \frac{\partial C}{\partial \theta} = \frac{\partial^2 C}{\partial \tilde{x}^2} \qquad (8) $$

For the considered conservative front, this equation has the solution,

$$ C(\tilde{x}, \theta) = \frac{1}{2} \left( 1 + \operatorname{erf} \left( \tilde{x} \frac{L}{\tilde{\delta} \sqrt{4\theta}} \right) \right), \qquad (9) $$

leading to the steady state solution (when $\exp(2\gamma t) \gg 1$),

$$ C(x) = \frac{1}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{s_B} \right) \right). \qquad (10) $$

This therefore implies that the fronts should stabilize asymptotically to a steady state width equal to the Batchelor scale, balancing the competing effects of compression and diffusion.
## 2.2. Reactive Mixing in Compressed Mixing Fronts at the Pore Scale

Analytical predictions of reactive mixing in linear shear flows were obtained by coupling the lamellar description with bi-molecular chemical reactions (Bandopadhyay et al., 2017; Larralde et al., 1992). Here we extend this theory for exponential compression, which is relevant to mixing fronts at converging flows and to chaotic flows. As for the conservative case above, the reactive transport equation for a bimolecular reaction $A + B \rightarrow C$ is expressed in the compression direction $x$ as (Bandopadhyay et al., 2017; Villermaux, 2019):

$$ \frac{\partial C_A}{\partial t} - \gamma x \frac{\partial C_A}{\partial x} = D \frac{\partial^2 C_A}{\partial x^2} - k C_A C_B \quad (11) $$

where $C_A$ and $C_B$ are concentrations of reactants normalized by the initial concentration $C_0$, $k$ is the reaction constant for dimensionless concentrations defined by $k = \tilde{k}C_0$, with the rate constant $\tilde{k}$ in units of the inverse of the product of concentration by time. Following the same change of variable as above (Equation 7), this equation can be transformed to a diffusion-reaction equation:

$$ \frac{L^2}{\delta^2} \frac{\partial C_A}{\partial \theta} = \frac{\partial^2 C_A}{\partial \tilde{x}^2} - D A C_A C_B \quad (12) $$

with

$$ D A = k L^2 / D \quad (13) $$

the Damköhler number. We define the normalized reaction rate $R(\tilde{x})$ and the normalized total reaction intensity $I$ as,

$$ R(\tilde{x}) = D A C_A C_B, \quad (14) $$

and the integrated reaction rate over the line across the mixing front,

$$ I = \int_{-\infty}^{+\infty} R(\tilde{x}) d\tilde{x}. \quad (15) $$

### 2.2.1. Mixing-Limited Regime Da < Pe

When $Da < Pe$, the reaction time is large compared to the compression time. The interpenetration of the reactants in the mixing zone is important since the decrease of the concentrations is slow (Figure 1c). In this regime, the concentration profiles for $C_A$ and $C_B$ are expected to be close to the conservative profiles (Equation 10):

$$ \begin{aligned} C_A &\approx \frac{1}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{s_B} \right) \right), \\ C_B &\approx \frac{1}{2} \left( 1 - \operatorname{erf} \left( \frac{x}{s_B} \right) \right) \end{aligned} \quad (16) $$

This leads to,

$$ R = \frac{1}{4} D A \left( 1 - \left[ \operatorname{erf} \left( \frac{x}{s_B} \right) \right]^2 \right) \quad (17) $$

Hence, the reaction width scales as the mixing scale $s_B$,

$$ s_R \sim s_B \sim Pe^{-1/2}, \quad (18) $$

and the maximum reaction rate and reaction intensity scale respectively as:

$$ R_{\max} \sim D A \quad (19) $$

and

$$ I \sim D A P e^{-1/2} \quad (20) $$
### 2.2.2. Reaction Limited Regime Da > Pe

For fast reaction (Da > Pe), the solution of Equation 12 can be approximated by a perturbation technique (Bandopadhyay et al., 2017; Larralde et al., 1992). The evolution equation for reactants $C_A$ and $C_B$ can be separated into a conserved quantity $F = C_A - C_B$ and a reactive quantity $g$ (Figure 1c). For identical initial concentrations of $C_A$ and $C_B$, this leads to:

$$ 
\begin{aligned}
C_A(x,t) &= F(x,t)H(x) + g(|x|,t) \\
C_B(x,t) &= F(-x,t)H(-x) + g(|x|,t),
\end{aligned}
\quad (21) $$

where $H$ is the Heaviside function. Inserting Equation 21 in Equation 13, the conservative variable $F$ follows the diffusion equation (Equation 12). Subjected to the boundary conditions of the mixing front $F(-\infty, \theta) = -F(\infty, \theta) = 1$, this equation has the solution

$$ F(\tilde{x}, \theta) = \operatorname{erf}\left(\tilde{x} \frac{L}{\delta \sqrt{4\theta}}\right), \quad (22) $$

which leads to the stationary solution

$$ F(\tilde{x}) = \operatorname{erf}(\tilde{x} \operatorname{Pe}^{1/2}), \quad (23) $$

where we have used $L/s_B = \operatorname{Pe}^{1/2}$. Inserting Equations 21 and 23 into Equation 12 and considering the stationary solution, we obtain a non-linear equation for the quantity $g$ (Larralde et al., 1992):

$$ \frac{\partial^2 g}{\partial \tilde{x}^2} = \operatorname{Da} g (\operatorname{erf}(\tilde{x} \operatorname{Pe}^{1/2}) + g) \quad (24) $$

with $\tilde{x} > 0$.

The case of instantaneous reaction (Da → ∞) is obtained by setting $g = 0$. In this limit, reactants do not overlap except at a middle point, and their concentrations follow

$$ 
\begin{aligned}
C_A &= F(\tilde{x})H(\tilde{x}) \\
C_B &= F(-\tilde{x})H(-\tilde{x}).
\end{aligned}
\quad (25) $$

In this case, the dimensionless reaction intensity $I_\infty$ is equal to the mass consumption of species $C_A$ or $C_B$ at $\tilde{x} = 0$ through the reaction $A + B \rightarrow C$. Thus, $I_\infty$ can be calculated from the diffusive flux at $\tilde{x} = 0$

$$ I_{\infty} = \left. \frac{\partial C_A}{\partial \tilde{x}} \right|_{\tilde{x}=0} = \frac{2}{\sqrt{\pi}} \operatorname{Pe}^{1/2}. \quad (26) $$

In the more general case where Da is finite, Equation 24 can be solved analytically if the interpenetration of reactants is limited to a small region near $x = 0$ (Figure 1c). In this case, an approximation of Equation 24 can be obtained by neglecting the $g^2$ term and by taking a linear approximation of the error function close to the origin (Larralde et al., 1992) (Figure 1c). This leads to the approximation:

$$ \frac{\partial^2 g}{\partial \tilde{x}^2} \approx \operatorname{Da} g \frac{\tilde{x} L}{s_B} = \operatorname{Da} \operatorname{Pe}^{1/2} g \tilde{x}. \quad (27) $$

This is an Airy differential equation, whose solution is

$$ g = C \operatorname{Ai}((\operatorname{Da}\operatorname{Pe}^{1/2})^{1/3} \tilde{x}), \quad (28) $$

where $C$ is a constant. Considering the stationary solution of Equation 12, the reaction rate is

$$ R(\tilde{x}) = \operatorname{Da} C_A C_B = \frac{\partial^2 C_A}{\partial \tilde{x}^2} = \frac{\partial^2 g}{\partial \tilde{x}^2}. \quad (29) $$

Since the solution is symmetric around $\tilde{x} = 0$, the reaction intensity is
$$I = 2 \int_{0}^{\infty} \frac{\partial^2 g}{\partial \tilde{x}^2} d\tilde{x} = 2 \left[ \frac{\partial g}{\partial \tilde{x}} \right]_{0}^{+\infty} = -2C (\text{DaPe}^{1/2})^{1/3} \text{Ai}'(0), \quad (30)$$

where $Ai'(0) = -1/(3^{1/3}\Gamma(1/3)) \approx -0.26$. Since Equation 30 is valid for arbitrary large Da, we equate Equations 26 and 30, which determines the constant C as:

$$C \approx \frac{1}{0.26\sqrt{\pi}} \text{Da}^{-1/3} \text{Pe}^{1/3}. \quad (31)$$

Therefore, the reaction intensity at finite Da is approximately the same as that at the infinite Da:

$$I = \frac{2}{\sqrt{\pi}} \text{Pe}^{1/2}, \quad (32)$$

At finite Da, the reactant interpenetrate more than at infinite Da as the reaction is slower. However, the integrated reaction rate is still equal to the flux of molecules A toward to the B domain, which ultimately all react. The local reaction rate $R(\tilde{x})$ is derived by inserting Equation 28 into Equation 29,

$$R(\tilde{x}) = \text{DaPe}^{1/2} g \tilde{x} = \text{Da}^{2/3} \text{Pe}^{5/6} \text{Ai}(\text{Da}^{1/3} \text{Pe}^{1/6} \tilde{x}), \quad (33)$$

which can be formulated as,

$$R(\tilde{x}) = R_{\max} \text{Ai} \left( \frac{\tilde{x}}{s_R} \right) \frac{\tilde{x}}{s_R}, \quad (34)$$

with the dimensionless reactive scale $s_R$ and the maximum reaction rate $R_{\max}$ defined as,

$$s_R = \text{Da}^{-1/3} \text{Pe}^{-1/6}, \quad (35)$$

$$R_{\max} = \frac{1}{0.26\sqrt{\pi}} \text{Da}^{1/3} \text{Pe}^{2/3}, \quad (36)$$

Note that these values are obtained with the hypothesis that $g(\tilde{x}) \ll \text{erf}(\tilde{x}\text{Pe}^{1/2})$, which is not fulfilled close to $x=0$. However, by numerically solving Equation 11 (see Figure S5 in Supporting Information S1) we found that this approximation is sufficient to predict the dependency of $s_R$, $R_{\max}$ and $I$ on the Pe and Da (Figure 2).

### 2.2.3. Synthesis of Analytical Results and Comparison With Numerical Simulations

We summarize our analytical derivations in Table 1. These scalings are strictly valid at the stagnation point where the compressions rate is constant. However, as discussed in the following, their relevance extends to a large part of the mixing front where the compression rate is sustained by acceleration along streamline due to fluid incompressibility. While the compression rate fluctuates along streamlines, it is on average equal to the compression rate at the stagnation point in this area. These theoretical results are also applicable to chaotic flows at the pore scale in 3D porous media, which are characterized by exponential elongation of material lines (Heyman et al., 2020; Lester et al., 2013, 2016; Souzy et al., 2020). Due to incompressibility, this leads to similar exponential compression and the development of a Batchelor scale controlling concentration gradients (Equation 4).

To test these scaling laws, we have solved numerically Equation 11 for both the Da < Pe and Da > Pe regimes (see details Figure S3 in Supporting Information S1). Figure 2 shows the results, where the subscript sim represents the results of simulations. Note that we have used dimensional parameters for the numerical simulations. The width was estimated from the standard deviation of the reaction rate profiles across the mixing front (see Equation 41). The results confirm the analytical scaling laws (Figure 2) and the transition between the two regimes at Da ≈ Pe.

## 3. Pore Scale Imaging of Mixing and Reaction at Stagnation Points

### 3.1. Millifluidic Experiments for Mixing and Reaction at Fluid Stagnation Points

To test our analytical predictions, we developed a millifluidic experiment to characterize mixing and reactions in two dimensions in mixing fronts creating by converging flows (Figures 3a and 3b, Figure S1 in Supporting Information S1). We characterized mixing and reaction both locally at the stagnation point and more globally at the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Reactive mixing properties, including reaction width $s_R$, maximum reaction rate $R_{max}$ and reaction intensity $I$, obtained from the numerical solution of Equation 11. (a)–(c) Dependency of reactive mixing properties on Pe at a fixed Da = 1.6. The predicted scaling laws are shown in red and yellow for respectively the low (Equations 32, 35, and 36) and large (Equations 18–20) Pe regimes. (d)–(f) Dependency of reactive mixing properties on Da at a fixed Pe = 1.6. The predicted scaling laws are shown in red and yellow for respectively the low (Equations 18–20) and large (Equations 32, 35, and, 36) Da regimes.</figcaption></figure>

scale of the compressed mixing front. The motivation for studying 2D flows, as in many other studies using 2D micromodels (de Anna, Jimenez-Martinez, et al., 2014; Izumoto et al., 2022; Jiménez-Martínez et al., 2015), is that it allows a detailed imaging of pore scale processes, here mixing-induced reactions. Imaging solute transport at the pore scale in 3D requires precise index matching with solutions of high glycerol content, which currently limits investigations to high Peclet numbers (Heyman et al., 2020; Souzy et al., 2020). It is likely that mixing in converging flows in 3D porous media would show some important differences as in 2D. However, the set up allows understanding important processes that are relevant to 3D, in particular the interplay between exponential fluid compression and pore scale heterogeneity. Exponential compression at a 2D converging flow is locally equivalent to the exponential elongation that occurs globally in 3D porous media (Heyman et al., 2020; Lester et al., 2013, 2016).

Two millifluidic flow cells were designed to study the effect of pore-scale heterogeneity on the reaction kinetics. A classical homogeneous cell with a thickness $h = 0.6$ mm between the parallel plates was first considered to

<table><caption>**Table 1**
Summary of Derived Scaling Laws for the Da > Pe and Da < Pe Regimes</caption><thead><tr><th></th><th>s<sub>R</sub></th><th>R<sub>max</sub></th><th>I</th></tr></thead><tbody><tr><td>Mixing-limited regime Da > Pe</td><td>s<sub>R</sub> ~ Pe<sup>-1/6</sup> Da<sup>-1/3</sup></td><td>R<sub>max</sub> ~ Pe<sup>2/3</sup> Da<sup>1/3</sup></td><td>I ~ Pe<sup>1/2</sup></td></tr><tr><td>Reaction-limited regime Da < Pe</td><td>s<sub>R</sub> ~ Pe<sup>-1/2</sup></td><td>R<sub>max</sub> ~ Da</td><td>I ~ Pe<sup>-1/2</sup> Da</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Illustration of the overall experimental setup. The area surrounded by the black dotted line is the observation zone. (b) Cross-section of the millifluidic setup in the direction normal to the mean flow. The 2D porous media (pillars attached on top of a transparent plate) is placed between the top plate and bottom plate. The water flows in the white space between the pillars. (c) The photo of the 2D porous medium, shown as gray part in (b), and its enlarged photo.</figcaption></figure>

validate the reactive mixing theory presented above. A second cell was designed to study reactive mixing in a two
dimensional porous media composed of cylindrical pillars. The cell with a thickness $h' = 1$ mm was created by
positioning vertical cylindrical pillars, spanning the cell thickness. The ~4,000 pillars were positioned randomly
according to a numerical design ensuring spatially-stationary disorder in the pillars' diameters and distances
between closest neighbors. Its porosity was 0.59. The cell was fabricated using a mold containing a plastic
plate of the required thickness where the pillars' space had been cut out using a laser-cutting machine (Izumoto
et al., 2022). The resulting assembly of pillars confined between two parallel plates constitute a 2D analog porous
medium with pore-scale heterogeneities but Darcy scale homogeneity. The pillar diameters ranged between 0.34
and 0.82 mm with a median size of 0.53 mm. The pore throats (median: 0.36 mm) were smaller than the height
of the cylinders (1 mm). This large aspect ratio was chosen to ensure that the flow field was independent of
the vertical coordinate (Soulaine et al., 2021), except in the vicinity of the top and bottom plates, to minimize
Taylor-Aris dispersion effects. Both cells were made of polydimethylsiloxane (PDMS) (SYLGARD 184 Silicone
Elastomer Kit, Dow, Midland, MI, USA) (Figure 3c). A uniform background flow of one reactant was imposed
in rectangular cuboid cells of length $l' = 150$ mm and width $w = 50$ mm. We injected the second reactant at equal
flow rate by a small needle located in the middle of the cell. This flow configuration creates a flow field similar
to so-called Rankine half-body flow with one stagnation point upstream of the point source (Figure 5).

In the following, we refer to the homogeneous and heterogeneous cells as “Hele-Shaw cell” and “porous
medium,” respectively. We performed numerical simulations of the flow field in the two cells to estimate the
velocity field and the compression rate in the mixing front. The flow field was simulated using a modified
Stokes equation implemented in the open-source software OpenFOAM (https://openfoam.org/), which calculates
the depth-averaged velocity between the two plates (Ferrari et al., 2015; Soulaine et al., 2021). For the porous
medium, we performed three simulations with different locations of the point source injection within a radius of
2 mm from the middle of the cell. Further details of the numerical methods are given in the Figure S2 in Support-
ing Information S1.

**3.2. Chemiluminescence Reaction for Reactive Mixing**

For imaging the reaction field, we used luminol chemiluminescence, where blue light is emitted by the chem-
ical reaction (Burdo & Seitz, 1975; Collaudin & Blum, 1997; Klopf & Nieman, 1983; Uchida et al., 2004). A
previous study applied chemiluminescence for investigating reactive transport in porous media (de Anna, Dentz,
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Reaction rate over time in a fully mixed system (beaker). The red line is the fit of Equation 39 to the data. Note that the fitted lines are different between early times (main plot) and later times (inset).</figcaption></figure>

et al., 2014), but the method required a highly corrosive chemical compound (1,8-diazabicyclo-[5,4,0]-undec-7-ene, corrosive to metal and skin). The luminol chemiluminescence is safer and easier to handle. The luminol chemiluminescence reaction involves a catalytic reaction of H₂O₂ with Co²⁺ followed by an oxidation reaction of luminol with OH· and O₂· radicals. This chain reaction can be written as:

$$ \text{Luminol} + \text{H}_2\text{O}_2 \xrightarrow{\text{Co}^{2+}} \text{3-Aminophthalatedianion} + \text{N}_2 + \text{H}_2\text{O} + \text{photon} \quad (37) $$

Luminol is thus oxidized to 3-Aminophthalatedianion with emission of a blue light photon (λ = 420–460 nm). The reaction rate is thus proportional to the intensity of the emitted blue light. This reaction can be approximated by the bimolecular second-order reaction (Matsumoto & Matsuo, 2015),

$$ A + B \rightarrow C + \text{photon} \quad (38) $$

where A and B are associated to the H₂O₂ and the luminol species respectively. Please note that the actual reaction process is not one-step bimolecular reaction, but it involves several reaction steps. We check that the reaction rate can be approximated by bimolecular reaction in the experiments.

### 3.3. Batch Reaction

We first tested the chemiluminescence reaction (Equation 37) in a batch where the two reactants were fully mixed. For this, two different solutions were prepared. One of them was a mixture of 1 mM Luminol, 7 mM NaOH, 0.01 mM CoCl₂ and 1.95 mM NaCl (called luminol solution in the following), and the other solution (called the H₂O₂ solution in the following) was a mixture of 0.5 mM H₂O₂ and 1.95 mM NaCl. The luminol solution and H₂O₂ solution were mixed and immediately poured in a beaker, and the emitted light was monitored for 3 hr. A rapid decay of the light by a factor 100 was observed over the first 5 s, followed by a much slower decay at long time from 100 s to 3 hr (Figure 4). The solution stayed transparent. We compare these batch kinetics to a second order reaction, for which the rate of reaction R follows (Atkins & De Paula, 2006):

$$ R = \frac{kC_dC_r^2\exp(kC_d t)}{(C_r\exp(kC_d t) - 1)^2} \quad (39) $$

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Examples of concentration fields measured in the conservative tracer experiments with Pe = 0.31. (b) Hele-Shaw cell and (a) porous medium. The concentration is normalized by the maximum concentration c₀. The blue circles indicate the position of the tubing used for the point injections.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Examples of concentration profiles in (a) the Hele-Shaw cell and (b) the porous medium along the line perpendicular to the mixing front at the stagnation points (see Figure 5). The concentration is normalized by the maximum concentration $c_0$. The relative position $x'$ is determined with respect to the point where the concentration is 0.5. Missing data points in porous medium cases are due to the overlapping with grains. The experimental data points (dots) are fitted with an error function (continuous lines, Equation 40) for each Pe. (c) Scaling of the conservative mixing width as a function of Pe for the Hele-Shaw cell (dots) and porous medium (triangles) experiments and comparison with the predicted scaling (black line, Equation 4).</figcaption></figure>

where $k$ is the reaction rate constant to be determined ($s^{-1} \text{ mM}^{-1}$), $C_d = C_{A,0} - C_{B,0} = 0.25 \text{ mM}$, and $C_v = C_{A,0}/C_{B,0} = 2$. We fitted the analytical solution for the bimolecular reaction term $R = kC_A C_B$ (Equation 39) to the data, with the initial solute concentrations of luminol, $C_{A,0} = 0.5 \text{ mM}$ and $H_2O_2$, $C_{B,0} = 0.25 \text{ mM}$. To compare experiments to the prediction Equation 39, the light intensity in the image $I$ was transformed to a reaction rate relative to maximum reaction rate $R = I/a$, where $a$ is the image intensity at time 0. The analytical solution captures the initial decay well with $k_1 \approx 2.56 \text{ s}^{-1} \text{ mM}^{-1}$, indicating that the luminol reaction can be well approximated as a bimolecular reaction at short times. The long-time decay is captured by taking a reaction constant $k_2 = 1.16 \times 10^{-3} \text{ s}^{-1} \text{ mM}^{-1}$, 2,200 times smaller than $k_1$. Although the luminol reaction shows different kinetics at early and late time, it is reasonable to assume that the fast kinetics is dominating in our experiments because the concentration of reactants stays high close to the interface and the light production by the fast kinetics is 100 times larger than by the slow kinetics. We thus take $k = k_1$. The cutoff light intensity below which the image intensity was not sensitive to the reaction rate was also determined for our experimental cell (see Figure S4 in Supporting Information S1).

### 3.4. Experimental Protocol

We first carried out conservative tracer tests in the two types of flow cells. We injected a fluorescein sodium salt solution ($12.5 \text{ mg L}^{-1}$ for the Hele-Shaw cell and $6.25 \text{ mg L}^{-1}$ for the porous medium) from the middle injection port and de-ionized water from the side injection port. Both solutions were injected with the same flow rate for all experiments. The tracer distribution reached a steady state, after which an image was taken. A digital camera (14-bit) with a macro lens was placed above the millifluidic setup. This resulted in an image resolution of 0.017 mm per pixel for the Hele-Shaw cell and 0.016 mm for the porous medium. These resolutions were selected to be below the mixing width (Equation 4) for all flow rates, to ensure that the mixing front can be sufficiently resolved (see Figure 6). A light source with a blue bandpass filter illuminated the millifluidic cell from below to excite the fluorescein sodium salt tracer. A bandpass filter (Narrow Green Bandpass Filter, Midwest Optical Systems, Palatine, IL, USA) was attached to the lens to pass the wave length emitted by the fluorescein. The obtained fluorescent images depict the macroscopic characteristics of the flow well, with a stagnation point upstream of the point source (Figure 5). We varied the Péclet number by increasing the background and point source flow rate simultaneously and step-wisely between 0.14 and $8.88 \text{ mm}^3 \text{ s}^{-1}$. Note that for the same flow rate, the average interstitial velocity has comparable values in the porous medium and in the Hele-Shaw cell, since the porosity of the porous media (0.59) approximately compensates for the smaller thickness of the Hele-Shaw cell (0.6 mm) as compared to that of the porous medium cell (1 mm).

We then carried out reactive transport experiments using the chemiluminescence reaction described above. We injected the luminol solution from the point source injection in the middle of the cell and the $H_2O_2$ solution as the background flow. As in the conservative tracer experiments, both solutions were injected with the same flow rate.
The flow rate was then varied in a similar manner as in the conservative experiments, and the light emission was recorded by the camera. Since the amount of light emitted by the chemiluminescence reaction is generally small, large exposure time and ISO are required. This induced some background noise in the images. To filter such noise, we applied a median filter using the open source image processing software ImageJ (https://imagej.net/software/imagej/) to all reactive images (5 pixels for the porous media and 50 pixels for the Hele-Shaw cell). The median filter selects a pixel with median intensity in a moving window. To obtain the spatial field of a normalized reaction rate, we subtracted the background light intensity from all reactive images, rescaled intensities by the bit depth (2¹⁴ − 1) and multiplied by 0.59 for the porous medium images to compensate for the difference in thickness between the Hele-Shaw cell and the porous medium cell.

The experimental Pe (Equation 5) was estimated with $D = 10^{-9}$ m² s⁻¹ for the molecular diffusion coefficient in water and $L = 0.363$ mm (the median pore throat). The compression rate $\gamma$ was estimated as $\gamma = v_l/l$, where $v_l$ is the mean flow velocity at distance $l$ from the stagnation point. We found that a good estimate of $\gamma$ could be obtained by taking $v_l$ as the average interstitial velocity of the uniform background flow $\bar{v} = Q/(hw)$ ($Q$ being the background flow discharge) and $l$ the distance from the injection point at which the radial flow velocity $v_r = Q/(2\pi hl)$ created by the point injection was equal to $\bar{v}$. This gave $l = w/(2\pi) \approx 8$ mm. This estimate was close to the experimental observation (≈5 mm, see Figure 5a). With these estimates, the Pe varied from $7.63 \times 10^{-2}$ to 4.88 in our experiments. The Da (Equation 13) was estimated from the initial concentration of the luminol in the mixed solution, $C_{A,0} = 0.5$ mM, and $k = 2.56$ s⁻¹ mM⁻¹ determined by batch tests, yielding Da ≈ 169. Hence, Da > Pe for all experiments. We thus explore experimentally the second regime discussed in the theoretical section, where reaction is fast compared to mixing but not instantaneous, leading to fully coupled mixing and reaction (Table 1). We first characterized mixing and reaction locally at the stagnation points and then investigated reactive mixing at the scale of the compressed mixing front.

# 4. Experimental Results

## 4.1. Conservative Mixing

In Figure 5 we show the conservative tracer concentration in the Hele-Shaw cell and the porous medium for Pe = 0.31. Concentration profiles across the mixing front were extracted from the 2D data in sections close to the stagnation point. In the Hele-Shaw cell, this point is easily identified in the middle of the cell (Figure 5a). In the porous medium, the exact location of the stagnation point is more difficult to assess. According to the simulated velocity fields (Figure 11b) and reactive images (Figure 7), the stagnation points are localized on a pillar (Figure 5b) because of the splitting of streamlines upstream of the stagnation points on the same pillar. The mixing interface is thus split between 2 stagnation points on each side of the pillar.

The concentration gradient across the mixing front increases with Pe, in both the Hele-Shaw cell and the porous medium. To estimate the conservative mixing width $s_c$ from the conservative tracer experiments, we fitted an error function to the concentration profiles $C_A(x)$ across the stagnation points:

$$C_A = \frac{1}{2} \left( 1 + \text{erf} \left( \frac{x - \bar{x}}{2s_c} \right) \right), \quad (40)$$

where $\bar{x}$ corresponds to the location where $C_A = 0.5$. Examples of measured concentration profiles and fitted error functions for different flow velocities are shown in Figures 6a and 6b. The mixing widths represented 3–30 pixels per standard deviation, ensuring a sufficient resolution. For the porous medium results, we averaged the fitted mixing widths of the two stagnation points for each Pe. The mixing width in the porous medium was generally smaller and more fluctuating than in the Hele-Shaw cell (Figure 6c). Consistently with the conservative theory (Equation 4), the estimated mixing width scales as predicted by the Batchelor scale $s_B \sim Pe^{-1/2}$ for both the Hele-Shaw cell and porous medium experiments. Note that the image resolution was sufficient to have 3–30 pixels for one standard deviation length. The vertical variation of velocity potentially impacts the measurements since we took images from the top. In that case, the diffusion term in Equation 4 should be replaced by the dispersion due to Taylor-Aris dispersion. This would lead to different scaling exponent over Pe than Equation 4. However, the experimental results were consistent with the prediction, suggesting that the vertical variation of velocity did not have major impact on the measurements. The main reasons is that (a) in the vicinity of the stagnation point velocities go to zero and therefore diffusion dominates over Taylor dispersion (b) downstream in the mixing front mixing occurs mostly transversally and thus Taylor dispersion does not influence it.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Examples of the reaction rate fields for reactive transport experiments (Pe = 3.5) in (a) the Hele-Shaw and (b) the porous medium, and close view of them in the vicinity of the stagnation point in (c) the Hele-Shaw and (d) the porous medium. The reaction rate is normalized by the reaction rate where the image gets saturated ($R_s$), corresponding to the upper limit of the detection range of the camera as explained in the method section. The reaction rate profiles were determined by averaging over the small width of the rectangular domains indicated by the white dotted lines.</figcaption></figure>

## 4.2. Reactive Mixing

We quantify the reaction rate distribution from the light intensity of the chemiluminescence reaction normalized by the upper limit of the detection range of the camera (Figure 7). Namely, the reference reaction rate $R_s$ is defined as the maximum reaction rate recorded by the camera without saturation. As expected, the reaction rate was highest along the mixing interface. In the porous medium, the light intensity was generally larger than in the Hele-Shaw cell, up to four times larger locally, and fluctuations over short length scales were also observed due to pore-scale flow heterogeneities. Figure 8 shows the reaction rate profiles along the line perpendicular to the mixing front at the stagnation points for all Pe. The reaction rate increased with Pe both for the Hele-Shaw cell and the porous medium. From these profiles, we estimated the width $s_R$, maximum reaction rate $R_{max}$ and integrated reaction rate $I$ as a function of Pe (Figure 9). The width was estimated from the standard deviation

$$s_R = \sqrt{\int x^2 P(x) dx - \left(\int x P(x) dx\right)^2}, \quad (41)$$

where $x$ is the position and $P(x) = R(x)/R(x)dx$ with $R(x)$ the normalized reaction rate at position $x$. The estimated widths ranged between 12 and 36 pixels per standard deviation, ensuring a sufficient resolution. For the porous medium, $s_R$, $R_{max}$, and $I$ were averaged over the two sides of the cylinder. In the Hele-Shaw cell at Pe < 0.61, the reactive zone began to overlap with the injection region, hence small flow-rates were not considered.

The scaling of $s_R$, $R_{max}$ and $I$ with Pe were similar for the Hele-Shaw cell and the porous medium (Figure 9). The reaction rate and reaction intensity globally increased with Pe, while the width of the reaction front tends to decrease with Pe. The reactive front in the porous medium yielded larger $R_{max}$ and $I$ than those in Hele-Shaw cell, while the reactive scale $s_R$ was smaller in the porous medium. In addition to the local reaction rate at the stagnation point, we estimated the total reaction rate over a larger area under compression (area shown in Figures 7a and 7b). We found that the total reaction rate obeyed the same Pe$^{1/2}$ scaling as observed close to the stagnation point for both the porous medium and the Hele-Shaw cell, with a total reaction rate approximately two times larger in the porous medium than in the Hele-Shaw cell.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Reaction rate profiles normalized by the reaction rate where the image gets saturated ($R_s$) across the stagnation points for different Pe in (a) the Hele-Shaw cell and (b) the porous medium stagnation point #1 and (c) porous medium stagnation point #2. The relative position $x'$ is determined with respect to the point where the reaction rate takes its maximum value at the highest Pe.</figcaption></figure>

# 5. Discussion

## 5.1. Enhanced Mixing and Reaction in Converging Flows

For the Hele-Shaw cell and the porous medium, the experimental results (Figure 9) were in relatively good agreement with the theoretical scaling laws derived for Da > Pe regime for the maximum reaction rate (Equation 36) and the reaction intensity (Equation 32). The evolution of the reaction intensity with Pe was consistent with the prediction (Equation 32) both at the local stagnation point scale (Figure 9c) and at the scale of the broader section of the mixing front that was under compression (Figure 9d). For the reaction front width, the scaling exponent of the experimental results was slightly larger than the expected scaling exponent (Figure 9a). Our hypothesis is that the width estimate may be sensitive to deviations from the pure bimolecular reaction at low concentrations (Figure 4 inset). In addition, this reflects the difficulty to precisely evaluate such small changes of the reaction front width. Overall, these findings show that both the conservative and reactive mixing dynamics can be well understood and quantified from the mechanisms of coupled compression, diffusion and reaction (Table 1). As

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Evolution of the reaction front properties as a function of Pe: (a) width of the reaction front, (b) maximum reaction rate normalized by the reaction rate where the image gets saturated ($R_s$), (c) reaction intensity and (d) total reaction rate integrated over the area shown in Figures 7a and 7b. The predicted scaling laws are shown as black lines for the reaction width (Equation 35), the maximum reaction rate (Equation 36) and the reaction intensity (Equation 32).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Calculated ratio of the compression rate $\alpha$ (Porous medium/Hele-Shaw cell) from experimental and simulation results. The conservative and reactive cases were calculated from the experimental results, with the reactive case calculated by $(\alpha_{R_{max}} + \alpha_{I_R})/2$. Results for “kinematics” are calculated from the velocity gradient in the numerical simulations illustrated in Figure 11. The solid line is the average $\alpha$ from the reactive experiments.</figcaption></figure>

in many reactive transport systems, incomplete mixing at the pore scale (de Anna, Jimenez-Martinez, et al., 2014) determines here macroscopic reaction rates. We have thus shown that the scaling of the reaction rate with Péclet at the scale of the mixing front (Figure 9d) follows the same law as the pore scale reaction rate at the stagnation point (Figure 9c).

The compression rate $\gamma$ induced by the converging flows is equal to the mean velocity gradient created by the converging flows $\gamma = \bar{v}/l$, where $l$ is the characteristic distance over which the velocities decay from the background flows to the mixing front. In 3D porous media, the exponential compression due to converging flows in the mixing front should combine with exponential elongation to chaotic pore scale flow (Lester et al., 2016). The compression rate induced by chaotic mixing is on the order of $\gamma = 0.2\bar{v}/d$, with $d$ the grain size (Heyman et al., 2020). Therefore, in 3D flows, compression induced by chaotic flows would dominate over that of converging flows for characteristic distances of velocity changes toward the front much larger than five grain sizes $l \gg d/0.2$. In all cases, mixing and reaction would obey the predictions presented here for exponential compression.

## 5.2. Enhanced Compression Rates in Porous Media

The experimental results (Figures 6c and 9a), suggest that the mixing and reaction widths in the porous media are smaller than in the Hele-Shaw cell.

This leads to larger concentration gradients and therefore larger mixing and reaction rates (Figures 7 and 9). This suggests that pore-scale heterogeneities enhance local fluid compression. To quantify the enhancement of compression in porous media, we define the enhancement ratio by $\alpha = \gamma_p/\gamma_H$, where $\gamma_p$ is the compression rate in the porous medium and $\gamma_H$ is the compression rate in the Hele-Shaw cell. Introducing $\alpha$ in Equations 5, 32, and 36, we obtain the following estimates of $\alpha$ from the ratio of variables measured in the porous medium (denoted by subscript $p$) and the Hele-Shaw cell (denoted by subscript $H$):

$$ \alpha_{s_c} = \left( \frac{S_{c,H}}{S_{c,P}} \right)^2, \qquad (42) $$

for the conservative mixing width,

$$ \alpha_{R_{max}} = \left( \frac{R_{max,P}}{R_{max,H}} \right)^{3/2} \qquad (43) $$

for the maximum reaction rate, and

$$ \alpha_{I_R} = \left( \frac{I_{R,P}}{I_{R,H}} \right)^2 \qquad (44) $$

for the integral of the reaction rate.

We plot in Figure 10 these estimates of $\alpha$ at the stagnation points as a function of Péclet number. In the conservative experiments, $\alpha_{s_c}$ ranges between 1 and 3. For the reactive experiments, we plot $(\alpha_{R_{max}} + \alpha_{I_R})/2$ for each Pe, which ranges between 3 and 8, with a mean value of 4.4 (see Text S6 in Supporting Information S1 for $\alpha_{R_{max}}$ and $\alpha_{I_R}$). These large fluctuations are likely due to a change in the location of the local stagnation points, leading to different local mixing rates. These experimental estimates are compared to the local compression rates obtained in numerical simulations (Figure 11, see details of computation in S3). From the simulated flow fields, we find a ratio $\alpha_N = 4.7$, close to the value obtained from the reactive experiments. The ratio of compression rates estimated from conservative experiments $\alpha_{s_c}$ is smaller than those estimated from reactive experiments and from numerical simulations. Our hypothesis is that this might be attributed to an overestimation of the mixing width in conservative experiments due to light scattering.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) Experimental results of reaction rate fields in the Hele-Shaw cell (top) and the porous medium (bottom). The reaction rate is normalized by the reaction rate where the image gets saturated ($R_c$). The yellow dotted circles represent the position of the stagnation points. (b) Simulated streamlines and velocity fields in the Hele-Shaw cell (top) and the porous medium (bottom) near the stagnation points for $Pe = 1.2$. The inflow is from the left and right boundaries, and the outflow is through the top and bottom boundaries. The white arrows indicate the flow directions. (c) Calculated compression rate along the x axis around stagnation points for the Hele-Shaw cell (top) and the porous medium (bottom) for $Pe = 1.2$. Black dotted lines indicate the positions to measure the local compression rate over the mixing width (results shown in Figure 10). The length of the line was 0.6 mm, which was $2s_R$ for the intermediate $Pe$ in the Hele-Shaw cell.</figcaption></figure>

### 5.3. Origin of Enhanced Local Compression in Porous Media

The larger reaction rate over the mixing front in the porous medium suggests that the local compression rates are about four times larger than in the Hele-Shaw cell although the mean interstitial velocities are identical. In the following, we propose a heuristic argument to explain this enhancement. We argue that stretching upstream and compression downstream of grains do not balance when a fluid element is accelerating as it travels in the flow field (i.e., in Lagrangian framework), leading to a net compression as the fluid travels around each grain. Along streamlines, the acceleration in the porous media fluctuates around an average positive acceleration (Figure 12a). The development of reactive tails downstream of grains in the mixing zone suggests that there is enhanced compression behind grains (Figure 12b). Let us consider such a grain in the compressed mixing zone where flow is accelerating. As the interface flows around the grain, its width is stretched upstream of the grain and compressed downstream. At a distance of about $d/2$ upstream of the grain, the velocity is on average similar to the one in the Hele Shaw cell at the same location $s$ on the mixing interface (Figure 12b). As it approaches the grain, the interface experiences a negative acceleration $\gamma^{up} \approx (0 - v_H(s))/(d/2)$ (Figure 12c). Owing to incompressibility, this leads to a stretching rate of $-\gamma^{up}$ in the direction perpendicular to the streamlines. Similarly, downstream of the grain, flow accelerates at a rate $\gamma^{down} \approx (v_H(s + 2d) - 0)/(d/2)$, leading to a compression rate of $\gamma^{down}$ in the direction perpendicular to the streamlines. The resulting net compression rate is $\gamma^{net} = \gamma^{down} + \gamma^{up} = (v_H(s + 2d) - v_H(s))/(d/2) = 4\gamma_H$ (Figure 12c). Locally at the fluid stagnation point, the same reasoning applies when taking
<figure><img src="image_14.png" /><figcaption>**Figure 12.** (a) Simulated acceleration along the streamlines from the stagnation points to downstream for Pe = 2.4, for the porous medium (red) and Hele-Shaw cell (blue). (b) Image of reaction rate showing the reactive tails behind grains. The white square focuses on a grain in the mixing zone with such enhanced mixing and reaction, motivating the sketch below. (c) Sketch of compression enhancement over a grain in an accelerating flow leading to reactive tails behind grains.</figcaption></figure>

$s = 0$ at the Hele-Shaw stagnation point and considering the downstream compression: $\gamma^{down} = (v_H(2d) - 0)/(d/2) = 4\gamma_H$. Note that this mechanism involves discontinuous streamlines which respectively end and start around the grains. Therefore, it is not captured when measuring acceleration along a specific continuous streamline in the mixing front (Figure 12a).

This mechanism may thus explain the enhancement of compression by a factor of about 4 between homogeneous and heterogeneous flows, and the consequences on reaction rates both locally at the stagnation and over the compression zone. Note that the proposed mechanism implies that mixing enhancement is independent of the grain diameter and would not occur in a regular coflow (two solutions flow in parallel), with no acceleration. This hypothesis remains to be tested using systems of different grain sizes and heterogeneity, and using other flow topologies, such as simple coflow systems. This would also allow investigating the effect of shear on this mixing enhancement.

## 6. Conclusions

We investigated the dynamics of mixing and reaction in mixing fronts developing at the interface of converging flows, both in homogeneous Hele-Shaw cells and porous media. Using the lamellar theory of mixing, we derived analytical solutions quantifying the dependence of mixing and reaction rates upon the Péclet and Damköhler numbers in such compressed mixing fronts. To test these theoretical predictions, we developed conservative and reactive transport millifluidic experiments in a homogeneous Hele-Shaw. A mixing front was created at the interface of opposing flows and pore scale mixing and reaction rates were quantified at the pore scale both at the stagnation point and through the mixing front. The evolution of the mixing and reaction properties with Péclet number were in good agreement with the theory both locally at the scale of the front. This demonstrates that the interplay between fluid compression, molecular diffusion and chemical reactions leads to enhanced mixing and reaction. Furthermore, we observed that mixing and

reaction rates in the porous medium were approximately four times larger than in the Hele-Shaw cell. To explain this enhancement, we proposed a simple mechanism based on the imbalance between stretching and compression upstream and downstream of grains, in accelerating flows. Future work should focus on validating this mechanism by varying grain sizes and distributions.

These findings provide new insights into the dynamics of mixing and reactive fronts, where fluids of different chemical compositions converge and react. Converging flows are regions of exponential fluid compression and enhanced concentration gradients commonly found in a large range of subsurface flows, where they act as mixing and reaction hot spots. These results are also relevant for mixing and reaction in chaotic flows, which produce similar exponential compression of mixing interfaces. The improvement of experimental techniques to visualize reaction rate field within 3D porous media should allow investigating the dynamics of mixing and reaction in 3D converging flows in the future. At larger scales, when the width of the mixing front is much larger than pore sizes, dispersion should play significant role, which should be integrated into the theory.

## Acknowledgments

Funding for this study was provided by European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No 722028 to the project ENIGMA (European training Network for In situ imaGing of dynaMic processes in heterogeneous subsurfAce environments) and by ANR JCJC SUCHY ANR-19-CE01-0013, as well as the CPER (Contrat de Plan État-Région) BUFFON for part of the equipment. We acknowledge funding by ERC project ReactiveFront 648377.

## Data Availability Statement

The data files used in this paper are available at Izumoto (2023).

## References

Atkins, P. W., & De Paula, J. (2006). *Physical chemistry for the life sciences*. W. H. Freeman and Company.

Bandopadhyay, A., Davy, P., & Le Borgne, T. (2018). Shear flows accelerate mixing dynamics in hyporheic zones and hillslopes. *Geophysical Research Letters*, 45(21), 11–659. https://doi.org/10.1029/2018gl079914
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. *Advances in Water Resources, 100*, 78–95. https://doi.org/10.1016/j.advwatres.2016.12.008

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. *Journal of Fluid Mechanics*, 5(01), 134. https://doi.org/10.1017/S0022112059000106

Benson, D. A., Aquino, T., Bolster, D., Engdahl, N., Henri, C. V., & Fernandez-Garcia, D. (2017). A comparison of Eulerian and Lagrangian transport and non-linear reaction algorithms. *Advances in Water Resources*, 99, 15–37. https://doi.org/10.1016/j.advwatres.2016.11.003

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

Bolster, D., Valdés-Parada, F. J., LeBorgne, T., Dentz, M., & Carrera, J. (2011). Mixing in confined stratified aquifers. *Journal of Contaminant Hydrology*, 120, 198–212. https://doi.org/10.1016/j.jconhyd.2010.02.003

Bresciani, E., Kang, P. K., & Lee, S. (2019). Theoretical analysis of groundwater flow patterns near stagnation points. *Water Resources Research*, 55(2), 1624–1650. https://doi.org/10.1029/2018WR023508

Briggs, M. A., Lautz, L. K., & Hare, D. K. (2014). Residence time control on hot moments of net nitrate production and uptake in the hyporheic zone. *Hydrological Processes*, 28(11), 3741–3751. https://doi.org/10.1002/hyp.9921

Burdo, T. G., & Seitz, W. R. (1975). Mechanism of cobalt catalysis of luminol chemiluminescence. *Analytical Chemistry*, 47(9), 1639–1643. https://doi.org/10.1021/ac60359a019

Burté, L., Cravotta, C. A., III, Bethencourt, L., Farasin, J., Pédrot, M., Dufresne, A., et al. (2019). Kinetic study on clogging of a geothermal pumping well triggered by mixing-induced biogeochemical reactions. *Environmental Science & Technology*, 53(10), 5848–5857. https://doi.org/10.1021/acs.est.9b00453

Collaudin, A. B., & Blum, L. J. (1997). Investigations of the enhancer effect of a high-salt concentration medium on the luminol chemiluminescent reaction. *Photochemistry and Photobiology*, 65(2), 303–308. https://doi.org/10.1111/j.1751-1097.1997.tb08562.x

de Anna, P., Dentz, M., Tartakovskiy, A., & Le Borgne, T. (2014). The filamentary structure of mixing fronts and its control on reaction kinetics in porous media flows. *Geophysical Research Letters*, 41(May), 4586–4593. https://doi.org/10.1002/2014GL060068

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. *Environmental Science & Technology*, 48(1), 508–516. https://doi.org/10.1021/es403105b

de Anna, P., Le Borgne, T., Dentz, M., Tartakovskiy, A. M., Bolster, D., & Davy, P. (2013). Flow intermittency, dispersion, and correlated continuous time random walks in porous media. *Physical Review Letters*, 110(18), 184502. https://doi.org/10.1103/PhysRevLett.110.184502

De Barros, F. P., Dentz, M., Koch, J., & Nowak, W. (2012). Flow topology and scalar mixing in spatially heterogeneous flow fields. *Geophysical Research Letters*, 39(8), 1–5. https://doi.org/10.1029/2012gl051302

Dentz, M., Hidalgo, J. J., & Lester, D. (2023). Mixing in porous media: Concepts and approaches across scales. *Transport in Porous Media*, 146(1–2), 5–53. https://doi.org/10.1007/s11242-022-01852-x

Dentz, M., Le Borgne, T., Englert, A., & Bijeljic, B. (2011). Mixing, spreading and reaction in heterogeneous media: A brief review. *Journal of Contaminant Hydrology*, 120–121(C), 1–17. https://doi.org/10.1016/j.jconhyd.2010.05.002

De Vriendt, K., Le Borgne, T., Pool, M., & Dentz, M. (2022). Subsurface mixing dynamics across the salt-freshwater interface. *Geophysical Research Letters*, 49(7), 1–10. https://doi.org/10.1029/2021GL0970

Engdahl, N. B., Benson, D. A., & Bolster, D. (2014). Predicting the enhancement of mixing-driven reactions in nonuniform flows using measures of flow topology. *Physical Review E*, 90(5), 051001. https://doi.org/10.1103/physreve.90.051001

Ferrari, A., Jimenez-Martinez, J., Le Borgne, T., Méheust, Y., & Lunati, I. (2015). Challenges in modeling unstable two-phase flow experiments in porous micromodels. *Water Resources Research*, 51(3), 1381–1400. https://doi.org/10.1002/2014wr016384

Fu, F., Dionysiou, D. D., & Liu, H. (2014). The use of zero-valent iron for groundwater remediation and wastewater treatment: A review. *Journal of Hazardous Materials*, 267, 194–205. https://doi.org/10.1016/j.jhazmat.2013.12.062

Gramling, C. M., Harvey, C. F., & Meigs, L. C. (2002). Reactive transport in porous media: A comparison of model prediction with laboratory visualization. *Environmental Science & Technology*, 36(11), 2508–2514. https://doi.org/10.1021/es0157144

Gu, C., Anderson, W., & Maggi, F. (2012). Riparian biogeochemical hot moments induced by stream fluctuations. *Water Resources Research*, 48(9), 1–17. https://doi.org/10.1029/2011WR011720

Hester, E. T., Cardenas, M. B., Haggerty, R., & Apte, S. V. (2017). The importance and challenge of hyporheic mixing. *Water Resources Research*, 53(5), 3565–3575. https://doi.org/10.1002/2016WR020005

Heyman, J., Lester, D. R., Turuban, R., Méheust, Y., & Le Borgne, T. (2020). Stretching and folding sustain microscale chemical gradients in porous media. *Proceedings of the National Academy of Sciences of the United States of America*, 117(24), 13359–13365. https://doi.org/10.1073/pnas.2002858117

Hidalgo, J. J., Dentz, M., Cabeza, Y., & Carrera, J. (2015). Dissolution patterns and mixing dynamics in unstable reactive flow. *Geophysical Research Letters*, 42(15), 6357–6364. https://doi.org/10.1002/2015GL065036

Hochstetler, D. L., & Kitanidis, P. K. (2013). The behavior of effective rate constants for bimolecular reactions in an asymptotic transport regime. *Journal of Contaminant Hydrology*, 144(1), 88–98. https://doi.org/10.1016/j.jconhyd.2012.10.002

Izumoto, S. (2023). [data]: Enhanced mixing and reaction at fluid stagnation points: Theory and pore scale imaging. https://doi.org/10.5281/zenodo.7950543

Izumoto, S., Huisman, J. A., Zimmermann, E., Heyman, J., Gomez, F., Tabuteau, H., et al. (2022). Pore-scale mechanisms for spectral induced polarization of calcite precipitation inferred from geo-electrical millifluidics. *Environmental Science & Technology*, 56(8), 4998–5008. https://doi.org/10.1021/acs.est.1c07742

Jiménez-Martínez, J., Anna, P. D., Tabuteau, H., Turuban, R., Borgne, T. L., & Méheust, Y. (2015). Pore-scale mechanisms for the enhancement of mixing in unsaturated porous media and implications for chemical reactions. *Geophysical Research Letters*, 42(13), 5316–5324. https://doi.org/10.1002/2015GL064513

Kitanidis, P. K., & McCarty, P. L. (2012). *Delivery and mixing in the subsurface: Processes and design principles for in situ remediation* (Vol. 4). Springer Science & Business Media.

Klopf, L. J., & Nieman, T. A. (1983). Effect of iron(II), cobalt(II), copper(II), and manganese(II) on the chemiluminescence of luminol in the absence of hydrogen peroxide. *Analytical Chemistry*, 55(7), 1080–1083. https://doi.org/10.1021/ac00258a023

Larralde, H., Araujo, M., Havlin, S., & Stanley, H. E. (1992). Reaction front for A + B → C diffusion-reaction systems with initially separated reactants. *Physical Review A*, 46(2), 855–859. https://doi.org/10.1103/PhysRevA.46.855

Le Borgne, T., Dentz, M., & Villermaux, E. (2015). The lamellar description of mixing in porous media. *Journal of Fluid Mechanics*, 770, 458–498. https://doi.org/10.1017/jfm.2015.117
Le Borgne, T., Ginn, T. R., & Dentz, M. (2014). Impact of fluid deformation on mixing-induced chemical reactions in heterogeneous flows. *Geophysical Research Letters*, *41*(22), 7898-7906. https://doi.org/10.1002/2014GL062038

Lee, S. H., & Kang, P. K. (2020). Three-dimensional vortex-induced reaction hot spots at flow intersections. *Physical Review Letters*, *124*(14), 144501. https://doi.org/10.1103/physrevlett.124.144501

Lester, D. R., Metcalfe, G., & Trefry, M. G. (2013). Is chaotic advection inherent to porous media flow? *Physical Review Letters*, *111*(17), 174101. https://doi.org/10.1103/PhysRevLett.111.174101

Lester, D. R., Metcalfe, G., & Trefry, M. G. (2014). Anomalous transport and chaotic advection in homogeneous porous media. *Physical Review E*, *90*(6), 063012. https://doi.org/10.1103/PhysRevE.90.063012

Lester, D. R., Trefry, M. G., & Metcalfe, G. (2016). Chaotic advection at the pore scale: Mechanisms, upscaling and implications for macroscopic transport. *Advances in Water Resources*, *97*, 175-192. https://doi.org/10.1016/j.advwatres.2016.09.007

Lysy, M., Ersland, G., & Fernø, M. (2022). Pore-scale dynamics for underground porous media hydrogen storage. *Advances in Water Resources*, *163*, 104167. https://doi.org/10.1016/j.advwatres.2022.104167

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. *Proceedings of the National Academy of Sciences*, *114*(17), 4330–4335. https://doi.org/10.1073/pnas.1617454114

Matsumoto, R., & Matsuo, R. (2015). Measurement of local chemical reaction rate in a microchannel by using luminol chemiluminescence. *Experimental Thermal and Fluid Science*, *67*, 88–95. https://doi.org/10.1016/j.expthermflusci.2014.11.013

Mays, D. C., & Neupauer, R. M. (2012). Plume spreading in groundwater by stretching and folding. *Water Resources Research*, *48*(7), W07501. https://doi.org/10.1029/2011wr011567

McClain, M. E., Boyer, E. W., Dent, C. L., Gergel, S. E., Grimm, N. B., Groffman, P. M., et al. (2003). Biogeochemical hot spots and hot moments at the interface of terrestrial and aquatic ecosystems. *Ecosystems*, *6*(4), 301–312. https://doi.org/10.1007/s10021-003-0161-9

Piscopo, A. N., Neupauer, R. M., & Mays, D. C. (2013). Engineered injection and extraction to enhance reaction for improved in situ remediation. *Water Resources Research*, *49*(6), 3618–3625. https://doi.org/10.1002/wrcr.20209

Ranz, W. E. (1979). Applications of a stretch model to mixing, diffusion, and reaction in laminar and turbulent flows. *AIChE Journal*, *25*(1), 41–47. https://doi.org/10.1002/aic.690250105

Rolle, M., & Le Borgne, T. (2019). Mixing and reactive fronts in the subsurface. *Reviews in Mineralogy and Geochemistry*, *85*(1), 111–142. https://doi.org/10.2138/rmg.2018.85.5

Sole-Mari, G., Fernández-Garcia, D., Sanchez-Vila, X., & Bolster, D. (2020). Lagrangian modeling of mixing-limited reactive transport in porous media: Multirate interaction by exchange with the mean. *Water Resources Research*, *56*(8), e2019WR026993. https://doi.org/10.1029/2019wr026993

Soler-Sagarra, J., Carrera, J., Bonet, E., Roig, C., & Becker, P. (2023). Modeling mixing in stratified heterogeneous media: The role of water velocity discretization in phase space formulation. *Transport in Porous Media*, *146*(1–2), 395–412. https://doi.org/10.1007/s11242-022-01795-3

Soulaine, C., Maes, J., & Roman, S. (2021). Computational microfluidics for geosciences. *Frontiers in Water*, *3*, 1–21. https://doi.org/10.3389/frwa.2021.643714

Souzy, M., Lhuissier, H., Méheust, Y., Le Borgne, T., & Metzger, B. (2020). Velocity distributions, dispersion and stretching in three-dimensional porous media. *Journal of Fluid Mechanics*, *891*, A16. https://doi.org/10.1017/jfm.2020.113

Szulczewski, M. L., MacMinn, C. W., Herzog, H. J., & Juanes, R. (2012). Lifetime of carbon capture and storage as a climate-change mitigation technology. *Proceedings of the National Academy of Sciences*, *109*(14), 5185–5189. https://doi.org/10.1073/pnas.1115347109

Tarkowski, R. (2019). Underground hydrogen storage: Characteristics and prospects. *Renewable and Sustainable Energy Reviews*, *105*, 86–94. https://doi.org/10.1016/j.rser.2019.01.051

Tennekes, H., & Lumley, J. L. (1972). *A first course in turbulence*. The MIT Press. https://doi.org/10.7551/mitpress/3014.001.0001

Turuban, R., Lester, D. R., Heyman, J., Borgne, T. L., & Méheust, Y. (2019). Chaotic mixing in crystalline granular media. *Journal of Fluid Mechanics*, *871*, 562–594. https://doi.org/10.1017/jfm.2019.245

Turuban, R., Lester, D. R., Le Borgne, T., & Méheust, Y. (2018). Space-group symmetries generate chaotic fluid advection in crystalline granular media. *Physical Review Letters*, *120*(2), 024501. https://doi.org/10.1103/physrevlett.120.024501

Uchida, S., Satoh, Y., Yamashiro, N., & Satoh, T. (2004). Determination of hydrogen peroxide in water by chemiluminescence detection, (II). *Journal of Nuclear Science and Technology*, *41*(9), 898–906. https://doi.org/10.1080/18811248.2004.9715562

Valocchi, A. J., Bolster, D., & Werth, C. J. (2019). Mixing-limited reactions in porous media. *Transport in Porous Media*, *130*(1), 157–182. https://doi.org/10.1007/s11242-018-1204-1

Villermaux, E. (2019). Mixing versus stirring. *Annual Review of Fluid Mechanics*, *51*(1), 245–273. https://doi.org/10.1146/annurev-fluid-010518-040306

Wang, Y., Fernández-Garcia, D., Sole-Mari, G., & Rodríguez-Escales, P. (2022). Enhanced NAPL removal and mixing with engineered injection and extraction. *Water Resources Research*, *58*(4), e2021WR031114. https://doi.org/10.1029/2021wr031114

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. *Environmental Science & Technology*, *42*(9), 3185–3193. https://doi.org/10.1021/es7022835

Ye, Y., Chiogna, G., Cirpka, O. A., Grathwohl, P., & Rolle, M. (2015). Experimental evidence of helical flow in porous media. *Physical Review Letters*, *115*(19), 194502. https://doi.org/10.1103/physrevlett.115.194502

Ziliotto, F., Basilio Hazas, M., Rolle, M., & Chiogna, G. (2021). Mixing enhancement mechanisms in aquifers affected by hydpeaking: Insights from flow-through laboratory experiments. *Geophysical Research Letters*, *48*(21), e2021GL095336. https://doi.org/10.1029/2021gl095336

Zoback, M. D., & Gorelick, S. M. (2012). Earthquake triggering and large-scale geologic storage of carbon dioxide. *Proceedings of the National Academy of Sciences*, *109*(26), 10164–10168. https://doi.org/10.1073/pnas.1202473109

## References From the Supporting Information

De, N., Meunier, P., Méheust, Y., & Nadal, F. (2021). Bi-dimensional plume generated by the convective dissolution of an extended source of CO₂. *Physical Review Fluids*, *6*(6), 063503. https://doi.org/10.1103/physrevfluids.6.063503

Le Borgne, T., Dentz, M., & Villermaux, E. (2013). Stretching, coalescence, and mixing in porous media. *Physical Review Letters*, *110*(20), 1–5. https://doi.org/10.1103/PhysRevLett.110.204501