Effective reaction kinetics of steady mixing fronts in porous media

Abstract

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.

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; 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\to 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\to 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 $\boldsymbol{v}$ , the dispersion tensor $\unicode{x1D63F}$ , 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:

(2.1) \begin{align} \frac {\partial {C_A}}{\partial {t}} + {\boldsymbol{{v}}}\cdot \boldsymbol{\nabla} {C_A} &= \boldsymbol{\nabla} \cdot \left (\unicode{x1D63F}\,\boldsymbol{\nabla} {C_A}\right )-{k}{C_A}{C_B}, \notag \\ \frac {\partial {C_B}}{\partial {t}} + {\boldsymbol{{v}}} \cdot \boldsymbol{\nabla} {C_B} &= \boldsymbol{\nabla} \cdot \left (\unicode{x1D63F}\,\boldsymbol{\nabla} {C_B}\right )-{k}{C_A}{C_B}. \end{align}

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)

(2.2) \begin{equation} \unicode{x1D63F} = \left (D_{{m}} + \alpha _T\, |\boldsymbol{v}| \right )\, \unicode{x1D644}+ \left (\alpha _L-\alpha _T\right ) {\left (\boldsymbol{v} \otimes \boldsymbol{v}\right )}/{|\boldsymbol{v}|}, \end{equation}

where $\alpha _L$ and $\alpha _T$ are the longitudinal and transverse dispersivities, respectively, $D_{{m}}$ is the molecular diffusivity, and $\unicode{x1D644}$ 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,

(2.3) \begin{equation} Pe = \frac {U L}{D_{{m}}} \end{equation}

and

(2.4) \begin{equation} Da= \frac {k C_A^0 L^2}{D_{{m}}}. \end{equation}

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\gt 0$ , and $(C_A,C_B)=(0,q C_A^0)$ for $y\lt 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 = k C_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

(2.5) \begin{equation} R_{{max}}(x) = \max _y {R(x,y)}, \end{equation}

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

(2.6) \begin{equation} I(x) = \int \text{d}y\, {R(x,y)}, \end{equation}

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

(2.7) \begin{equation} s^2_R= \frac {\int \text{d}x\, y^2\, R(x,y)}{\int \text{d}x\, R(x,y)}. \end{equation}

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 $\boldsymbol{v} \equiv (U,0)$ (figure 1) and mixing transversely through diffusion only (low $Pe$ ), such that (2.2) reduces to $\unicode{x1D63F} \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

(2.8) \begin{align} \frac {\partial C_A}{\partial t} &= D_{{m}}\frac {\partial ^2 C_A}{\partial y^2}- k C_A C_B,\notag \\ \frac {\partial C_B}{\partial t} &= D_{{m}}\frac {\partial ^2 C_B}{\partial y^2}- k C_A C_B. \end{align}

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\gt 0$ . The colour scale indicates local reaction rates obtained in numerical simulations of the hydrodynamic dispersion approximation (see § A.4).

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,

(2.9) \begin{equation} 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}. \end{equation}

In contrast, at times larger than the reaction time ( $t\gt t_r$ , $x\gt 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

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

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)$ :

(2.11) \begin{equation} \boldsymbol{v} = \gamma \left (\begin{array}{c} x \\ - y \end{array} \right ), \end{equation}

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 ( $\boldsymbol{y}$ ) and one unstable manifold ( $\boldsymbol{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 :

(2.12) \begin{equation} \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. \end{equation}

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

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

(2.13) \begin{equation} 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 k L\, Pe^{-1/2}, \end{equation}

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

(2.14) \begin{align} s_R &\sim L\, Da^{-1/3} \, Pe^{-1/6}, \nonumber \\ R_{{max}} &\sim \frac { C_A^0 D_{{m}}}{L^2}\, Da^{1/3} \, Pe^{2/3},\nonumber \\ I &\sim \frac { C_A^0 D_{{m}}}{L}\, Pe^{1/2}. \end{align}

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

(2.15) \begin{equation} D_T=\alpha _T U = D_m\, Pe\, \frac {\alpha _T}{d}, \end{equation}

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

(2.16) \begin{equation} Pe_D=\frac {L U}{D_T}=\frac {d}{\alpha _T}, \end{equation}

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

(2.17) \begin{equation} Da_D= \frac {Da}{Pe} \frac {d}{\alpha _T}. \end{equation}

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

(2.18) \begin{equation} s_r \sim {\left ( {\alpha _T x}\right )}^{1/2}, \quad R_{{max}} \sim (C_A^0)^2 k \quad \text{and}\quad I \sim (C_A^0)^2 k \, {\left ( {\alpha _T x} \right )}^{1/2}, \end{equation}

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

(2.19) \begin{align} 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}, \nonumber \\ \nonumber 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}. \end{align}

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 $\boldsymbol{v}(x,0)=(\gamma x,0)$ . Thus the transverse dispersivity is spatially dependent:

(2.20) \begin{equation} D_T=\alpha _T \gamma x = D_{{m}}\, Pe\, x \alpha _T/d^2 . \end{equation}

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

(2.21) \begin{equation} s_r \sim {\left ( {\alpha _T x} \right )}^{1/2}, \quad R_{{max}} \sim (C_A^0)^2 k\quad \text{and} \quad I \sim (C_A^0)^2 k {\left ( {\alpha _T x} \right )}^{1/2}. \end{equation}

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

(2.22) \begin{align} s_r &\sim Da^{-1/3} Pe^{1/3} {\left ( {\alpha _T x} \right )}^{1/2}, \nonumber \\ \nonumber R_{{max}} &\sim \frac {C_A^0 D_{{m}}}{d^2} \, Da^{1/3} Pe^{2/3} , \\ I &\sim \frac {C_A^0 D_{{m}}}{d^2} \, Pe\, \left ({\alpha _T x }\right )^{1/2}, \end{align}

where $Pe=\gamma d^2/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.

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.

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 $\mathrm{H_2 O_2}$ with $\mathrm{Co^{2+}}$ followed by an oxidation reaction of the luminol compound with $\text{OH}\cdot$ and $\mathrm{O_2}\cdot$ radicals (Uchida et al. 2004). These reactions can be written as

(3.1) \begin{equation} \mathrm{luminol + H_2 O_2 + 2 OH- \xrightarrow{\,\,\,Co^2+\,\,\,} \text{3-aminophthalatedianion} + N_2 + H_2 O + {photon}} \end{equation}

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:

(3.2) \begin{equation} \mathrm{A + B \longrightarrow C + {photon}}, \end{equation}

where $A$ and $B$ represent the $\mathrm{H_2 O_2}$ 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 $\mathrm{CoCl_2}$ (termed luminol solution), and the other was a mixture of 0.5 mM $\mathrm{H_2 O_2}$ and 3.9 mM NaCl (termed $\mathrm{H_2 O_2}$ solution). Izumoto et al. (2023) have estimated the reaction constant $k=0.08$ $\mathrm{s^{-1}\ 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/(k C_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 $\mathrm{H_2 O_2}$ solution may contain stabiliser compounds leading to a departure from the expected bimolecular kinetic. This did not occur for our 30 % $\mathrm{H_2 O_2}$ 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 3 a), 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 3 b), 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$ 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 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, $\gamma ^{-1}\gt t_r$ . Hence they are in the mixing-limited regime. A summary of the experimental runs and associated parameters is provided in table 2.

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 $\mathrm{mm}^2$ for the Hele-Shaw cell, and 811 $\mathrm{mm}^2$ for the cell for porous media. The thick arrows indicate the flow direction.

Table 2. Summary of experimental runs and conditions. Constant parameters are the Damköhler number $Da=1600$ , the molecular diffusion coefficient $D_{{m}}=10^{-9}$ $\text{m}^2$ $\text{m}^2\ \text{s}^{-1}$ , and the characteristic length $L=d=2$ mm. For co-flow, $Pe=U L/D_{{m}}$ . For saddle flow, $U=\gamma L$ , leading to $Pe=\gamma L^2/D_{{m}}$ . For each flow configuration (columns) the different runs are indicated inside brackets.

3.1.3. Protocol

Before each experiment, we filled the cell with deionised water. For experiments with porous media, we injected $\mathrm{CO_2}$ gas before filling the cell with water. Hence any trapped $\mathrm{CO_2}$ gas dissolved in the water, and there were no remaining bubbles. Then we injected luminol and $\mathrm{H_2 O_2}$ 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/\nu$ , with $\nu =10^{-6}\ \text{m}^2\ \text{s}^{-1}$ the kinematic viscosity, $d=2$ mm the grain size, and $U=Q/A\phi$ , where $Q$ is the volumetric injection rate, $A$ is the cross-sectional area of the tank at injection, and $\phi$ 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). We measured this weak light emission as follows. We first mixed the luminol solution and the $\mathrm{H_2 O_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

(3.3) \begin{equation} s_R(x)=\left (\int y^2\, P(x,y)\,{\rm d}y - \left (\int y\, P(x,y)\,{\rm d}y\right )^2\right )^{1/2}, \end{equation}

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)\,{\rm d}y$ .

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 4 a). 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 5 a) 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\lt t_r$ ), probably because the cell geometry did not allow for a perfectly 2-D co-injection of reactants.

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.

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.

In the saddle-flow configuration (figure 4 b), 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 5 b). The data agree well with the scaling laws expected for the mixing-limited regime (2.14), which is expected since $\gamma ^{-1}\gt t_r$ for all compression rates used experimentally. The consistency of these observations with 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 6 b), 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.

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.

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.

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 7 a1). 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\gt 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).

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=10 d$ in figure 8. At this position, all experiments are in the mixing-limited regime. For the intensity of the reaction front intensity (figure 8 a1,b1) and the maximum reaction rate (figure 8 a2,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 8 a3,b3) is consistent with the scaling predicted by hydrodynamic dispersion, since this spatially integrated measure is not modified by the presence of pore-scale gradients, as also found in conservative experiments (Rousseau et al. 2023) .

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.

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

Due to stretching, the mixing front forms an ensemble of lamellae that are elongated in one direction and compressed in the other (figure 9 a). 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 ={\rm e}^{\lambda t}$ , with

(4.1) \begin{equation} \lambda \sim \frac {U}{d} \sim \frac {D_m}{d^2} Pe, \end{equation}

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

(4.2) \begin{equation} \frac {\textrm{d} s}{\textrm{d} t} = -\lambda s + \frac {D_{{m}}}{s}. \end{equation}

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),

(4.3) \begin{equation} s\to s_B = \sqrt {\frac {D_m}{\lambda }} \sim d \, Pe^{-1/2}. \end{equation}

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

(4.4) \begin{equation} t_m \sim \frac {\log Pe}{2 \lambda }, \end{equation}

which is obtained when the lamella width reaches the Batchelor scale by compression from the initial width, $s_0\, {\rm 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 9 a). 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 9 b). 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):

(4.5) \begin{equation} \rho _m=\rho (t_m)={\rm e}^{\lambda t_m} \sim \sqrt {Pe}, \end{equation}

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 9 b). 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:

(4.6) \begin{equation} j_A \sim C_A \lambda \sim C_A \frac {D_m}{d^2}\, Pe, \end{equation}

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

(4.7) \begin{equation} n = \beta \sqrt {Pe}, \end{equation}

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

(4.8) \begin{equation} J_A = n j_A = \beta \frac {D_m}{d^2} C_A\, Pe^{3/2}. \end{equation}

Initially, reactants are fully segregated, and $C_A$ and $C_B$ remain close to their initial concentrations $C_A^0$ and $q C_A^0$ . At a further distance downstream, called $x_\varDelta$ , dispersion inevitably smooths scalar gradients at the interface such that $C_A\lt 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\lt x_{\varDelta }$ and $C_A \approx C_A^0$ ( $C_B \approx q C_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:

(4.9) \begin{equation} \frac {{\rm d}m_{A|B}}{{\rm d}t}= -k C_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}, \end{equation}

where we used $C_B \approx q C_A^0$ for $x\lt x_{\varDelta }$ and for the general case of an asymmetric reaction front. This equation leads to

(4.10) \begin{equation} m_{A|B}=\beta \frac {C_A^0 D_{{m}} t_r }{q d^2}\, Pe^{3/2} \left ( 1-{\rm e}^{- q t/t_r} \right ) . \end{equation}

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

(4.11) \begin{equation} m_{B|A}= \beta \frac {C_A^0 D_{{m}} t_r q}{d^2}\, Pe^{3/2} \left ( 1-{\rm e}^{- t/t_r} \right ) . \end{equation}

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

(4.12) \begin{equation} I = k C_B^0 m_{A|B} + k C_A^0 m_{B|A} = \beta q \frac {C_A^0 D_{{m}}}{d^2}\, Pe^{3/2} \left ( 2-{\rm e}^{- t/t_r}-{\rm e}^{- q t/t_r} \right ). \end{equation}

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

(4.13) \begin{equation} 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 ). \end{equation}

This equation predicts a linear growth of reaction intensity followed by a plateau (figure 10 a). 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

(4.14) \begin{equation} m_{A|B} = J_A t, \end{equation}

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

(4.15) \begin{equation} 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}. \end{equation}

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

(4.16) \begin{equation} I = J_A = \beta q \frac {D_{{m}} C_A^0}{d^2}\, Pe^{3/2}. \end{equation}

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.

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}$ m $^2$ s⁻¹, $C_A^0=1$ 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).

The model of (4.13) is in relatively good agreement with experimental data for the range of distances and Péclet numbers considered (figure 10 a). 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 10 a), and the dependence of the reaction intensity and reaction maximum with the Péclet number (figure 8 a1,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.2 U/d$ measured in bead packs, Heyman et al. (2020) estimated that the Batchelor scale should reach the grain size for $Pe \leqslant 5$ . For the reaction fronts considered, this appears to occur at larger Péclet numbers, of the order of $Pe \leqslant 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 10 a). 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 $\varDelta$ 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 $\varDelta$ , 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\gt x_{\varDelta }$ , 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):

(4.17) \begin{equation} s = \sqrt {4 \alpha _T x}. \end{equation}

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

(4.18) \begin{equation} x_{\varDelta }=\frac {\Delta ^2}{4\alpha _T}. \end{equation}

Considering the observation $x_{\varDelta }\approx 20 d$ leads to $\varDelta \approx 2.4 d$ , 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

(4.19) \begin{equation} Pe_x \sim \frac {\gamma x d}{D_{{m}}}, \end{equation}

in (4.16), we obtain

(4.20) \begin{equation} I = \beta ' q \frac {C_A^0}{(d \, D_{{m}})^{1/2}}{(\gamma x)}^{3/2}. \end{equation}

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 10 b). 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 10 b) and Péclet number (figure 8 b1). 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 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.