# Resolving pore-scale concentration gradients for transverse mixing and reaction in porous media

Paiman Shafabakhsh, Tanguy Le Borgne, François Renard, Gaute Linga

► To cite this version:

Paiman Shafabakhsh, Tanguy Le Borgne, François Renard, Gaute Linga. Resolving pore-scale concentration gradients for transverse mixing and reaction in porous media. Advances in Water Resources, 2024, 192, pp.104791. 10.1016/j.advwatres.2024.104791. insu-04672449

HAL Id: insu-04672449

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

Submitted on 19 Aug 2024

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

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

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

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

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

# Resolving pore-scale concentration gradients for transverse mixing and reaction in porous media

Paiman Shafabakhsh<sup>a,*</sup>, Tanguy Le Borgne<sup>b,a</sup>, François Renard<sup>a,c</sup>, Gaute Linga<sup>a,d,e,*</sup>

## ARTICLE INFO

**Keywords:**

Mixing

Porous media

Reactive transport

## ABSTRACT

Mixing-limited reactions are central to a wide range of processes in natural and engineered porous media. Recent advances have shown that concentration gradients sustained by flow at the pore-scale influence macroscopic reaction rates over a large range of reactive transport regimes. Yet, resolving concentration gradients driven by fluid mixing at the pore-scale is challenging with current simulation methods. Here, we introduce a computational methodology to resolve concentration gradients at the pore scale in mixing-limited reactions. We consider a steady-state reactive transport problem characterized by reactive fluids flowing in parallel in a porous material. Given a mesh representation of the pore space and a steady velocity field, we solve the steady advection-diffusion equation for conservative scalar transport using a stabilized finite-element method combined with mesh refinement adapted to local scalar gradients. Based on this solution and assuming instantaneous reaction kinetics in the fluid, we infer the distribution of species involved in an irreversible bi-molecular reaction. We validate the method by comparing our results for uniform flow with analytical solutions and then apply it to simulate mixing-limited reactions in a three-dimensional random bead pack and Berea sandstone sample. Chaotic flow within the pore space leads to sustained concentration gradients, which are captured by our numerical framework. The results underscore the ability of the methodology to simulate transverse mixing and mixing-limited reactions in complex porous media and to provide bottom-up numerical data to improve the prediction of effective reaction rates at larger scales.

## 1. Introduction

Fluid flow through porous media controls the transport of dissolved chemicals, pollutants, and biological compounds in a wide range of natural and industrial porous systems (Dentz et al., 2011; Rolle and Le Borgne, 2019). The complex network of interconnected pores can lead to highly non-uniform flow patterns and, consequently, heterogeneous distributions of chemical concentrations (Dentz et al., 2023). Understanding the dynamics of mixing at the pore level and its influence on reaction kinetics on larger scales is important for applications such as the prediction of contaminant migration in groundwater (Rolle and Le Borgne, 2019), production of oil-wet reservoirs by CO₂ injection (Jiménez-Martínez et al., 2016), designing effective remediation strategies in subsurface environments (e.g., Kitanidis and McCarty, 2012), and the optimization of transport processes in engineered porous media (e.g., Johnson and Locascio, 2002).

Solute mixing brings reactants together in a fluid, allowing chemical reactions to occur. Since this process can be relatively slow in porous media, a range of chemical reactions have kinetics with characteristic times significantly shorter than typical mixing times (Li et al., 2006; Valocchi et al., 2019). Such mixing-limited reactions result in a reaction zone confined to an interface that acts as a boundary between segregated reactants (e.g., Rolle et al., 2009; Bauer et al., 2008; Willingham et al., 2008; Cirpka et al., 1999; Eckert et al., 2012; Hidalgo et al., 2015; Perez et al., 2019). A range of reactive transport problems are characterized by steady-state or slowly moving plumes, where mixing occurs predominantly transverse to the flow direction in narrow mixing fronts, typically on the millimeter scale (Rahman et al., 2005; Benekos et al., 2006; Shafabakhsh et al., 2024a; Nambi et al., 2003; Jose and Cirpka, 2004). The localization and size of such transverse mixing zones are greatly influenced by the heterogeneous nature of geological formations (Cirpka et al., 2015).

<footnote><sup>a</sup> The Njord Centre, Departments of Geosciences and Physics, University of Oslo, Norway</footnote>

<footnote><sup>b</sup> Geosciences Rennes, University of Rennes, France</footnote>

<footnote><sup>c</sup> ISTERRE, University Grenoble Alpes, University Savoie Mont Blanc, University Gustave Eiffel, CNRS, IRD, France</footnote>

<footnote><sup>d</sup> PoreLab, The Njord Centre, Department of Physics, University of Oslo, Norway</footnote>

<footnote><sup>e</sup> PoreLab, Department of Physics, Norwegian University of Science and Technology, Norway</footnote>

<footnote>* Corresponding author at: The Njord Centre, Departments of Geosciences and Physics, University of Oslo, Norway. E-mail addresses: paimans@uio.no (P. Shafabakhsh), gaute.linga@mn.uio.no (G. Linga).</footnote>
At Darcy scale, most numerical studies of reactive transport employ a continuum-based approach in which solute concentrations are integrated within a representative elementary volume (e.g., Katz et al., 2011; Molins and Knabner, 2019; Porta et al., 2013, 2016; Boso and Battiato, 2013). This approach uses hydrodynamic dispersion to represent the effect of unresolved pore-scale fluctuations in fluid velocities. Although this macrodispersion framework is useful in describing the spread of solute plumes, it tends to overpredict the rates of mixing-limited reactions because of the persistence of incomplete mixing at the pore scale (Raje and Kapoor, 2000; Gramling et al., 2002; Liedl et al., 2005). This highlights the need for methodologies that can accurately capture these pore-scale processes.

Substantial progress has been achieved to simulate the resulting pore-scale concentration gradients by using different Lagrangian and Eulerian simulation methods (Benson et al., 2017; Oostrom et al., 2016). Lagrangian methods, which involve tracking individual solute particles, have been applied in investigating mixing phenomena in single- and multiphase flow systems (Guida et al., 2012; Minier et al., 2014; Voronov et al., 2011; Sternagel et al., 2021; Klimenko, 2009; Benson and Bolster, 2016; Ding et al., 2017; Benson et al., 2019; Engdahl et al., 2019). A key benefit of Lagrangian methods is their ability to accurately capture advection-dominated transport processes, characterized by high Péclet number Pe (the ratio of diffusive to advective fluxes), without introducing numerical dispersion (Salamon et al., 2006; Boso et al., 2013). However, the limitation is the requirement for a large number of particles to effectively resolve concentration gradients, leading to high computational costs (Noetinger et al., 2016). Moreover, interpolation methods for computing concentrations based on particle positions introduce a risk of over-smoothing concentration fields, potentially degrading accuracy (Fernàndez-Garcia and Sánchez-Vila, 2011). Eulerian simulation methods, on the other hand, solve the advection–diffusion equation directly in a discretized domain Sole-Mari et al. (e.g., 2022), Baek and Seo (e.g., 2017). To avoid spurious oscillations in the resulting concentration fields, various stabilization approaches have been proposed (Hughes, 1979; Brooks and Hughes, 1982). These may again introduce non-physical numerical dispersion and smear out concentration gradients unless the spatial and temporal resolution is high enough.

While mixing dynamics can be well resolved at the pore scale in two-dimensional simulations (e.g., Acharya et al., 2007; Du et al., 2023; Hejazi and Azaiez, 2013; Oostrom et al., 2016), this task is much more challenging in three-dimensional porous media (e.g., Sole-Mari et al., 2022). The additional degree of freedom induces chaotic stretching and folding at the pore scale, which sustains concentration gradients and incomplete mixing (Lester et al., 2014, 2016; Heyman et al., 2020; Souzy et al., 2020; Heyman et al., 2021; Aquino et al., 2023; Sanquer et al., 2024). Therefore, the transition towards three-dimensional simulations is essential to capture the spatial variability and heterogeneity of porous media, which cannot be adequately represented in two-dimensional models. The effects of chaotic mixing on reactive processes are best evidenced when considering transverse mixing in steady state mixing fronts (e.g. (Sanquer et al., 2024)). Yet, there is, to our knowledge, currently no method or studies that directly address how to properly and practically resolve chaotic mixing in 3D pore-scale flows. Eulerian methods suffer from numerical dispersion and numerical instabilities, while Lagrangian methods require a large amount of particles to accurately resolve concentration gradients in three dimensions. Methods that rely on temporal integration of the governing equations (both Lagrangian and Eulerian) require a long time to converge to the steady state of transverse mixing between two continuously injected solutions. Using time-dependent solvers such as OpenFoam, with a spatial resolution $\Delta x$ and time step $\Delta t$, requires $\Delta x \sim \Delta t \sim \text{Pe}^{-1/2}$ for stability as discussed e.g. by Martínez-Ruiz et al. (2018). Considering that the diffusion time in the transverse direction is $t_D \sim t_D/\Delta t \sim \text{Pe}^{3/2}$ and problem size $N \sim \Delta x^{-3}N_f \sim \text{Pe}^3$, resulting in high computational

costs for advection-dominated transport. Hence, there is a critical need for methodologies that can balance computational efficiency with the ability to accurately capture the pore-scale gradients for simulating mixing and reactions in 3D porous media.

To effectively resolve sharp concentration gradients, adaptive mesh refinement strategies have been widely used for the transport of solutes in porous media (Mansell et al., 2002), primarily in the context of transient flow at the Darcy scale (Nilsson et al., 2005; Pau et al., 2009; Cusini et al., 2019; Dell'Oca et al., 2018). Refinement criteria may be based on both local gradients (Cusini et al., 2019) and a *posteriori* error analysis (Dell'Oca et al., 2018). For advection-dominated conditions, robust refinement strategies based on the Streamline-Upwind Petrov-Galerkin method (Hughes, 1979; Brooks and Hughes, 1982) have been developed (John and Novo, 2013). These techniques may provide new opportunities to resolve sharp chemical gradients in complex three-dimensional flow fields that arise at the pore scale in realistic porous media.

Here, we introduce a new methodology for simulating steady-state mixing and reaction at the pore-scale in three-dimensional porous media. The approach relies on accurately solving steady conservative scalar transport in a porous domain, given a mesh representation of the pore space and a steady velocity field within it. To simultaneously achieve numerical stability and physically valid concentration fields, a key step is the combination of numerical stabilization (using the Streamwise-Upwind Petrov-Galerkin approach) and iterative refinement based on local concentration gradients. This methodology offers enhanced spatial resolution to capture detailed transverse mixing processes, and accurate representation of concentration gradients, which the existing modeling frameworks fail to adequately resolve at the pore-scale. Based on the proposed refinement criterion, we provide scaling estimates for the number of computational nodes needed to resolve the mixing zone, which is qualitatively different in two and three dimensions. We evaluate instantaneous reactions by post-processing the resulting conservative concentration fields. The adaptability of this methodology allows us to resolve pore-scale mixing in a wide range of porous media systems, including chaotic flow through a bead pack and flow in a natural rock sample of Berea sandstone.

## 2. Physical model

In the following, we present the considered physical model and the modeling assumptions used. We describe both conservative and reactive transport, as well as to how the two modes of transport can be linked assuming instantaneous reaction kinetics.

### 2.1. Porous domain and velocity field

We consider a fluid domain $\Omega \subset \mathbb{R}^3$, schematically displayed in Fig. 1, which represents the pore space within a solid porous matrix. This pore space is filled with single-phase fluid. $\Omega$ is contained within a box of size $L_x \times L_y \times L_z$. The boundary of $\Omega$ is comprised of the inlet boundary $\Gamma_{\text{in}}$ (at $z=0$, ref. Fig. 1), the outlet boundary $\Gamma_{\text{out}}$ (at $z=L_z$), and the pore wall boundary $\Gamma_{\text{wall}}$, i.e. the surface of $\Omega$ that separates fluid and solid.

A steady fluid velocity field $\mathbf{u}(\mathbf{x})$ is defined within $\Omega$. Here, $\mathbf{x} = (x, y, z) \in \Omega$ is the spatial coordinate and the velocity field is that of an incompressible fluid, such that $\nabla \cdot \mathbf{u} = 0$. Furthermore, we assume that the solid is impermeable, i.e. $\mathbf{u} \cdot \hat{\mathbf{n}} = 0$ for $\mathbf{x} \in \Gamma_{\text{wall}}$. Here, $\hat{\mathbf{n}}$ the outward normal of $\Omega$. For creeping flows typical to porous media, the velocity field $\mathbf{u}(\mathbf{x})$ will generally be well represented by the solution of the Stokes equations (Feder et al., 2022) (see Appendix B).
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Schematic representation of the single-phase fluid domain $\Omega$ within a solid porous matrix, and its boundaries. Here, the two fully miscible single-phase solutions A (red solution) and B (blue solution) are separated by a sharp interface, and are co-injected from the surface $\Gamma_{in}$. This inlet condition corresponds to the Heaviside step function Eq. (16) for the conserved composite species $\delta$, as indicated. $\Gamma_{wall}$ is the boundary surface that separates fluid and solid, and $\Gamma_{out}$ is the surface at the bottom.</figcaption></figure>

## 2.2. Advection–diffusion equation

We now consider the concentration $c_i$ of a given chemical species $C_i$ dissolved in a fluid that flows through a porous medium. The evolution of $c_i(x, t)$ in space $x \in \Omega$ and time $t$, is described at the pore-scale by the advection–diffusion equation

$$ \frac{\partial c_i}{\partial t} + \mathbf{u} \cdot \nabla c_i - D_i \nabla^2 c_i = 0, \quad (1) $$

where $D_i$ is the diffusivity of the species. A Dirichlet boundary condition is supplied at the inlet boundary as

$$ c_i(x, t) = C_{i,0}(x) \quad \text{for } x \in \Gamma_{\text{in}}, \quad (2) $$

where the precise function $C_{i,0}(x)$ will be specified later. For the pore walls, we also take the boundaries to be impermeable for the chemical species, which when $\hat{\mathbf{n}} \cdot \mathbf{u} = 0$ is expressed by $\hat{\mathbf{n}} \cdot \nabla c_i = 0$ for $x \in \Gamma_{wall}$. To close the equations, we also choose the boundary condition $\hat{\mathbf{n}} \cdot \nabla c_i = 0$ for $x \in \Gamma_{out}$, in this case representing no diffusive solute flux at the outlet.

For constant injection rates, a steady transport regime will eventually be established. Here we consider a co-injection of two fluids with different solute concentrations (e.g., Acharya et al., 2007) leading to a non-trivial steady-state concentration $c_i(x)$ described by removing the time-dependence in Eq. (1):

$$ \mathbf{u} \cdot \nabla c_i - D_i \nabla^2 c_i = 0, \quad \text{for } x \in \Omega, \quad (3a) $$

$$ c_i(x) = C_{i,0}(x), \quad \text{for } x \in \Gamma_{\text{in}}, \quad (3b) $$

$$ \hat{\mathbf{n}} \cdot \nabla c_i = 0, \quad \text{for } x \in \Gamma_{\text{wall}} \cup \Gamma_{\text{out}}. \quad (3c) $$

## 2.3. System of advection–diffusion–reaction equations

When chemical reactions take place within the fluid, the concentration fields $c_i$ of the different species $C_i$ become coupled to each other. Eq. (1) then generalizes to the advection–diffusion–reaction equation:

$$ \frac{\partial c_i}{\partial t} + \mathbf{u} \cdot \nabla c_i - D_i \nabla^2 c_i = R_i. \quad (4) $$

Here, the chemical species $C_i$, $i \in \{1, \dots, N\}$ interchange mass through the reaction terms $R_i$. In principle, these terms could result from several

reactions taking place in a complex reaction network (De Simoni et al., 2007a). However, to demonstrate our approach, we consider a simple bimolecular reaction for which pore-scale mixing plays an important role (Valocchi et al., 2019).

## 2.4. Irreversible bi-molecular reaction

The simplest fluid–fluid reaction involving more than one species is the irreversible bimolecular reaction:

$$ A + B \rightarrow C. \quad (5) $$

This reaction has been well studied in porous media both in simulations (e.g., Borgne et al., 2014; Kim, 2019) and in experiments (e.g., Willingham et al., 2008; Izumoto et al., 2023; d. Anna et al., 2014). Here, mixing between two reactants A and B yields product C upon contact between the two miscible fluids. The local reaction rate expressing the creation of product C can be expressed by:

$$ R = k_{ab}, \quad (6) $$

where $k$ is the reaction rate constant. The steady advection–diffusion–reaction system corresponding to reaction (5) can be described by Eq. (4) and becomes:

$$ \mathbf{u} \cdot \nabla a - D \nabla^2 a = -R \quad (7a) $$

$$ \mathbf{u} \cdot \nabla b - D \nabla^2 b = -R \quad (7b) $$

$$ \mathbf{u} \cdot \nabla c - D \nabla^2 c = R, \quad (7c) $$

where $a$, $b$, and $c$ are the concentrations of the involved species, and we have taken the diffusion coefficient $D$ to be equal for all species.

### 2.4.1. Inlet boundary conditions

We consider the idealized scenario of two miscible fluids containing, respectively, species A and B, which mix within the porous domain $\Omega$ (see Fig. 1). This scenario fixes the inlet boundary condition (2) for the reaction (5). Both species have the same initial concentration (normalized to 1, for simplicity), and the fluids are injected from two different locations at the inlet boundary $\Gamma_{in}$ (i.e., co-flow injection). The concentration of species C is equal to zero at the inlet boundary, meaning that all production of C by construction occurs within the domain. The boundary $\Gamma_{in}$ can be split into two subdomains separated by a sharp interface, where either

$$ (a, b, c) = (1, 0, 0) \quad \text{or} \quad (0, 1, 0). \quad (8) $$

In particular, in this work, we will approximate the inlet boundary condition

$$ a(x, y, 0) = \Theta(x), \quad (9) $$

$$ b(x, y, 0) = 1 - \Theta(x), \quad (10) $$

where $\Theta(x)$ is the Heaviside step function,

$$ \Theta(x) = \begin{cases} 0 & \text{for } x < 0, \\ 1 & \text{for } x \ge 0. \end{cases} \quad (11) $$

### 2.4.2. Simplifications by identifying conserved and non-conserved quantities

From Eq. (7), we identify two conserved quantities $\delta$, $\psi$ and a non-conserved quantity $\rho$:

$$ \delta = a - b, \quad (12a) $$

$$ \rho = a + b, \quad (12b) $$

$$ \psi = c - \rho/2 \quad (12c) $$

The reaction rate (6) can then be expressed as

$$ R = \frac{k}{4} (\rho^2 - \delta^2). \quad (13) $$
Accordingly, Eqs. (7a) and (7b) can be expressed as:

$$ \mathbf{u} \cdot \nabla \delta - D \nabla^2 \delta = 0, \qquad (14a) $$

$$ \mathbf{u} \cdot \nabla \rho - D \nabla^2 \rho = \frac{k}{2} (\delta^2 - \rho^2), \qquad (14b) $$

$$ \mathbf{u} \cdot \nabla \psi - D \nabla^2 \psi = 0. \qquad (14c) $$

The evolution Eqs. (14a) and (14c) for $\delta, \psi$ are linear and do not depend on any other field but themselves, and can be solved independently. Further, the evolution Eq. (14b) for $\rho$ depends only on the found solution for $\delta$. Hence, Eqs. (14a)–(14c) can be solved sequentially, with the only constraint that Eq. (14b) must be solved after Eq. (14a). This approach is computationally beneficial, as it reduces the size of the linear system resulting from a finite element discretization (see Section 3.2) compared to a coupled solution of Eq. (7).

Following Eqs. (8) and (12), we find for the composite concentrations at the inlet boundary:

$$ \delta \in \{-1, 1\}, \quad \rho = 1, \quad \psi = 1/2 \quad \text{for } \mathbf{x} \in \Gamma_{\text{in}}. \qquad (15) $$

For the specific inlet conditions considered herein, which divides the inlet into two halves along the $x$ direction, we obtain from Eq. (10):

$$ \delta(x, y, 0) = 2\Theta(x) - 1, \qquad (16) $$

where $\Theta$ is the Heaviside step function defined in Eq. (11).

Eq. (14a) shows that $\delta$ will stay within the interval $[-1, 1]$ in the entire domain, while Eq. (14b) shows that $\rho$ will tend to decay with the distance away from the inlet boundary (as $\rho \ge |\delta|$). Finally, with the inlet condition $\psi = 1/2$, the only solution that can satisfy Eq. (14c) is a trivial one, $\psi = 1/2$ everywhere in the domain. Hence, the concentration of product C can be calculated explicitly by solving

$$ c = \frac{1 - \rho}{2}, \qquad (17) $$

and henceforth Eq. (14c) can be neglected.

### 2.4.3. Instantaneous reaction

Mixing-limited reactions are characterized by a high Damköhler number (Da), defined as the ratio between the characteristic time of fluid motion to that of the chemical reaction. Considering a reaction with a first-order-kinetic constant ($k$), the Damköhler number is

$$ Da = \frac{kd^2}{D} \qquad (18) $$

where $d$ is a characteristic pore size. We now assume that the reaction occurs on time scales much smaller than the diffusive and advective time scales, which corresponds to the limit of large Damköhler number $Da \to \infty$, i.e. instantaneous reaction. Taking $k \to \infty$ in Eq. (14b), we then obtain:

$$ \rho = |\delta|. \qquad (19) $$

The concentration of reactants A, B, and product C can be readily calculated using Eqs. (12) and (19) as:

$$ a = \frac{|\delta| + \delta}{2}, \qquad (20a) $$

$$ b = \frac{|\delta| - \delta}{2}, \qquad (20b) $$

$$ c = \frac{1 - |\delta|}{2}. \qquad (20c) $$

Thus, given that the reaction (5) is instantaneous and irreversible, the system expressed by Eq. (7) is completely specified once the solution $\delta$ of Eq. (14a) is known. This property will be exploited in the ensuing sections; in particular, we will solve the steady advection-diffusion Eq. (14a) and use Eq. (20) to calculate the distribution of the species involved in the reaction (5).

### 2.5. Metrics of mixing and reaction

We investigate mixing and reaction at different Péclet numbers, which characterizes the ratio of the diffusive and advective time scales at the pore scale,

$$ Pe = \frac{Ud}{D}, \qquad (21) $$

where $d$ is a characteristic pore size and $U$ the average flow speed. Fluid mixing can be quantified by measuring different quantities related to the concentration field, flux, transported mass, and reaction rate of the product, as explored in previous studies (e.g., de Dreuzy et al., 2012; Rolle et al., 2009; d. Anna et al., 2014; Alhashmi et al., 2015). Here we summarize the main quantities used for measuring mixing and reaction in the present study.

For quantifying conservative mixing, we use the scalar dissipation rate, which has been used in several studies of fluid mixing in porous media as a proxy for the mixing rate (e.g., Luo et al., 2008; Le Borgne et al., 2010; Bolster et al., 2010; Chiogna et al., 2011; Jha et al., 2011; Engdahl et al., 2013). Using $\delta$ as a conserved quantity, the conservative scalar dissipation rate is defined as

$$ I = \int_{\Omega_{\perp}(z)} D|\nabla \delta|^2 dxdy, \qquad (22) $$

where the integral is taken over a cross section $\Omega_{\perp}(z)$ of $\Omega$.

We quantify effective reaction rates using the total reaction rate in each cross-section. The latter can be estimated as (see Appendix D),

$$ R_c = \frac{dJ_c}{dz}. \qquad (23) $$

where $J_c$ is the total flux of $c$ through a cross-section $\Omega_{\perp}(z)$ of $\Omega$,

$$ J_c(z) = \int_{\Omega_{\perp}(z)} \hat{\mathbf{z}} \cdot \mathbf{j}_c \, dxdy, \qquad (24) $$

and the local flux $\mathbf{j}_c$ is the sum of the advective and diffusive fluxes,

$$ \mathbf{j}_c = \mathbf{u}c - D\nabla c. \qquad (25) $$

The total produced mass of C in each cross-section can be calculated as

$$ M_c = \int_{\Omega_{\perp}(z)} c(x) \, dxdy, \qquad (26) $$

which is directly proportional to the evolving volume of the mixing zone.

## 3. Computational method

### 3.1. Overall workflow

The workflow consists of four distinct steps, as shown schematically in Fig. 2. The prerequisites for the workflow are a computational mesh representing the fluid domain and a well-resolved velocity field within it. Specifically, we assume that we have a three-dimensional computational mesh representation of the pore space $\Omega$ (as discussed in Section 2.1 and shown in Fig. 1). Computing pore-scale velocity fields is in itself an active research topic and is not the focus of this paper, and we refer to other resources for more specialized literature on this topic (Indelman, 2001; Siena et al., 2015; Aramideh et al., 2018). In this work, we primarily computed the velocity field by solving the Stokes equations with no-slip conditions on the solid boundaries and slip conditions on the lateral side walls (see Appendix B).

The main step of the workflow is to solve the advection-diffusion Eq. (14a) to calculate the concentration field of conserved species $\delta$ (Fig. 2). As a consequence of the nonuniform concentration field that we impose on the inlet boundary in the advection-diffusion equation, we obtain nontrivial steady concentration distributions within the porous media. These distributions form the basis for subsequent analyses. Accurate determination of $\delta$ is achieved through an iterative
<figure><img src="image_5.png" /><figcaption>Fig. 2. Workflow for simulation of transverse mixing and reaction in a porous geometry. (a) Starting with a computational mesh and a steady velocity field, (b) we solve the advection–diffusion equation to compute the concentration field of the conserved species δ. (c) An optional step of mesh refinement adapted to the mixing zone (interface between reactive fluids) before (b) can be repeated iteratively. (d) Finally, we may compute spatial distributions of chemical species by post-processing analysis.</figcaption></figure>

adaptive refinement procedure, which is described in Section 3.3. Finally, in the post-processing phase, the conserved scalar field δ is used to compute the spatial distribution of chemical species involved in Eq. (5) within the porous media by solving the local Eq. (20).

## 3.2. Finite element scheme

The finite element method (FEM) offers a powerful approach for solving partial differential equations like Eq. (14a). In the limit of large Pe, Eq. (14a) becomes singularly perturbed (Roos, 2008) and naive discretization approaches will typically fail to converge or lead to un-physical oscillatory solutions unless very fine meshes are used (Elman et al., 2014). To obtain meaningful and convergent results, various stabilization techniques have been proposed (Roos, 2008). One such technique is the Streamwise-Upwind Petrov–Galerkin method (Brooks and Hughes, 1982; Hughes, 1979), which modifies the standard Galerkin finite element method by adding numerical diffusion in the streamline direction (John et al., 2018; Elman et al., 2014).

The variational formulation corresponding to Eq. (14a) with Streamwise-Upwind Petrov–Galerkin stabilization is given by the following. Find $\delta \in \mathcal{V}$ such that for all test functions $v \in \mathcal{V}$:

$$ \int_{\Omega} (\mathbf{u}_h \cdot \nabla \delta) v \, d^3 \mathbf{x} + \int_{\Omega} D \nabla \delta \cdot \nabla v \, d^3 \mathbf{x} + \int_{\Omega} \tau(\mathbf{u}_h \cdot \nabla \delta)(\mathbf{u}_h \cdot \nabla v) \, d^3 \mathbf{x} = 0. \quad (27) $$

Here, we consider the function space $\mathcal{V} \subset H^1(\Omega)$ and $\tau$ is a space-dependent stabilization parameter. The latter is modeled according to Elman et al. (2014, Eq. 3.44):

$$ \tau = \frac{h}{2|\mathbf{u}_h|} \Theta(1 - \mathrm{Pe}_h^{-1}), \quad (28) $$

where $h$ represents the local mesh size, $|\mathbf{u}_h|$ is the local magnitude of the velocity vector $\mathbf{u}_h$ interpolated from the input velocity field $\mathbf{u}$, and $\Theta(x)$ is the Heaviside step function defined in Eq. (11). Further, the element Péclet number is defined by $\mathrm{Pe}_h = |\mathbf{u}_h|h/(2D)$.

As discussed below in Section 3.3, $\mathbf{u}$ is typically defined on a coarser mesh than $\mathbf{u}_h$. In practice, we solved $\mathbf{u}$ once with a higher-order ($P_2$) finite element basis and used a piecewise constant representation ($P_0$) of $\mathbf{u}_h$, which we found to give the most robust convergence. With adaptive mesh refinement in the interface region, the $P_0$ representation $\mathbf{u}_h$ on the finer mesh approaches the $P_2$ representation $\mathbf{u}$ on the coarser mesh. In our simulations, we used piecewise linear ($P_1$) elements for $\delta$ for robustness. Note that the form of Eq. (27) assumes this; for higher-order elements, a diffusive term must be added to Eq. (27) (Elman et al., 2014).

## 3.3. Refinement criteria

To obtain accurate solutions of the scalar field $\delta$, it is critical to combine the stabilization discussed above with mesh refinement. In our case, the interface between the fluids is of particular interest as it will form a layer of sharp change in $\delta$, and to extract physically valid results it is imperative to add as little numerical diffusion as possible to this zone. Thus we need a refinement strategy where we identify the regions where $\delta$ exhibits significant variations and iteratively refine this region. The large fraction of the mesh where $\delta$ is relatively constant, even though the velocity field $\mathbf{u}$ is complex, may remain coarse. In the regime of pore-scale flows we are interested in, corresponding to the Batchelor regime (Haynes and Vanneste, 2005; Heyman et al., 2023), the length scale of velocity fluctuations is typically larger than the length scale of concentration variations. Hence, we only need to refine the mesh used for computing $\delta$, which justifies computing $\mathbf{u}$ only once and using an interpolated $\mathbf{u}_h$ (on a progressively finer mesh) as discussed in Section 3.2.

We adopt a simple refinement criterion employed for mesh refinement based on the product of the local mesh size $h$ and the magnitude of the concentration gradient, $|\nabla\delta|$, expressed as

$$ h|\nabla\delta| < \Delta_{\delta}. \quad (29) $$

where $\Delta_{\delta}$ is a user specified tolerance that limits how much $\delta$ can change over an element. In the simulations presented in Section 4, we used a tetrahedral mesh representation of the domain and chose $h$ as the longest edge in the associated tetrahedron. Every cell of the mesh wherein the criterion Eq. (29) is not met, is marked for refinement and subsequently divided according to the algorithm of Plaza and Carey (1996) as implemented in FEniCS (Alnæs et al., 2015) (see Section 3.4). We found the value $\Delta_{\delta} = 0.2$ to strike a good balance between computational efficiency and the accuracy of the solution, effectively placing most of the degrees of freedom in regions where the concentration field changes drastically over short length scales.

## 3.4. Implementation

We solve the advection–diffusion Eq. (14a) in the form of Eq. (27) using the FEniCS/DOLFIN framework (Alnæs et al., 2015). FEniCS is a software collection for automating the solution of differential equations using the FEM, and DOLFIN functions as the primary user interface to FEniCS. Our FEniCS-based solver, named AddICTIF,<fnref n="1" /> is implemented

<footnote n="1">¹ Advection-Diffusion-Chemistry in a Time-Independent Framework.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Flowchart of the steps to use the software package AdDiCTIF for steady-state simulation of transverse mixing and reaction.</figcaption></figure>

in Python and C++, and available for use and further development in a dedicated Git repository (Shafabakhsh and Linga, 2024). The code is structured into distinct modules as shown in Fig. 3 for easy code reuse:

* The module `ade_steady` solves the advection-diffusion equation given the mesh and velocity field. To solve the linear system resulting from discretizing Eq. (27), we use the (iterative) bi-conjugate gradient method (bicgstab), with the Hypre Euclid implementation of the ILU preconditioner (hypre_euclid).

* After solving the advection-diffusion equation and obtaining the concentration field of $\delta$, the mesh may be refined by using the module `refine`. This task may be repeated by executing `ade_steady` and `refine` iteratively to progressively refine the solution. The iteration is stopped at the latest when the criterion Eq. (29) is met everywhere, and thus no cells are marked for refinement.

* Subsequently, post-processing steps can be performed to analyze and visualize the results. With the module `postprocess_abc`, we compute the concentration fields of the components `a`, `b`, and `c`, based on sets of Eq. (20).

* Following this step, the module `analyze_data` interpolates the solution into a Cartesian grid. This step enables easy analysis of different cross-sections.

* Lastly, with module `compute_averages` we can compute the averages of different indicators presented in Section 2.5 in each cross-section.

## 4. Numerical simulations

To assess the performance of the method, we apply it to three representative cases. To establish the validity of the approach, we first consider the case of uniform flow, where approximate analytical solutions can be obtained. This benchmark reveals differences between refinement when going from two to three dimensions. Second, we consider mixing and reaction through a three-dimensional monodisperse spherical bead pack, which is a well-studied model system for porous granular media (Heyman et al., 2020). Lastly, we consider mixing and reaction in a three-dimensional natural porous rock sample. These simulations highlight the ability of the model to effectively resolve and quantify mixing in a broad range of domains. Notably, variations in pore-scale geometry are known to impact the mixing dynamics, making such evaluations crucial.

### 4.1. Uniform flow

To establish the convergence of the method with mesh refinement, we first consider the case of uniform flow, for which approximate solutions are available (Abramowitz et al., 1968; Sale and McWhorter, 2001; Kumar et al., 2009). For clarity, we give a brief derivation of these solutions in the Appendix A. The uniform flow field is given by

$u = U\hat{z}$, where $U$ is a constant. To numerically approximate the inlet boundary condition Eq. (16) on a discrete mesh, we use

$$ \delta(x, y, 0) = \operatorname{erf}\left(\frac{x}{\epsilon}\right), \qquad (30) $$

where $\epsilon$ is a smoothing length. In our simulations, we used $\epsilon = 10^{-3}$.

As presented in Appendix A, for a simple one-dimensional diffusion equation, where $z$ assumes the role of a time coordinate, the solution subject to Eq. (30) is given by

$$ \delta \approx \operatorname{erf}\left(\frac{x}{\sqrt{4\epsilon z}}\right), \qquad (31) $$

where $\alpha = D/U$ and the approximation is valid for $\epsilon \to 0$.

We consider a cubic domain of side lengths $L = 1$, with an initial mesh contained 2400 nodes. In Fig. 4a and Fig. 4b, we show the calculated $\delta$ obtained for the uniform flow case when $\alpha = 0.01$ in a refined mesh containing 489,024 nodes, achieved after eight refinement iterations. The comparison between different iterations of refinement is presented below. With the choice of $\alpha$ and system size, we expect longitudinal diffusion to be negligible everywhere except for a small neighborhood around the inlet ($x = 0$) so that the analytical expression (31) above should be applicable.

Fig. 4a shows the steady-state concentration distribution of conserved species $\delta$. In Fig. 4b, the concentration profiles along the $x$-direction at various heights are shown. The numerically obtained profiles (indicated by points) are compared to the approximate theory (31) (represented by continuous lines). The agreement between the numerical data and theory is very good, which serves as a first validation of our simulations.

As a second test, we consider the conservative scalar dissipation rate (Le Borgne et al., 2010) which is approximated in Appendix A for the depth-averaged expression

$$ I_x(z) \approx \frac{2\sqrt{2D}}{\sqrt{\pi(4Dz + \epsilon^2)}}. \qquad (32) $$

In Fig. 4c, we compare the calculated depth-averaged dissipation rate for the uniform flow scenario with $\alpha \in \{0.01, 0.0025, 0.001, 0.00025\}$ for both our numerical models (represented by points) and the analytical solution (32) (represented by continuous lines). The numerical data closely follow the theoretical curves, particularly as $\alpha$ is lowered. The agreement between the model and the analytical solution can be attributed to the fact that under conditions of low diffusivity (characterized by low $\alpha$ values), the longitudinal component of diffusion has a vanishing contribution to the transport, in accordance with the assumptions made in the analytical calculations above.

To achieve the desired level of precision, we performed an iterative refinement process until no element were marked for refinement according to our selected criterion Eq. (29). For all values of $\alpha$, convergence was reached after eight iterations, as shown in Fig. C.13a. The detailed refinement process and its dependence on the smoothing scale of the inlet profile, as well as the differences observed between two and three dimensions, are discussed in detail in Appendix C.
<figure><img src="image_7.png" /><figcaption>Fig. 4. (a) Numerically obtained steady-state concentration distribution of conserved scalar $\delta$ for uniform flow ($\alpha = 0.01$). (b) The concentration profiles along the x-direction at different heights $z$ comparing numerical data (dots) and the approximate analytical expression (lines, Eq. (31)). (c) Dissipation rate integrated over the cross-section for uniform flow computed from numerical simulations (dots) compared to the analytical solutions (lines, Eq. (32)) for different dispersivities $\alpha$.</figcaption></figure>

Fig. 5a shows the convergence of the numerical dissipation rate towards the analytical value with different iterations for the case of $\alpha = 0.001$, demonstrating the convergence of our numerical solution towards the analytical solution with mesh refinement. We note that the points for iterations larger than four are virtually indistinguishable, despite an order of magnitude smaller mesh size. This result is due to excessive refinement near the inlet having negligible effects downstream, and suggests that our refinement criterion could be limited, for example, by increasing the inlet smearing scale $\epsilon$ or stopping the refinement at a specified smallest mesh size $h_{min}$.

In Fig. 5b and Fig. 5c, we analyze the cross-sectional properties of the product C for uniform flow with $\alpha = 0.001$ for different mesh refinement iterations. The integrated mass, $M_c$, is shown in Fig. 5b and Fig. 5c shows the reaction rates $R_c$ calculated by Eq. (23). Both cases are shown to approach the analytical solutions presented in Appendix E. The results confirm the importance of sufficient refinement, showing that in this case sufficient resolution is achieved after around 4 refinement iterations. A slight decrease is noticeable near $z = 1$ in the numerical solution for reaction rate compared to the analytical solution (Fig. 5c), which we interpret as a relic from an initially coarse mesh, where the refinement procedure fails to identify that further refinement is needed near the outlet.

The excellent match between the numerical results of this section and analytical solutions under uniform flow conditions shows the robustness of our approach in resolving concentration gradients. This validation supports expanding our simulations to more complex porous media.

## 4.2. Flow through a bead pack

We investigate conservative and reactive transport through a random bead pack, which is an archetypical example of a granular porous medium (Heyman et al., 2020; Turuban et al., 2019). The bead pack was generated using the discrete element method code Yade (Smilauer et al., 2023) to approximate a random close packing. 149 beads of radius $d = 0.2$ were placed inside a periodic box (in order to suppress effects that may cause order in the bead pack), which was shrunk until all the beads were locked in place and the linear size of the box was $L = 1 \pm 0.001$. This procedure resulted in a porosity of 38%. To avoid resolving extremely small elements near the cusp-shaped contacts between beads (Turuban et al., 2019), smaller beads of diameter $0.2d$ were also inserted at bead contacts, mainly overlapping with the large beads. The pore space between these grains was then meshed with tetrahedral elements using gmsh (Geuzaine and Remacle, 2024), resulting in a mesh containing approximately 465,000 nodes.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Effect of mesh refinement for the case of uniform flow, with $\alpha = 0.001$. (a) Convergence of the numerical dissipation rate towards the analytical value with each iteration. (b) Mass per length of product C for different iterations of mesh refinement. (c) Integrated reaction rate of C for different iterations.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. (a) Steady-state velocity field in a bead pack computed from the Stokes equations Appendix B. (b) The distribution of conserved species $\delta$ within the bead pack based on the velocity field, with increasing Péclet number, Pe, from left to right. (c) The concentration $c$ of product C resulting from a first-order bimolecular reaction (5), computed using the conservative solution of $\delta$, with the same Pe values as above.</figcaption></figure>

To approximate the sharp inlet boundary condition (16) by a smooth function to avoid over-refinement in the inlet region, we used

$$ \delta(x, y, 0) = \tanh \left( \frac{x}{\sqrt{2\epsilon}} \right), \qquad (33) $$

with $\epsilon = 10^{-2}$. This profile closely resembles that used for the uniform flow (Eq. (30)). We verified that results are insensitive to the choice of inlet profile as long as the gradient at the inlet is significantly higher than gradients downstream.

Fig. 6a shows the steady velocity field obtained by solving the Stokes equations (see Appendix B) within the domain. We computed the concentration distribution of conservative species (Fig. 6b), and the reactive product (Fig. 6c) for $Pe \in \{1.3 \cdot 10^1, 1.3 \cdot 10^2, 1.3 \cdot 10^3\}$, estimated from Eq. (21) where the characteristic length scale is taken as the bead diameter $d$ and the mean velocity $U$ is normalized to 1. From Fig. 6b it is clear that the mixing region within the bead pack, indicated by intermediate values of $\delta$, narrows as Pe increases. This result highlights the relatively higher influence of diffusion at lower Pe.

A representative part of (a) the initial mesh and (b) the mesh at the final iteration for the flow in a sub-sample of the bead pack is displayed in Fig. C.14, showing how mesh refinement captures the narrow mixing zone. Fig. 7a shows the evolution of the total number of computational degrees of freedom (mesh nodes) as a function of mesh refinement iteration. For each value of Pe, the mesh refinement process culminates at a specific iteration, i.e. when the criterion (29) is met in all cells. The refinement stops after three, five, and six iterations for the three considered Péclet numbers in increasing order. The number of nodes at

which the refinement process saturates ($N_{sat}$) as a function of Pe values is shown in Fig. 7b, consistent with a linear scaling $N_{sat} \sim Pe$.

This relation can be rationalized in light of the Batchelor scale that sets the mixing scale in chaotic flows (Heyman et al., 2020). Stretching and folding of fluid elements in the pore space lead to an exponential elongation of fluid elements (Lester et al., 2016; Turuban et al., 2019; Souzy et al., 2020) $\rho \sim e^{\lambda t/t_a}$, where $\lambda$ is the dimensionless Lyapunov exponent and $t_a = d/U$ is the advection time. The resulting constant stretching rate implies that diffusion and compression equilibrate at the Batchelor scale $s_B \sim \sqrt{Dd/(\lambda U)}$. Thus, the width of the mixing interface is expected to scale with $s_B \sim dPe^{-1/2}$. As detailed in Appendix G, a fixed number of nodes $N_{\perp}$ is required to satisfy the refinement criterion Eq. (29) in the direction perpendicular to the interface, leading to $h \sim s_B/N_{\perp}$. Assuming roughly isotropic elements yields a total number of nodes $N_{sat} \sim N_{\perp}A/h^2$, where $A$ is the area of the interface, which for large Pe will only increase weakly with Pe. Thus, we find $N_{sat} \sim N_{\perp}^3Pe$, in particular $N_{sat} \sim Pe$, explaining the result of Fig. 7b. This scaling highlights the effect of pore scale chaotic mixing on the enhancement of concentration gradients and of the resulting required mesh refinement.

On a physical level, the steady Eq. (14a) reflects mass conservation, i.e. the flux of a conserved species $\delta$ is the same through each cross-section $\Omega_{\perp}(z)$ of the domain transverse to the flow direction. To ensure that our results are valid, this fundamental property should be respected also numerically. Fig. 8 shows the average flux of $\delta$, denoted as $J_{\delta}$, for different iterations of mesh refinement for the three considered Péclet numbers. The flux $J_{\delta}$ is here relatively constant, fluctuating around $\bar{J}_{\delta} \approx -0.04$, which is different from 0 because the domain is not symmetric. The fluctuations in $J_{\delta}$ are of the order
<figure><img src="image_10.png" /><figcaption>Fig. 7. (a) Total number of nodes following mesh generation of each refinement iteration for flow through the bead pack. (b) Number of nodes $N_{sat}$ at which the refinement saturates for different Pe, indicating a scaling $N_{sat} \sim Pe$. The red dashed line is the estimated linear relation as presented in Eq. (G.5) where $N_{\perp} = 23$, $\lambda = 0.2$, and $d = 0.2$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Average flux of conserved species $\delta$ ($J_{\delta}$) at Pe equal to (a) $1.3 \cdot 10^1$, (b) $1.3 \cdot 10^2$, and (c) $1.3 \cdot 10^3$ for different iterations of mesh refinement in the case of flow through a bead pack.</figcaption></figure>

$< 0.005$, which must be compared to the total flux of velocity, which by construction is $J \approx 1$ (since the mean velocity $U = 1$). The relatively consistent values of the average $J_{\delta}$ for each value of Pe demonstrate that mass conservation is well respected. Furthermore, the fluctuations in $J_{\delta}$ decay with increasing refinement for all Pe values. This result suggests that the deviations from mass conservation originate from an imperfect input velocity field due to the fact that the solution method (FEM) only weakly satisfies the divergence-free condition (Eq. (B.2)). Altogether, these results establish confidence in our simulations.

To further investigate mixing-limited reactions inside the bead pack, we examine the distribution of the reaction product C for different refinement iterations within the bead pack for flow with $Pe = 1.3 \cdot 10^3$. Figs. 9a and 9b show the integrated mass $M_c$ (Eq. (26)) and reaction rate (Eq. (23)) of product C in each cross-section along the z-direction, respectively. We verify that refinement iterations increase the accuracy of the distribution of c and its reaction dynamics in the bead pack. In Supplementary Information (Movies S1–S6), we provide videos of the spatial distribution of $\delta$ and c in cross sections moving along the z-direction in the bead pack.

The simulations of this section provide insights into the behavior of solute mixing under conditions of chaotic flow through a bead pack. The adaptive mesh refinement method successfully captures the intricate patterns of mixing and reaction at the pore scale, demonstrating its

efficiency in handling complex geometries. These results advance our understanding of how chaotic flow regimes influence solute dispersion and reaction rates, providing a robust framework for future studies in granular porous media.

### 4.3. Application in natural rock: Berea sandstone sample

Finally, we demonstrate the applicability of the method to scenarios involving a natural rock with a more complex and rough pore space and show the ability of the method to predict and analyze solute transport in a real-world heterogeneous environment. We consider flow through a Berea sandstone sample with a porosity of 24% and with a dimensionless pore size $d = 0.15$, a value close to the size of the beads, $d = 0.2$, in the bead pack simulations. A subsample of the digitized rock was taken from the “11 sandstones” dataset available through the Digital Rocks Portal (Neumann et al., 2020), and meshed using CGAL (The CGAL Project, 2023). The initial tetrahedral mesh contained approximately 1.3 million nodes. Fig. 10a shows the steady-state velocity field obtained by numerically solving the Stokes equations (see Appendix B) in this domain. We added thin slabs of fluid at the inlet and outlet of the rock sample (see 10b) to make sure that the mesh consisted of a single connected entity. For the solute transport simulation, we chose $Pe = 100$, calculated using Eq. (21), in which
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Mass of reaction product C integrated over the cross-section at varying distances *z* from the inlet for flow through a bead pack with Pe = 1.3 · 10³ for different iterations. (b) Corresponding integrated reaction rates of *c* over the cross-section for different iterations of mesh refinement.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** (a) Steady-state velocity field **u** in a Berea sandstone sample obtained by solving Stokes equations. (b) Concentration of conserved species **δ** based on the velocity field **u** for Pe = 100. (c) The concentration *c* of product C is calculated assuming an instantaneous first-order bimolecular reaction (5), and computed based on the conservative solution of **δ**.</figcaption></figure>

*d* = 0.15, corresponding to the average pore size, and *U* = 1 is the mean normalized velocity. Fig. 10b shows the distribution of conserved species *δ* obtained using our numerical simulations. Fig. 10c shows the concentration *c* of product C assuming the first-order instantaneous reaction (5).

Analyzing the mass (Eq. (26)) and flux (Eq. (25)) of the different species within cross-sections perpendicular to the direction of flow is useful to understand the dynamics of reactive mixing through rock. Fig. 11a shows the integrated cross-sectional mass derived from the averaged concentration profiles for species A, B, and C. The concentration distributions were calculated using Eqs. (20a)–(20c), respectively. As expected, due to the reaction, we observe a general decreasing trend in reactant species concentrations *a* and *b*, while product concentration *c* increases. Fig. 11b shows the average cross-sectional flux of the species, including the conservative component *δ*. The fluxes of *a* and *b* decrease as they are consumed in the reaction, while accordingly the flux of *c* shows an increase. The flux for the conservative component *δ* remains constant, indicating that mass conservation is satisfied numerically.

The iterative mesh refinement process for the Berea sandstone required four iterations to satisfy the criterion Eq. (29), as shown in Fig. H.16a. It is important to note that a finer initial mesh reduces the number of iterations required to satisfy the mesh criteria Eq. (29). For example, in the case of uniform flow in Section 4.1, a higher number of iterations was required for the refinement process to saturate as the initial mesh was coarser. However, the final mesh size is relatively insensitive to the initial mesh size, as it is typically at least an order of magnitude higher.

The analysis of the parallel performance of our method, including the scaling of Berea sandstone simulations on an in-house computing

server, is provided in Appendix H. This analysis examines the solution time as a function of the number of CPU cores and evaluates the impact of mesh size on scalability. For the mixing and transport beyond the Representative Elementary Volume (REV) scale, we can estimate the number of nodes $N_{sat}$ required to reach the mesh refinement criterion in the entire domain as the function of Pe. The computational costs increase proportionally with $N_{sat}$, depending on the flow regime. The relation between $N_{sat}$ and Pe for a 2D and 3D uniform flow is estimated as Eqs. (F.9) and (F.10) respectively, and for a 3D chaotic flow is estimated as Eq. (G.5).

## 5. Conclusion

We have introduced a computational workflow for quantifying solute mixing and mixing-limited reactions in porous media, with a focus on numerically resolving transverse mixing and concentration gradients at the pore scale. The presented simulations advance knowledge in solute transport and mixing by validating the accuracy of the numerical method through comparisons with analytical solutions, capturing the complex dynamics of chaotic flow using adaptive mesh refinement, and demonstrating the practical application of the methodology to heterogeneous porous media. Dynamic refinement of the mesh in regions with sharp concentration gradients, ensures high resolution that is computationally efficient and feasible. In contrast to the time-dependent models where the number of required nodes evolves as $N \sim Pe^3$ (see above), we find $N \sim Pe$ with our method. The key elements of our method are to combine stabilization techniques and an iterative refinement procedure to obtain highly resolved concentration fields in steady state mixing fronts. Our method takes as input a mesh representing the
<figure><img src="image_14.png" /><figcaption>Fig. 11. Cross-sectional (a) total mass per length and (b) integrated flux of different species for flow through a Berea sandstone sample.</figcaption></figure>

discretized fluid domain and a steady-state velocity field. Through several numerical tests and examples, we have demonstrated that it constitutes an accurate and robust tool for investigating mixing and reaction phenomena at the pore scale, including in natural rock. We have thus illustrated and rationalized how the degree of required mesh refinement depends on the complexity and topology of the flow. In particular, the transition to chaotic flow in three dimensions leads to a significant enhancement of concentration gradients, which is captured by our mesh refinement technique. Our implementation of the methodology, based on the finite element method through the open-source FEniCS framework (Alnæs et al., 2015), is openly available from our Git repository (Shafabakhsh and Linga, 2024).

The presented methodology and accompanying code could become a useful tool for the porous media community in uncovering how microstructure influences effective mixing and reaction rates across different types of porous media. These highly resolved simulations open new perspectives to investigate how chaotic mixing dynamics observed in rock and granular media (Lester et al., 2016; Heyman et al., 2020, 2021; Aquino et al., 2023) influences fluid-fluid and fluid-solid reactions. To address this point, future study should consider larger system sizes and different injection scenarios to robustly extract asymptotic quantities.

Future work should include modeling the effects of finite reaction kinetics (Bandopadhyay et al., 2017), which will require numerically solving the nonlinear Eq. (14b). In this context, for fast, but finite reactions, the concentration distributions for infinite reaction kinetics (20) will provide good initial guesses for iterative nonlinear solvers (e.g. using Newton's method). Finally, our simulation methodology can be extended to arbitrarily complex chemical reaction networks, as long as we are assuming instantaneous reactions in the fluid phase. For example, an implementation of the reaction network considered by De Simoni et al. (2007b) is already available in the module post-process_crn of AdDiCTIF (Shafabakhsh and Linga, 2024). The influence of fluid mixing on more complex reactions (e.g., Steefel et al., 2013; Huppert and Neufeld, 2014; Soltanian et al., 2016; Kampman et al., 2014) and heterogeneous reactions occurring at fluid-solid interfaces (Roller, 2002; Liu and Mostaghimi, 2018) are important prospects for further investigation.

### CRediT authorship contribution statement

**Paiman Shafabakhsh:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Investigation. Tanguy Le

**Borgne:** Writing – review & editing, Supervision, Formal analysis, Conceptualization. **François Renard:** Writing – review & editing, Supervision, Project administration, Funding acquisition. **Gauté Linga:** Writing – review & editing, Writing – original draft, Supervision, Software, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

The meshes and velocity fields of the bead pack (Section 4.2) and the Berea sample (Section 4.3) are shared in a dedicated data repository (Shafabakhsh et al., 2024b). The simulation code is openly accessible at our Git repository (Shafabakhsh and Linga, 2024).

### Acknowledgments

The authors thank Joachim Mathiesen (University of Copenhagen) and Joris Heyman (University of Rennes) for stimulating discussions. The project received funding from the Akademiaavtalen between Equinor and the University of Oslo (project Modiflow). GL acknowledges support from the Research Council of Norway through grants number 325819 (M4) and 262644 (PoreLab Center of Excellence). FR acknowledges support from the FricFrac project funded by the Center for Advanced Study (CAS) at the Norwegian Academy of Science and Letters during the academic year 2023–2024. Simulations were performed on resources provided by Sigma2 – the National Infrastructure for High-Performance Computing and Data Storage in Norway (projects NN8108K, NN10091K).

## Appendix A. Analytical solutions for uniform flow

The advection-diffusion equation (1) for $\delta$ (i.e., Eq. (14a)) simplifies to

$$U \frac{\partial \delta}{\partial z} = D \nabla^2 \delta = D \left[ \frac{\partial^2 \delta}{\partial x^2} + \frac{\partial^2 \delta}{\partial z^2} \right], \quad (A.1)$$

where in the last equality, we used the fact that the system is invariant in the y direction through the boundary condition Eq. (30). Neglecting
diffusion along z (which is acceptable for advection-dominated flows) results in

$$ \frac{\partial \delta}{\partial z} = \alpha \frac{\partial^2 \delta}{\partial x^2}, \qquad (A.2) $$

where $\alpha = D/U$. Eq. (A.2) is easily identified as a simple one-dimensional diffusion equation, where z assumes the role of a time coordinate. The solution subject to Eq. (30) is given by

$$ \delta = \operatorname{erf} \left( \frac{x}{\sqrt{4\alpha z + \epsilon^2}} \right) \approx \operatorname{erf} \left( \frac{x}{\sqrt{4\alpha z}} \right), \qquad (A.3) $$

where the last approximation is valid for $\epsilon \to 0$.

We consider the conservative scalar dissipation rate defined by Eq. (22), and using Eq. (31), we approximate

$$ I \approx D \left| \frac{\partial \delta}{\partial x} \right|^2 = \frac{4D}{\pi(4Dz + \epsilon^2)} \exp \left( -\frac{2x^2}{4Dz + \epsilon^2} \right). \qquad (A.4) $$

This equation yields to the depth-averaged expression

$$ I_x(z) = D \langle |\nabla \delta|^2 \rangle_x \approx D \left\langle \left| \frac{\partial \delta}{\partial x} \right|^2 \right\rangle_x \qquad (A.5) $$

$$ = \frac{4D}{\pi(4Dz + \epsilon^2)} \sqrt{\frac{\pi}{2}(4Dz + \epsilon^2)} = \frac{2\sqrt{2D}}{\sqrt{\pi(4Dz + \epsilon^2)}}. \qquad (A.6) $$

## Appendix B. Stokes equations for creeping flow

The Stokes equations governing steady, viscosity-dominated incompressible fluid flow in the pore space are given by:

$$ \mu \nabla^2 \mathbf{u} = -\nabla p + \mathbf{f} \qquad (B.1) $$

$$ \nabla \cdot \mathbf{u} = 0, \qquad (B.2) $$

defined on a domain $\Omega$ (see Fig. 1). Here, $\mu$ is the dynamic viscosity of the fluid, $p$ is the fluid pressure, and $\mathbf{f} = f\hat{\mathbf{z}}$ is a body force (e.g., gravity). In our computations, we impose the boundary conditions

$$ p = 0 \qquad \text{for } \mathbf{x} \in \Gamma_{\text{in}} \cup \Gamma_{\text{out}} \qquad (B.3a) $$

$$ \mathbf{u} = 0 \qquad \text{for } \mathbf{x} \in \Gamma_{\text{wall}} \setminus \Gamma_{\text{box}}, \qquad (B.3b) $$

$$ \mathbf{u} \cdot \hat{\mathbf{n}} = 0, \hat{\mathbf{n}} \cdot (\nabla \mathbf{u} + \nabla \mathbf{u}^T) \times \hat{\mathbf{n}} = 0 \qquad \text{for } \mathbf{x} \in \Gamma_{\text{box}}. \qquad (B.3c) $$

Here, $\Gamma_{\text{box}}$ refers to the lateral boundary planes of the domain $\Omega$, where we through Eq. (B.3c) apply a free-slip condition (instead of a no-slip condition Eq. (B.3b)) to limit their effects on the flow (Sole-Mari et al., 2022). The velocity field is obtained by solving Eqs. (B.1)–(B.3).

## Appendix C. Progression of mesh refinement

Fig. C.12 shows an example of two-dimensional mesh refinement progression for the case of uniform flow (Section 4.1). After the fourth iteration, the refinement modifies only the mesh near the inlet ($x \approx z \approx 0$) due to the sharp gradient in concentration in that area.

To achieve the desired level of precision, we performed the iterative refinement process described in Section 3.3, until no elements were marked for refinement according to the selected criterion of Eq. (29). Fig. C.13a shows how the total number of mesh nodes depends on the refinement iteration for the uniform flow, comparing discretizations in two and three dimensions. For all four values of $\alpha$, the convergence criterion is fully reached after eight iterations. A longitudinal two-dimensional cross-section of the mesh throughout the iterative refinement process is shown in Fig. C.12 of Appendix C.

The refinement process can be strongly influenced by the choice of a finite smoothing scale $\epsilon$ of the inlet profile and we demonstrate that the dependence is qualitatively different in three dimensions than in two dimensions. With a very short transition length $\epsilon$ in the inlet condition Eq. (30), relatively more refinement is required near the inlet (as supported in Appendix C). The number of mesh nodes $N_{\text{sat}}$ at which the refinement process saturates in Fig. C.13a is plotted as a

function of $\alpha$ in Fig. C.13b. In particular, we find that $N_{\text{sat}}$ scales with $\alpha$ differently in two and three dimensions, as shown by the superimposed lines in Fig. C.13b. This result is in accordance to theoretical estimates Eqs. (F.9) and (F.10) derived in Appendix F based on the known approximate solution (31). The difference in the slope underscores that refinement in three dimensions needs special care compared with two dimensions, especially as the prefactor in three dimensions Eq. (F.10) depends on $\epsilon$, whereas it does not in two dimensions Eq. (F.9).

## Appendix D. Calculation of the reaction rate integrated over a cross-section

Here, we calculate the total rate $R_c$ integrated over a cross-section of product formation for the instantaneous reaction (5). In the steady state, the advection-diffusion-reaction Eq. (7c) can be written as

$$ R = \nabla \cdot \mathbf{j}_c, \qquad (D.1) $$

where $\mathbf{j}_c$ is defined by Eq. (25). We may decompose $\mathbf{j}_c = \hat{\mathbf{z}}j_{c,z} + \mathbf{j}_{c,\perp}$, where $j_{c,z} = \hat{\mathbf{z}} \cdot \mathbf{j}_c$ represents the flux component in the z-direction and $\mathbf{j}_{c,\perp}$ is the flux in the plane perpendicular to the z-axis. Now, using $\nabla \cdot \mathbf{j}_c = \partial_z j_{c,z} + \nabla_{\perp} \cdot \mathbf{j}_{c,\perp}$, where $\nabla_{\perp} = \hat{\mathbf{x}}\partial_x + \hat{\mathbf{y}}\partial_y$, we can determine the total reaction rate of the product C within each cross-section perpendicular to the flow:

$$ R_c = \int_{\Omega_{\perp}} R dxdy = \int_{\Omega_{\perp}} (\partial_z j_{c,z} + \nabla_{\perp} \cdot \mathbf{j}_{c,\perp}) dxdy \qquad (D.2) $$

$$ = \frac{dJ_c}{dz} + \oint_{\Gamma_{\perp}} \hat{\mathbf{n}}_{\perp} \cdot \mathbf{j}_{c,\perp} ds \qquad (D.3) $$

Here, $\Gamma_{\perp}$ is the boundary of the cross-sectional domain $\Omega_{\perp}(z)$, and $\hat{\mathbf{n}}_{\perp}$ is the outward normal vector to $\Omega_{\perp}$ in the plane perpendicular to the z-direction. The last term of the surface integral was converted to a line integral using the divergence theorem. Using Eq. (25), along $\Gamma_{\perp}$, we have $\hat{\mathbf{n}} \cdot \mathbf{j}_c = -D\hat{\mathbf{n}} \cdot \nabla c = -Dn_z \partial_z c - D\sqrt{1-n_z^2} \hat{\mathbf{n}}_{\perp} \cdot \nabla_{\perp} c = -Dn_z \partial_z c + \sqrt{1-n_z^2} \hat{\mathbf{n}}_{\perp} \cdot \mathbf{j}_{c,\perp} = 0$, that is, $\hat{\mathbf{n}}_{\perp} \cdot \mathbf{j}_{c,\perp} = Dn_z \partial_z c / \sqrt{1-n_z^2}$, where $n_z = \hat{\mathbf{n}} \cdot \hat{\mathbf{z}}$. The integral is then given by

$$ \oint_{\Gamma_{\perp}} \hat{\mathbf{n}}_{\perp} \cdot \mathbf{j}_{c,\perp} ds = D \oint_{\Gamma_{\perp}} \frac{n_z}{\sqrt{1-n_z^2}} \partial_z c \, ds. \qquad (D.4) $$

For the advection-dominated flows we are interested in, the last integral is not expected to contribute significantly, which we also verified numerically (see below). Therefore, the total reaction rate of product C in each perpendicular cross-section can be estimated by Eq. (23).

Fig. D.15 shows the comparison of the calculated reaction rates of product C, computed using Eq. (23) or Eq. (D.2) directly, for flow through a bead pack with $Pe = 1.3 \cdot 10^3$ at the sixth iteration. Both plots exhibit a consistent trend, indicating agreement in the overall behavior predicted by the two approaches. The reaction rate derived from Eq. (D.2) shows higher fluctuations, attributed to direct integration of the concentration gradient. However, the reaction rate computed from Eq. (23) is calculated from the derivative of the average flux at each interface derived from the interpolated data, resulting in a smoother profile.

## Appendix E. Mass and reaction rate of product C for uniform flow

Here, we present the analytical solution of the mass and the reaction rate of product formation for uniform flow. Following Eq. (26), the total mass $M_c$ of product C per cross-section normal to the z direction is:

$$ M_c(z) = \int_{-\infty}^{\infty} c(x, z) dx \qquad (E.1) $$

where $c(x, z)$ can be obtained from Eqs. (20c) and (31) as

$$ c(x, z) = \frac{1 - |\delta|}{2} = \frac{1}{2} \left[ 1 - \operatorname{erf} \left( \frac{|x|}{\sqrt{4\alpha z}} \right) \right]. \qquad (E.2) $$
<figure><img src="image_15.png" /><figcaption>Fig. C.12. Generated meshes during iterative refinement for uniform flow with $\alpha = 0.01$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. C.13. (a) Total number of nodes following mesh generation of each iteration for uniform flow in two and three dimensions. (b) Number of nodes $N_{sat}$ at which the refinement saturates for different $\alpha$. The red and blue curves are the best fit for two and three dimensions cases respectively based on estimates in Appendix F. The rightmost points do not follow the predicted scaling, as $\alpha$ and thus longitudinal gradients are too high for the analytical solution to apply.</figcaption></figure>

Exploiting the symmetry $c(x, z) = c(-x, z)$, we may consider only $x > 0$:

$$M_c(z) = 2 \int_{0}^{\infty} c(x, z) dx \qquad (E.3)$$

$$= \int_{0}^{\infty} \left[ 1 - \operatorname{erf} \left( \frac{x}{\sqrt{4\alpha z}} \right) \right] dx. \qquad (E.4)$$

By defining $\xi = x/\sqrt{4\alpha z}$ and substituting for Eq. (E.4), the analytical solution for the mass of C is:

$$M_c(z) = \sqrt{4\alpha z} \int_{0}^{\infty} (1 - \operatorname{erf}(\xi)) d\xi \qquad (E.5)$$

$$= \sqrt{\frac{4\alpha z}{\pi}} \qquad (E.6)$$

To obtain the reaction rate in uniform flow from Eq. (23), we first need to calculate the flux of product C along the z direction. Following Eq. (25), we have:

$$J_c = \int_{-\infty}^{\infty} j_{c,z} dx = \int_{-\infty}^{\infty} (u_z c - D \partial_z c) dx = U M_c - D \frac{dM_c}{dz}, \qquad (E.7)$$

where $U$ is the constant velocity ($u = U\hat{z}$). By substituting Eq. (E.6) for $M_c$, the flux can be estimated as

$$J_z \approx U \sqrt{\frac{4\alpha z}{\pi}} - D \sqrt{\frac{\alpha}{\pi z}}. \qquad (E.8)$$

The analytical solution of reaction rate for this case is:

$$R_c = \frac{dJ_c}{dz} = U \sqrt{\frac{\alpha}{\pi z}} \left(1 + \frac{\alpha}{2z}\right). \qquad (E.9)$$

## Appendix F. Estimating $N_{sat}$ for uniform flow

Based on the known analytical solution, we can estimate the number of nodes $N_{sat}$ required to reach the mesh refinement criterion in the entire domain. The analytical solution, including an initial smearing length scale $\epsilon$, is given by modifying Eq. (31):

$$\delta(x, z) = \operatorname{erf} \left( \frac{x}{\sqrt{4\alpha(z + z_0)}} \right). \qquad (F.1)$$
<figure><img src="image_17.png" /><figcaption>Fig. C.14. The same zoom-in at the (a) initial mesh and (b) the mesh at the final iteration of refinement for flow with Pe = 1.3 · 10³ in a bead pack.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. D.15. Comparison of reaction rates of product C computed from Eqs. (23) and (D.2) at varying distances z from the inlet for flow through a bead pack with Pe = 1.3 · 10³ for mesh refinement at the sixth iteration.</figcaption></figure>

Here, we have assumed that the interface is centered at $x_0 = 0$. Moreover, $z_0$ is related to $\epsilon$ by

$$ \delta(x, 0) = \operatorname{erf}\left(\frac{x}{\sqrt{4\alpha z_0}}\right) = \operatorname{erf}\left(\frac{x}{\epsilon}\right) \quad \Rightarrow \quad \epsilon = \sqrt{4\alpha z_0}. \qquad (F.2) $$

Assuming that gradients in $\delta$ are always contained in the domain $x \in [-L_x/2, L_x/2]$, we know that a fixed number of cells $N_x \approx 2/\Delta\delta$ is required across the interface to satisfy the criterion Eq. (29). (For our numerical examples, we have $N_x \approx 2/0.2 = 10$.) On the other hand, the width of the interface depends on z:

$$ w_x(z) \propto \sqrt{\alpha(z + z_0)}. \qquad (F.3) $$

### F.1. Scaling

We approximate the x-dependence of $\delta(x, z)$ by a linear function for $x \in (-w_x(z)/2, w_x(z)/2)$ and constant ($\pm 1$) otherwise. The local element density, at a given height position $z$, will be set by the local grid size limited by gradients in the x direction:

$$ \phi(z) = \frac{1}{\Delta x(z)} \approx \frac{N_x}{w_x} \sim \frac{N_x}{\sqrt{\alpha(z + z_0)}}. \qquad (F.4) $$

To estimate the prefactor of $w_x$, we may use the maximum slope of $\delta$, i.e. $|\partial_x \delta| = (\pi\alpha(z + z_0))^{-1/2}$, which yields

$$ w_x(z) \approx \frac{N_x \Delta \delta}{|\partial_x \delta|} \approx \sqrt{4\pi\alpha(z + z_0)} \equiv \beta \sqrt{\alpha(z + z_0)}. \qquad (F.5) $$

Now, since the grid size should be equal in all directions, we may estimate the number of nodes $N_{sat}$ by integrating over d spatial dimensions:

$$ N_{\text{sat}} \approx \int d^d \mathbf{x} \phi^d(z) \approx N_x^d \prod_{i=2}^d \int_0^{L_{x_i}} dx_i \int_{-w_x(z)/2}^{w_x(z)/2} dx w^{-d}(z) \qquad (F.6) $$

$$ = N_x^d \prod_{i=2}^d \int_0^{L_{x_i}} dx_i w^{-d+1}(z) = \frac{N_x^d L_{y,d}}{(\beta \sqrt{\alpha})^{d-1}} \int_0^{L_z} dz (z + z_0)^{-(d-1)/2} \qquad (F.7) $$

$$ = \frac{N_x^d L_{y,d}}{(\beta \sqrt{\alpha})^{d-1}} \int_{z_0}^{L_z + z_0} dz z^{-d(d-1)/2}, \qquad (F.8) $$

where $L_{y,d} = 1$ if $d = 2$ and $L_y$ if $d = 3$. The integral behaves differently depending on dimensionality $d$. We will handle these cases separately.

* In the case $d = 2$, we find

$$ N_{\text{sat}} = \frac{2N_x^2}{\beta \sqrt{\alpha}} (\sqrt{L_z + z_0} - \sqrt{z_0}) \approx \frac{2N_x^2}{\beta} \sqrt{\frac{L_z}{\alpha}}, \qquad (F.9) $$

where we have taken the limit $\epsilon \to 0$ in the last step.

* For $d = 3$, on the other hand, the latter limit is divergent:

$$ N_{\text{sat}} = \frac{N_x^3 L_y}{\beta^2 \alpha} \ln \left| 1 + \frac{L_z}{z_0} \right| \approx \frac{N_x^3 L_y}{\beta^2 \alpha} \ln \left| \frac{4\alpha L_z}{\epsilon^2} \right|. \qquad (F.10) $$

We observe that the scaling both with $\epsilon$ and $\alpha$ is very different depending on $d$.
<figure><img src="image_19.png" /><figcaption>Fig. H.16. (a) Total number of degrees of freedom (mesh nodes) as a function of refinement iteration for the case of flow through a Berea sandstone sample. (b) Scaling test: Simulation time as a function of the number of CPU cores for different iterations.</figcaption></figure>

## Appendix G. Estimating $N_{sat}$ for chaotic flow

We consider an interface (narrow mixing zone) in three dimensions, whose length can be described by $ℓ(z)$. The area of this interface is approximately

$$A \approx \int_{0}^{L_z} \ell(z) dz \qquad (G.1)$$

Across the interface, the number of nodes $N_⊥$ needed to satisfy the refinement criterion Eq. (29), will be independent of the value of Pe,

$$N_{\perp} \sim \frac{2}{\Delta_{\delta}}. \qquad (G.2)$$

However, the thickness of the interface is approximately given by the Batchelor scale,

$$s_B \sim d \sqrt{\frac{D}{U d \lambda}} \sim \frac{d}{\sqrt{\lambda}} \text{Pe}^{-1/2}, \qquad (G.3)$$

where $λ$ is the non-dimensional Lyapunov exponent and $U$ is a characteristic velocity. Neglecting over-refinement in the inlet region (which can be achieved by choosing a reasonable $ε$) and coalescence of the interface with itself, the local mesh size is thus everywhere approximately given by

$$h \sim \frac{s_B}{N_{\perp}}. \qquad (G.4)$$

The estimated total number of nodes is thus

$$N_{sat} \sim \frac{N_{\perp} A}{h^2} \sim \frac{N_{\perp}^3 A}{s_B^2} \sim \frac{N_{\perp}^3 \lambda A}{d^2} \text{Pe}. \qquad (G.5)$$

We expect that coalescence, increasingly at lower $Pe$, will contribute to lower (absolute) concentration gradients and effectively smaller $A$. Hence, the effective slope will be slightly higher than $N \sim Pe$.

## Appendix H. Computational time and CPU cores

To test the parallel performance of the method, we consider the scaling of the Berea sandstone simulations on an in-house computing server with 4 × 18 × 2 Intel(R) Xeon(R) Gold 6254 CPUs 3.10 GHz running Ubuntu 20.04 and FEniCS version 2019.2.0.13.dev0. Fig. H.16b represents the solution time as a function of the number of CPU cores employed for different iterations. We observe that the simulation time decreases for all iterations as we increase the number of CPU cores. However, the scalability of our code is limited for iterations three and

four, in which the mesh size is bigger, especially when we increase the number of CPU cores beyond 32. This result may reflect a reduction in mesh quality as many refinement iteration steps are performed.

## Appendix I. Supplementary material

Supplementary material includes the videos showing the spatial distribution of $δ$ and $c$ in cross sections moving along the $z$-direction in the bead pack and Berea sample (movies S1–S7).

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.advwatres.2024.104791.

## References

Abramowitz, M., Stegun, I.A., et al., 1968. Handbook of Mathematical Functions, Volume 10. Dover, New York.

Acharya, R.C., Valocchi, A.J., Werth, C.J., Willingham, T.W., 2007. Pore-scale simulation of dispersion and reaction along a transverse mixing zone in two-dimensional porous media. Water Resour. Res. 43.

Alhashmi, Z., Blunt, M., Bijeljic, B., 2015. Predictions of dynamic changes in reaction rates as a consequence of incomplete mixing using pore scale reactive transport modeling on images of porous media. J. Contam. Hydrol. 179, 171–181.

Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M.E., Wells, G.N., 2015. The fenics project version 1.5. Arch. Numer. Softw. 3.

d. Anna, P., Jimenez-Martinez, J., Tabuteau, H., Turuban, R., Le Borgne, T., Derrien, M., Méheust, Y., 2014. Mixing and reaction kinetics in porous media: An experimental pore scale quantification. Environ. Sci. Technol. 48, 508–516.

Aquino, T., Le Borgne, T., Heyman, J., 2023. Fluid-solid reaction in porous media as a chaotic restart process. Phys. Rev. Lett. 130, 264001.

Aramideh, S., Vlachos, P.P., Ardekani, A.M., 2018. Pore-scale statistics of flow and transport through porous media. Phys. Rev. E 98, 013104.

Baek, K.O., Seo, I.W., 2017. Estimation of the transverse dispersion coefficient for two-dimensional models of mixing in natural streams. J. Hydro-Environ. Res. 15, 67–74.

Bandopadhyay, A., Le Borgne, T., Méheust, Y., Dentz, M., 2017. Enhanced reaction kinetics and reactive mixing scale dynamics in mixing fronts under shear flow for arbitrary damköhler numbers. Adv. Water Resour. 100, 78–95.

Bauer, R.D., Maloszewski, P., Zhang, Y., Meckenstock, R.U., Griebler, C., 2008. Mixing-controlled biodegradation in a toluene plume—results from two-dimensional laboratory experiments. J. Contam. Hydrol. 96, 150–168.

Benekos, I.D., Cirkpa, O.A., Kitanidis, P.K., 2006. Experimental determination of transverse dispersivity in a helix and a cochlea. Water Resour. Res. 42.

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. Adv. Water Resour. 99, 15–37.

Benson, D.A., Bolster, D., 2016. Arbitrarily complex chemical reactions on particles. Water Resour. Res. 52, 9190–9200.
Benson, D.A., Pankavich, S., Bolster, D., 2019. On the separate treatment of mixing and spreading by the reactive-particle-tracking algorithm: An example of accurate upscaling of reactive poiseuille flow. Adv. Water Resour. 123, 40–53.

Bolster, D., Benson, D.A., Le Borgne, T., Dentz, M., 2010. Anomalous mixing and reaction induced by superdiffusive nonlocal transport. Phys. Rev. E 82, 021119.

Borgne, T.L., Ginn, T.R., Dentz, M., 2014. Impact of fluid deformation on mixing-induced chemical reactions in heterogeneous flows. Geophys. Res. Lett. 41, 7898–7906.

Boso, F., Battiato, I., 2013. Homogenizability conditions for multicomponent reactive transport. Adv. Water Resour. 62, 254–265.

Boso, F., Bellin, A., Dumbser, M., 2013. Numerical simulations of solute transport in highly heterogeneous formations: A comparison of alternative numerical schemes. Adv. Water Resour. 52, 178–189.

Brooks, A.N., Hughes, T.J., 1982. Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations. Comput. Methods Appl. Mech. Engrg. 32, 199–259.

Chiogna, G., Cirpka, O.A., Grathwohl, P., Rolle, M., 2011. Relevance of local compound-specific transverse dispersion for conservative and reactive mixing in heterogeneous porous media. Water Resour. Res. 47.

Cirpka, O.A., Chiogna, G., Rolle, M., Bellin, A., 2015. Transverse mixing in three-dimensional nonstationary anisotropic heterogeneous porous media. Water Resour. Res. 51, 241–260.

Cirpka, O.A., Frind, E.O., Helmig, R., 1999. Numerical simulation of biodegradation controlled by transverse mixing. J. Contam. Hydrol. 40, 159–182.

Cusini, M., Gielisse, R., Groot, H., van Kruisdijk, C., Hajibeygi, H., 2019. Incomplete mixing in porous media: Todd-longstaff upscaling approach versus a dynamic local grid refinement method. Comput. Geosci. 23, 373–397.

de Dreuzy, J.-R., Carrera, J., Dentz, M., Le Borgne, T., 2012. Time evolution of mixing in heterogeneous porous media. Water Resour. Res. 48.

De Simoni, M., Sanchez-Vila, X., Carrera, J., Saaltink, M., 2007a. A mixing ratios-based formulation for multicomponent reactive transport. Water Resour. Res. 43.

De Simoni, M., Sanchez-Vila, X., Carrera, J., Saaltink, M., 2007b. A mixing ratios-based formulation for multicomponent reactive transport. Water Resour. Res. 43.

Dell'Oca, A., Porta, G.M., Guadagnini, A., Riva, M., 2018. Space-time mesh adaptation for solute transport in randomly heterogeneous porous media. J. Contam. Hydrol. 212, 28–40.

Dentz, M., Hidalgo, J.J., Lester, D., 2023. Mixing in porous media: concepts and approaches across scales. Transp. Porous Media 146, 5–53.

Dentz, M., Le Borgne, T., Englert, A., Bijeljic, B., 2011. Mixing, spreading and reaction in heterogeneous media: A brief review. J. Contam. Hydrol. 120, 1–17.

Ding, D., Benson, D.A., Fernández-Garcia, D., Henri, C.V., Hyndman, D.W., Phanikumar, M.S., Bolster, D., 2017. Elimination of the reaction rate scale effect: Application of the lagrangian reactive particle-tracking method to simulate mixing-limited, field-scale biodegradation at the schoolcraft (mi, usa) site. Water Resour. Res. 53, 10411–10432.

Du, Z., Chen, J., Ke, S., 2023. Transverse mixing zone under dispersion in porous media: Effects of medium heterogeneity and fluid rheology. Phys. Fluids 35.

Eckert, D., Rolle, M., Cirpka, O.A., 2012. Numerical simulation of isotope fractionation in steady-state bioreactive transport controlled by transverse mixing. J. Contam. Hydrol. 140, 95–106.

Elman, H.C., Silvester, D.J., Wathen, A.J., 2014. Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics. Oxford University Press.

Engdahl, N.B., Ginn, T.R., Fogg, G.E., 2013. Scalar dissipation rates in non-conservative transport systems. J. Contam. Hydrol. 149, 46–60.

Engdahl, N.B., Schmidt, M.J., Benson, D.A., 2019. Accelerating and parallelizing lagrangian simulations of mixing-limited reactive transport. Water Resour. Res. 55, 3556–3566.

Feder, J., Flekkøy, E.G., Hansen, A., 2022. Physics of Flow in Porous Media. Cambridge University Press.

Fernández-Garcia, D., Sánchez-Vila, X., 2011. Optimal reconstruction of concentrations, gradients and reaction rates from particle distributions. J. Contam. Hydrol. 120, 99–114.

Geuzaine, Christophe, Remacle, Jean-Francois, 2024. Gmsh. URL: http://http://gmsh.info/.

Gramling, C.M., Harvey, C.F., Meigs, L.C., 2002. Reactive transport in porous media: A comparison of model prediction with laboratory visualization. Environ. Sci. Technol. 36, 2508–2514.

Guida, A., Nienow, A.W., Barigou, M., 2012. Lagrangian tools for the analysis of mixing in single-phase and multiphase flow systems. AIChE J. 58, 31–45.

Haynes, P.H., Vanneste, J., 2005. What controls the decay of passive scalars in smooth flows? Phys. Fluids 17.

Hejazi, S., Azaiez, J., 2013. Nonlinear simulation of transverse flow interactions with chemically driven convective mixing in porous media. Water Resour. Res. 49, 4607–4618.

Heyman, J., Borgne, T.L., Davy, P., Villermaux, E., 2023. Mixing as a correlated aggregation process. arXiv preprint arXiv:2309.14040.

Heyman, J., Lester, D., Le Borgne, T., 2021. Scalar signatures of chaotic mixing in porous media. Phys. Rev. Lett. 126, 034505.

Heyman, J., Lester, D.R., Turuban, R., Méheust, Y., Le Borgne, T., 2020. Stretching and folding sustain microscale chemical gradients in porous media. Proc. Natl. Acad. Sci. 117, 13359–13365.

Hidalgo, J.J., Dentz, M., Cabeza, Y., Carrera, J., 2015. Dissolution patterns and mixing dynamics in unstable reactive flow. Geophys. Res. Lett. 42, 6357–6364.

Hughes, T.J., 1979. A multidimensional upwind scheme with no crosswind diffusion. Finite Element Methods for Convect. Dominated Flows, AMD 34.

Huppert, H.E., Neufeld, J.A., 2014. The fluid mechanics of carbon dioxide sequestration. Annu. Rev. Fluid Mech. 46, 255–272.

Indelman, P., 2001. Steady-state source flow in heterogeneous porous media. Transp. Porous Media 45, 105–127.

Izumoto, S., Heyman, J., Huisman, J.A., De Vriendt, K., Soulaine, C., Gomez, F., Tabuteau, H., Méheust, Y., Le Borgne, T., 2023. Enhanced mixing and reaction in converging flows: Theory and pore-scale imaging. Water Resour. Res. 59, e2023WR034749.

Jha, B., Cueto-Felgueroso, L., Juanes, R., 2011. Quantifying mixing in viscously unstable porous media flows. Phys. Rev. E 84, 066312.

Jiménez-Martínez, J., Porter, M.L., Hyman, J.D., Carey, J.W., Viswanathan, H.S., 2016. Mixing in a three-phase system: Enhanced production of oil-wet reservoirs by CO₂ injection. Geophys. Res. Lett. 43, 196–205.

John, V., Knobloch, P., Novo, J., 2018. Finite elements for scalar convection-dominated equations and incompressible flow problems: a never ending story? Comput. Vis. Sci. 19, 47–63.

John, V., Novo, J., 2013. A robust SUPG norm a posteriori error estimator for stationary convection-diffusion equations. Comput. Methods Appl. Mech. Engrg. 255, 289–305.

Johnson, T.J., Locascio, L.E., 2002. Characterization and optimization of slanted well designs for microfluidic mixing under electroosmotic flow. Lab Chip 2, 135–140.

Jose, S.C., Cirpka, O.A., 2004. Measurement of mixing-controlled reactive transport in homogeneous porous media and its prediction from conservative tracer test data. Environ. Sci. Technol. 38, 2089–2096.

Kampman, N., Bickle, M., Wigley, M., Dubacq, B., 2014. Fluid flow and CO₂-fluid-mineral interactions during CO₂-storage in sedimentary basins. Chem. Geol. 369, 22–50.

Katz, G.E., Berkowitz, B., Guadagnini, A., Saaltink, M.W., 2011. Experimental and modeling investigation of multicomponent reactive transport in porous media. J. Contam. Hydrol. 120, 27–44.

Kim, M.C., 2019. Effect of the irreversible A + B → C reaction on the onset and the growth of the buoyancy-driven instability in a porous medium: Asymptotic, linear, and nonlinear stability analyses. Phys. Rev. Fluids 4, 073901.

Kitanidis, P.K., McCarty, P.L., 2012. Delivery and Mixing in the Subsurface: Processes and Design Principles for in Situ Remediation, Volume 4. Springer Science & Business Media.

Klimenko, A., 2009. Lagrangian particles with mixing, ii. sparse-lagrangian methods in application for turbulent reacting flows. Phys. Fluids 21.

Kumar, A., Jaiswal, D.K., Kumar, N., 2009. Analytical solutions of one-dimensional advection-diffusion equation with variable coefficients in a finite domain. J. Earth Syst. Sci. 118, 539–549.

Le Borgne, T., Dentz, M., Bolster, D., Carrera, J., De Dreuzy, J.-R., Davy, P., 2010. Non-fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media. Adv. Water Resour. 33, 1468–1475.

Lester, D.R., Dentz, M., Le Borgne, T., 2016. Chaotic mixing in three-dimensional porous media. J. Fluid Mech. 803, 144–174.

Lester, D., Dentz, M., LeBorgne, T., 2014. Topological mixing in three-dimensional porous media. In: Proceedings of the 19th Australasian Fluid Mechanics Conference 2014 Australasian Fluid Mechanics Society. AFMS 2014, Australasian Fluid Mechanics Society (AFMS), pp. 1–4.

Li, L., Peters, C.A., Celia, M.A., 2006. Upscaling geochemical reaction rates using pore-scale network modeling. Adv. Water Resour. 29, 1351–1370.

Liedl, R., Valocchi, A.J., Dietrich, P., Grathwohl, P., 2005. Finiteness of steady state plumes. Water Resour. Res. 41.

Liu, M., Mostaghimi, P., 2018. Numerical simulation of fluid-fluid-solid reactions in porous media. Int. J. Heat Mass Transfer 120, 194–201.

Luo, J., Dentz, M., Carrera, J., Kitanidis, P., 2008. Effective reaction parameters for mixing controlled reactions in heterogeneous media. Water Resour. Res. 44.

Mansell, R., Ma, L., Ahuja, L., Bloom, S., 2002. Adaptive grid refinement in numerical models for water flow and chemical transport in soil: a review. Vadose Zone J. 1, 222–238.

Martínez-Ruiz, D., Meunier, P., Favier, B., Duchemin, L., Villermaux, E., 2018. The diffusive sheet method for scalar mixing. J. Fluid Mech. 837, 230–257.

Minier, J.-P., Chibbaro, S., Pope, S.B., 2014. Guidelines for the formulation of lagrangian stochastic models for particle simulations of single-phase and dispersed two-phase turbulent flows. Phys. Fluids 26.

Molins, S., Knabner, P., 2019. Multiscale approaches in reactive transport modeling. Rev. Mineral. Geochem. 85, 27–48.

Nambi, I.M., Werth, C.J., Sanford, R.A., Valocchi, A.J., 2003. Pore-scale analysis of anaerobic halorespiring bacterial growth along the transverse mixing zone of an etched silicon pore network. Environ. Sci. Technol. 37, 5617–5624.

Neumann, R., Andreeta, M., Lucas-Oliveira, E., 2020. 11 Sandstones: raw, filtered and segmented data. http://dx.doi.org/10.17612/f4h1-w124, http://www.digitalrocksportal.org/projects/317.
Nilsson, J., Gerritsen, M., Younis, R., 2005. A novel adaptive anisotropic grid framework for efficient reservoir simulation. In: SPE Reservoir Simulation Conference?. SPE, pp. SPE-93243.

Noetinger, B., Roubinet, D., Russian, A., Le Borgne, T., Delay, F., Dentz, M., De Dreuzy, J.-R., Gouze, P., 2016. Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale. Transp. Porous Media 115, 345-385.

Oostrom, M., Mehmani, Y., Romero-Gomez, P., Tang, Y., Liu, H., Yoon, H., Kang, Q., Joekar-Niasar, V., Balhoff, M., Dewers, T., et al., 2016. Pore-scale and continuum simulations of solute transport micromodel benchmark experiments. Comput. Geosci. 20, 857-879.

Pau, G.S., Almgren, A.S., Bell, J.B., Lijewski, M.J., 2009. A parallel second-order adaptive mesh algorithm for incompressible flow in porous media. Phil. Trans. R. Soc. A 367, 4633-4654.

Perez, L.J., Hidalgo, J.J., Dentz, M., 2019. Upscaling of mixing-limited bimolecular chemical reactions in poiseuille flow. Water Resour. Res. 55, 249-269.

Plaza, A., Carey, G., 1996. A New Refinement Algorithm for Tetrahedral Grids Based on Skeleton. Texas Institute for Computational and Applied Mathematics, University of Texas at Austin.

Porta, G.M., Ceriotti, G., Thovert, J.-F., 2016. Comparative assessment of continuum-scale models of bimolecular reactive transport in porous media under pre-asymptotic conditions. J. Contam. Hydrol. 185, 1-13.

Porta, G.M., Chaynikov, S., Thovert, J.-F., Riva, M., Guadagnini, A., Adler, P.M., 2013. Numerical investigation of pore and continuum scale formulations of bimolecular reactive transport in porous media. Adv. Water Resour. 62, 243-253.

Rahman, M.A., Jose, S.C., Nowak, W., Cirpka, O.A., 2005. Experiments on vertical transverse mixing in a large-scale heterogeneous model aquifer. J. Contam. Hydrol. 80, 130-148.

Raje, D.S., Kapoor, V., 2000. Experimental study of bimolecular reaction kinetics in porous media. Environ. Sci. Technol. 34, 1234-1239.

Rolle, M., Eberhardt, C., Chiogna, G., Cirpka, O.A., Grathwohl, P., 2009. Enhancement of dilution and transverse reactive mixing in porous media: Experiments and model-based interpretation. J. Contam. Hydrol. 110, 130-142.

Rolle, M., Le Borgne, T., 2019. Mixing and reactive fronts in the subsurface. Rev. Mineral. Geochem. 85, 111-142.

Roller, P.S., 2002. The physical and chemical relatics in fluid phase heterogeneous reaction. J. Phys. Chem. 39, 221-238.

Roos, H.-G., 2008. Robust Numerical Methods for Singularly Perturbed Differential Equations. Springer.

Salamon, P., Fernández-Garcia, D., Gómez-Hernández, J.J., 2006. A review and numerical assessment of the random walk particle tracking method. J. Contam. Hydrol. 87, 277-305.

Sale, T.C., McWhorter, D.B., 2001. Steady state mass transfer from single-component dense nonaqueous phase liquids in uniform flow fields. Water Resour. Res. 37, 393-404.

Sanquer, H., Heyman, J., Hanna, K., Le Borgne, T., 2024. Microscale chaotic mixing as a driver for chemical reactions in porous media. Environ. Sci. Technol. 58, 8899-8908.

Shafabakhsh, P., Cordonnier, B., Pluymakers, A., Le Borgne, T., Mathiesen, J., Linga, G., Hu, Y., Kaestner, A., Renard, F., 2024a. 4D neutron imaging of solute transport and fluid flow in sandstone before and after mineral precipitation. Water Resour. Res. 60, e2023WR036293.

Shafabakhsh, P., Le Borgne, T., Renard, F., Linga, G., 2024b. Data repository. http://dx.doi.org/10.5281/zenodo.10908276, (Accessed 27 March 2024).

Shafabakhsh, P., Linga, G., 2024. AdDiCTIF: Git repository for simulation code. http://github.com/gautelinga/addictif. (Accessed 27 March 2024).

Siena, M., Hyman, J., Riva, M., Guadagnini, A., Winter, C.L., Smolarkiewicz, P., Gouze, P., Sadhukhan, S., Inzoli, F., Guédon, G., et al., 2015. Direct numerical simulation of fully saturated flow in natural porous media at the pore scale: a comparison of three computational systems. Comput. Geosci. 19, 423-437.

Smilauer, V., Angelidakis, V., Catalano, E., Caulk, R., Chareyre, B., Chevremont, W., Dorofeenko, S., Duriez, J., Dyck, N., Elias, J., et al., 2023. Yade documentation. arXiv preprint arXiv:2301.00611.

Sole-Mari, G., Bolster, D., Fernandez-Garcia, D., 2022. A closer look: High-resolution pore-scale simulations of solute transport and mixing through porous media columns. Transp. Porous Media 1-27.

Soltanian, M.R., Amooie, M.A., Dai, Z., Cole, D., Moortgat, J., 2016. Critical dynamics of gravito-convective mixing in geological carbon sequestration. Sci. Rep. 6, 35921.

Souzy, M., Lhuissier, H., Méheust, Y., Le Borgne, T., Metzger, B., 2020. Velocity distributions, dispersion and stretching in three-dimensional porous media. J. Fluid Mech. 891, A16.

Steefel, C.I., Molins, S., Trebotich, D., 2013. Pore scale processes associated with subsurface CO₂ injection and sequestration. Rev. Mineral. Geochem. 77, 259-303.

Sternagel, A., Loritz, R., Klaus, J., Berkowitz, B., Zehe, E., 2021. Simulation of reactive solute transport in the critical zone: a lagrangian model for transient flow and preferential transport. Hydrol. Earth Syst. Sci. 25, 1483-1508.

The CGAL Project, 2023. CGAL User and Reference Manual, 5.6 ed. CGAL Editorial Board, URL: https://doc.cgal.org/5.6/Manual/packages.html.

Turuban, R., Lester, D.R., Heyman, J., Le Borgne, T., Méheust, Y., 2019. Chaotic mixing in crystalline granular media. J. Fluid Mech. 871, 562-594.

Valocchi, A.J., Bolster, D., Werth, C.J., 2019. Mixing-limited reactions in porous media. Transp. Porous Media 130, 157-182.

Voronov, R.S., VanGordon, S.B., Sikavitsas, V.I., Papavassiliou, D.V., 2011. Efficient lagrangian scalar tracking method for reactive local mass transport simulation through porous media. Internat. J. Numer. Methods Fluids 67, 501-517.

Willingham, T.W., Werth, C.J., Valocchi, A.J., 2008. Evaluation of the effects of porous media structure on mixing-controlled reactions using pore-scale modeling and micromodel experiments. Environ. Sci. Technol. 42, 3185-3193.