# Reactive transport in porous media: Pore-network model approach compared to pore-scale model

Clément Varloteaux, Minh Tan Vu, Samir Bekri, Pierre M. Adler

► To cite this version:

Clément Varloteaux, Minh Tan Vu, Samir Bekri, Pierre M. Adler. Reactive transport in porous media: Pore-network model approach compared to pore-scale model. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2013, 87 (2), 10.1103/PhysRevE.87.023010. [https://doi.org/10.1103/PhysRevE.87.023010](https://doi.org/10.1103/PhysRevE.87.023010)

HAL Id: hal-01196671

https://hal.science/hal-01196671

Submitted on 10 Jul 2020

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
Reactive transport in porous media: Pore-network model approach compared to pore-scale model

Clément Varloteaux,¹,* Minh Tan Vu,² Samir Békri,¹,† and Pierre M. Adler²,‡

¹IFP Energies nouvelles, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France

²Sisyphe, UPMC, 4 place Jussieu, 75252 Paris, France

(Received 1 October 2012; published 15 February 2013)

Accurate determination of three macroscopic parameters governing reactive transport in porous media, namely, the apparent solute velocity, the dispersion, and the apparent reaction rate, is of key importance for predicting solute migration through reservoir aquifers. Two methods are proposed to calculate these parameters as functions of the Péclet and the Péclet-Dahmköhler numbers. In the first method called the pore-scale model (PSM), the porous medium is discretized by the level set method; the Stokes and convection-diffusion equations with reaction at the wall are solved by a finite-difference scheme. In the second method, called the pore-network model (PNM), the void space of the porous medium is represented by an idealized geometry of pore bodies joined by pore throats; the flow field is computed by solving Kirchhoff's laws and transport calculations are performed in the asymptotic regime where the solute concentration undergoes an exponential evolution with time. Two synthetic geometries of porous media are addressed by using both numerical codes. The first geometry is constructed in order to validate the hypotheses implemented in PNM. PSM is also used for a better understanding of the various reaction patterns observed in the asymptotic regime. Despite the PNM approximations, a very good agreement between the models is obtained, which shows that PNM is an accurate description of reactive transport. PNM, which can address much larger pore volumes than PSM, is used to evaluate the influence of the concentration distribution on macroscopic properties of a large irregular network reconstructed from microtomography images. The role of the dimensionless numbers and of the location and size of the largest pore bodies is highlighted.

DOI: 10.1103/PhysRevE.87.023010

PACS number(s): 47.56.+r, 47.11.—j, 47.15.G—

I. INTRODUCTION

When a reactive fluid is injected into a porous medium, chemical reactions can modify the petrophysical properties of the porous media, such as the porosity $\varepsilon$ and the permeability $K$. This physical phenomenon consists of two coupled major processes, namely, solute transport and reaction at the fluid-solid interface, which depend on the initial conditions, on the geometry, and also on the nature of the components. An important application of such studies is carbon dioxide storage in saline aquifers where dissolution of $CO_2$ in brine may cause its acidification and thus mineral dissolution [1,2].

The main objective of this paper is to determine the distribution of the solute by two different numerical methods and to compare the results.

In the pore-scale model (PSM) developed by Bekri et al. [3], the porous medium is represented by void and solid voxels (voxel method). Local equations governing the solute concentration are solved by the finite-difference method. Then, the evolution of rock-fluid interfaces is calculated and the porosity and permeability modifications are determined. In order to accurately simulate the complex surface motions, the voxel method is replaced by the level set method (LSM) [4]. The advantage of the LSM is that it can deal with curves and surfaces on a fixed Cartesian grid without having to parametrize these objects. Also, the LSM can follow shapes that change topology, for instance, when it splits into two or, reversely, develops holes [5].

The PSM combined with the LSM is an accurate method, but it is time consuming and only limited pore volumes can be addressed. An alternative method is the pore-network model (PNM) which allows one to study reactive transport phenomena in much larger pore volumes with the same computational resources. This approach is based on a simplified microstructure of the porous medium, which is schematized by pore bodies connected by pore throats [6–11].

The PNM is versatile and can account for various phenomena occurring on the pore scale. It was originally developed by Fatt [12] to calculate multiphase flow properties of porous media. Over the last decades, it has been extensively used to simulate basic phenomena such as capillarity and multiphase flow through porous media [8,13–16]. This approach was extended to study pore evolution and changes in petrophysical properties due to particle capture [17], asphalt precipitation [18], deposition and dissolution in diatomite [19], and filtration combustion [20]. Recently, adsorption and reaction processes were tentatively integrated into the PNM. Raoof et al. [21] quantified the effective kinetics of adsorption processes whereas Li et al. [22] and Kim [23] concentrated their research on effective reaction rates in porous media using the PNM and its possible implementation on the reservoir scale. Algive et al. [24,25] proposed the PNM approach to study mineral dissolution and precipitation caused by $CO_2$ sequestration.

In this paper, the PSM is used to validate the PNM for reactive transport phenomena since unexpected phenomena were observed in the PNM at the beginning of this work. Sections II and III describe the PSM and the PNM, respectively. In Sec. IV, the two models are compared on synthetic cases. Then, PNM hypotheses are validated and the accuracy of the model is assessed by PSM calculations; various reaction regimes which occur in a case study are illustrated and explained. Section V

<footnote>*clement.varloteaux@ifpen.fr</footnote>

<footnote>†samir.bekri@ifpen.fr</footnote>

<footnote>‡pierre.adler@upmc.fr</footnote>
is devoted to the application of the PNM to a large irregular network derived from microtomography measurements.

## II. PORE-SCALE MODEL

The pore-scale model is based on the resolution of the Stokes equation and of the convection-diffusion equation supplemented by conditions on the deposition or dissolution flux at the walls. These equations are discretized by a finite-difference scheme which is detailed by Bekri et al. [3].

The LSM is used to smooth out the fluid-rock interfaces. The most important improvement provided by the LSM is that the boundary conditions are written at the level set surface and not at the discretized surface by voxels. Furthermore, the flux at the interface is precisely computed by using the unit vector which is normal to the smooth surface.

### A. Governing equations

Consider reactive transport in a porous medium $\Omega$ which consists of a fluid phase $\Omega_F$ and a solid phase $\Omega_S$ separated by an interface $\Gamma$.

When the flow is steady and when inertial effects are negligible, the fluid motion is governed by the Stokes equation. Hence,

$$ \mu \nabla^2 \mathbf{v} + \nabla p = 0 \quad \text{in } \Omega_F, \quad (1a) $$

$$ \nabla \cdot \mathbf{v} = 0 \quad \text{in } \Omega_F, \quad (1b) $$

where $\mu$ is the viscosity of the fluid, which is assumed to be constant, $p$ is the pressure, and $\mathbf{v}$ is the fluid velocity.

The no-slip condition should be satisfied at the fluid-solid interface $\Gamma$:

$$ \mathbf{v} = 0 \quad \text{on } \Gamma. \quad (1c) $$

The solute flux $\mathbf{J}$ can be written as

$$ \mathbf{J} = c\mathbf{v} - D\nabla c, \quad (2a) $$

where $D$ is the solute molecular diffusion and $c$ is the solute concentration. $D$ is assumed to be constant.

When there is no bulk chemical reaction, $c$ obeys the local convection-diffusion equation

$$ \frac{\partial c}{\partial t} + \nabla \cdot (c\mathbf{v} - D\nabla c) = 0 \quad \text{in } \Omega_F. \quad (2b) $$

The boundary condition for $c$ at the wall $\Gamma$ is assumed to be a first-order surface reaction:

$$ \mathbf{n} \cdot \mathbf{J} = \kappa(c - \tilde{c}) \quad \text{on } \Gamma, \quad (2c) $$

where $\kappa$ is the local reaction rate constant and $\tilde{c}$ is the equilibrium concentration of the solute.

This reaction causes a displacement $W$ normal to the wall, which is proportional to the solute flux at the wall [3],

$$ \frac{\partial W}{\partial t} = -K_c \rho_F \kappa (c - \tilde{c}) \quad \text{on } \Gamma, \quad (3) $$

where $\rho_F$ is the fluid density and $K_c$ is the stoichiometric coefficient of the reaction. Of course, the velocity of this displacement is assumed to be very small with respect to the fluid velocity.

### B. Dimensionless formulation

In order to derive the parameters which control the problem, the previous equations can be made dimensionless by introducing a characteristic length scale $l_c$ for the porous medium and a characteristic velocity chosen as the interstitial velocity $\langle v \rangle$. Similarly, the average concentration $\langle c \rangle$ is defined as the concentration scale. Thus, a new system of dimensionless variables indicated by primes can be defined:

$$ \begin{gathered} \nabla' = l_c \nabla, \quad \mathbf{v}' = \frac{\mathbf{v}}{\langle v \rangle}, \quad p' = \frac{pl_c}{\mu \langle v \rangle}, \quad c' = \frac{c - \tilde{c}}{\langle c \rangle - \tilde{c}}, \quad (4) \\ t' = \frac{t}{T} \quad \text{with} \quad T = \frac{l_c^2}{K_c \rho D (\langle c \rangle - \tilde{c})}. \end{gathered} $$

The dimensionless equations are

$$ \nabla'^2 \mathbf{v}' - \nabla' p' = 0 \quad \text{in } \Omega_F, \quad (5a) $$

$$ \nabla' \cdot \mathbf{v}' = 0 \quad \text{in } \Omega_F, \quad (5b) $$

$$ \mathbf{v}' = 0 \quad \text{on } \Gamma, \quad (5c) $$

$$ \text{Pe} \nabla' c' \cdot \mathbf{v}' - \nabla'^2 c' = - \frac{L^2}{DT} \frac{\partial c'}{\partial t} \quad \text{in } \Omega_F, \quad (5d) $$

$$ \mathbf{n} \cdot \nabla' c' = -\text{PeDa} c' \quad \text{on } \Gamma, \quad (5e) $$

$$ \frac{\partial W'}{\partial t'} = \text{PeDa} c' \quad \text{on } \Gamma, \quad (5f) $$

where the dimensionless Péclet and Dahmköhler numbers are defined as

$$ \text{Pe} = \frac{\langle v \rangle l_c}{D}, \quad \text{Da} = \frac{\kappa}{\langle v \rangle}. \quad (6) $$

Pe compares convection and diffusion while Da compares the speed of the chemical reaction and the fluid velocity. The product of these two numbers, PeDa, is often used; it compares reaction to diffusion characteristic times.

In addition, the system is assumed to be not very far from chemical equilibrium and that the rate of deformation of the solid surface is very slow; hence, the velocity field in the fluid can be determined at any time by solving (5a)–(5c). For a given geometry and velocity field, the dimensionless equations governing $c'$ possess a solution of the form $e^{-\gamma t} c'(\mathbf{x}')$. We shall assume that the time required to reach the asymptotic regime is small compared to the wall evolution characteristic time. Consequently, the geometrical changes mainly occur during the asymptotic regime and the asymptotic concentration field is used to determine the wall evolution rate.

### C. Level set method

The solid-liquid interface is tracked by means of the LSM. In this method, the real surface is defined by a distance function based on the usual fixed Cartesian grid. The interface is represented by a triangulated surface at the zero level of this distance function. The numerical codes using this method solve transport and flow fields more accurately and topology changes are effectively handled [5].

The LSM represents the interface $\Gamma(\mathbf{x})$ as the zero-level contour of a function $\phi(\mathbf{x}, t)$ [26]:

$$ \Gamma = \{\mathbf{x} | \phi(\mathbf{x}, t) = 0\}. \quad (7) $$
Moreover, the level set function $\phi(\mathbf{x},t)$ satisfies the properties that $\phi > 0$ for phase 1 and $\phi < 0$ for phase 2. In practice, $\phi(\mathbf{x},t)$ is a signed distance function.

The chain derivation rule applied to $\phi(\mathbf{x},t)$ yields

$$ \frac{\partial \phi}{\partial t} + \nabla \phi \cdot \frac{d\mathbf{x}'}{dt} = 0. \qquad (8) $$

Let $\mathbf{n}$ be the normal to $\Gamma$ pointing outward of the solid phase:

$$ \mathbf{n} = \frac{\nabla \phi}{\|\nabla \phi\|}. \qquad (9) $$

The velocity of the interface is defined along $\mathbf{n}$ as

$$ V_i = \frac{d\mathbf{x}'}{dt} \cdot \mathbf{n}. \qquad (10) $$

This propagation velocity $V_i$ is related to (5f) by

$$ V_i = \frac{\partial W'}{\partial t} = \text{PeDa } c' \quad \text{on } \Gamma. \qquad (11) $$

Finally, the evolution of the level set function $\phi$ is governed by

$$ \frac{\partial \phi}{\partial t} + V_i \|\nabla \phi\| = 0 \qquad (12) $$

for a given initial geometry $\phi(\mathbf{x},t = 0)$.

### D. Algorithm description

During the simulation, the coupled Stokes (5a) and convection-diffusion (5d) equations are solved by the same algorithm as used by Bekri et al. [3] which comprises five steps. (i) The velocity field $v^n$ and the concentration field $c^n$ are calculated for the current interface $\phi^n$. (ii) The concentration at the interface is extrapolated from the field $c^n$ to calculate the interface propagation velocity $F$ (11). (iii) The new interface $\phi^{n+1}$ is determined at time $t^{n+1}$ by using the interface velocity and (12). (iv) The new interface is updated and the medium is visualized. (v) This process is repeated until the end condition is verified.

This algorithm is schematized in Fig. 1.

## III. PORE-NETWORK MODEL

The PNM describes the flow and the transport on the pore scale. It can address larger pore volumes than can the PSM with the same computational resources. This section describes the porous medium representation and the resolution of flow and concentration in the asymptotic regime defined in Sec. II B; therefore, only long-term phenomena are studied.

### A. Geometry

The PNM is based on a simplified representation of the void space, which is approximated by a network of bonds (pore throats) and nodes (pore bodies) with an idealized geometry; the pore bodies are spherical while the pore throats are cylindrical channels with a circular, square, or triangular cross section. The distinction between pore bodies and pore throats and their simplified geometry makes complex problems easier to solve by using analytical or semianalytical solutions [27].

<figure><img src="image_2.png" /><figcaption>FIG. 1. General scheme of the reactive transport resolution by the PSM combined with the LSM.</figcaption></figure>

The pore network can be a regular or an irregular three-dimensional lattice structure (Fig. 2). In Fig. 2(a), the pore space is defined on a cubic lattice where each pore body is assumed to be connected to six pore throats; therefore, the coordination number is equal to 6. The ratio between the pore-body and the pore-throat diameters (aspect ratio) is constant. The pore-throat diameters are randomly generated according to a given probability density function. Of course, the coordination number and the aspect ratio can be variable.

In order to construct a representative pore network of a porous medium, the probability density function has to be chosen in order to reproduce some petrophysical parameters

<figure><img src="image_3.png" /><figcaption>FIG. 2. Pore-network models (a) reconstructed with a regular lattice in order to reproduce petrophysical properties of real porous media [28] and (b) extracted from microtomography measurements [29].</figcaption></figure>
such as porosity and permeability. In addition, Bekri and Vizika [28] recommend that the formation factor and the capillary pressure curve be equal to those of the considered porous medium since they are very sensitive to its structure. The choice of a compatible pore-throat size distribution is a key parameter for the construction of a representative pore network using this cubic structure.

An alternative to the regular lattice pore network has been recently developed in order to get closer to the real medium geometry [7,29–31]. This method made important progress because of synchrotron computed microtomography, which generates three-dimensional (3D) data sets on the micrometer scale.

The first step of this method is to measure the exact 3D pore space of the porous medium. Then, a 3D image in gray levels is reconstructed using x-ray microtomography. A threshold is chosen to distinguish between pores and rock. Then, the skeleton of the pore space is computed by a hybrid algorithm which combines thinning and a distance map such as the one derived by Thovert et al. [32]. Additionally, the pore space is partitioned into pore bodies and pore throats according to the conceptual description of the pore-network model [Fig. 2(b)]. Finally, geometrical parameters are extracted from the 3D pore-space images [see [29], for more details].

## B. Flow field

The fluid flow is governed by the Stokes equations (1). The fluid velocity in the capillary tubes which compose the pore network is given by the Poiseuille parabolic profile

$$v(\rho) = 2v_z \left(1 - \frac{\rho^2}{r^2}\right)z \quad \text{with} \quad v_z = -\frac{r^2}{8\mu} \frac{\partial P}{\partial z}, \qquad (13)$$

where $r$ and $\rho$ are the radius of the tube and the radial coordinate, respectively, and $z$ is the unit vector parallel to the tube axis. The pressure drop between two neighbor pore bodies, $\Delta P$, is related to the flow rate $Q$ passing through the capillary tube by

$$Q = v_z S = \frac{\pi r^4 \Delta P}{8\mu L_z}, \qquad (14)$$

where $L_z$ is the length of the capillary tube.

The pore-body flow field is conceptually more difficult to define analytically. For the sake of simplicity, pressure is assumed to be constant within a pore body; in other words, the velocity in pore bodies is assumed to be zero and therefore negligible compared to the velocity in the pore throats.

Then, a mass balance is performed over each pore body and a linear system for the pressures can be derived,

$$G \cdot P = b, \qquad (15)$$

where $G$ is the conductivity matrix, which only depends on the geometric properties of the network, $P$ is the unknown pressure vector, and $b$ is a vector related to the external boundary conditions.

By using a linear solver such as a conjugate gradient technique [33], the entire pressure field can be evaluated for an imposed pressure drop over the network. Thus, Eqs. (13) provide velocity in every pore throat.

When the pressure and velocity fields are known, the permeability $K$ of the porous medium can be computed as

$$K = \frac{Q_{\text{tot}} \mu L_{\text{tot}}}{S_{\text{tot}} \Delta P}, \qquad (16)$$

where $Q_{\text{tot}}$ is the total flow rate passing through the cross section $S_{\text{tot}}$ of the porous medium; $L_{\text{tot}}$ and $\Delta P$ are the total length and the pressure drop over the network.

## C. Concentration field

A comparable approach to flow is used to compute the concentration field within a pore network. An analytical solution of the local problem is provided for the simplified geometry used in the PNM. The flux at the pore body–pore throat interface and the concentration are related by an analytical solution. Mass balance yields a nonlinear system solved by an optimization algorithm which provides the concentration distribution within the pore network. This section presents the analytical solutions for the pore bodies and the pore throats. Then, the concentration field within the whole pore network is determined.

### 1. Pore throats

The resolution of the reactive transport in a pore throat is divided into two steps. First, the transverse profile is calculated by assuming that the flow and transport transverse profiles are established. Second, a macroscopic equation governing the average concentration in the cross section is deduced.

The reactive transport within simple geometries such as parallel plates, infinite tubes, or closed spheres has already been studied [24,34,35]. For example, Bekri et al. [3] provide the analytical solution between two infinite parallel plates.

As explained at the end of Sec. II B, the system is assumed to be close to chemical equilibrium and the rate of deformation of the surface is assumed to be very slow; the transitional regime is short compared to the asymptotic regime and the asymptotic regime is assumed to be reached. For a first-order reaction, this assumption is equivalent to supposing that the normalized concentration $c'$ (4) undergoes an exponential decay with time characterized by the decrease rate $\lambda$:

$$c'(\rho,t) \approx X(\rho) \exp(-\lambda t), \qquad (17)$$

where $X$ is the transverse profile of $c'$.

For an infinite cylinder, $X$ can be written as a function of the Bessel functions $J_0$ and $J_1$:

$$X(\rho') = \frac{\omega^2}{2\text{PeDa}} \frac{J_0(\omega\rho')}{J_0(\omega)} \bar{X}^t \qquad (18a)$$

$$\text{with } \bar{X}^t = \frac{1}{S} \int_S X d^2s, \qquad (18b)$$

where $\rho'$ is the dimensionless radial coordinate ($\rho' = \rho/r$). $\omega$ is the first positive solution of the equation

$$\omega J_1(\omega)/J_0(\omega) = \text{PeDa}. \qquad (19)$$

Calculation details are provided in Appendix A 1.
Using (18a), Shapiro and Brenner [35] showed that the average normalized concentration $\bar{c}'(z)$ is governed by

$$ \frac{\partial \bar{c}'}{\partial t} + \frac{\partial}{\partial z} \left( v^* \bar{c}' - D^* \frac{\partial \bar{c}'}{\partial z} \right) + \gamma^* \bar{c}' = 0 \quad (20a) $$

$$ \text{with } \bar{c}' = \frac{1}{S} \int_S c' d^2s, \quad (20b) $$

where $\gamma^*, v^*$, and $D^*$ are the apparent volume reaction rate, the apparent velocity, and the dispersion of the solute, respectively. These coefficients, which characterize the solute behavior in the pore throats, can be derived from the spatial global moments $\mathbf{m}_i$ of the concentration [36]:

$$ \gamma^* = -\frac{1}{m_0} \frac{d}{dt}(m_0), \quad (21a) $$

$$ \mathbf{v}^* = \frac{d}{dt} \left( \frac{\mathbf{m}_1}{m_0} \right), \quad (21b) $$

$$ \mathbf{D}^* = \frac{1}{2} \frac{d}{dt} \left[ \frac{\mathbf{m}_2}{m_0} - \left( \frac{\mathbf{m}_1}{m_0} \right)^2 \right], \quad (21c) $$

$$ \text{with } \mathbf{m}_i(t) = \int c(\mathbf{r},t) \mathbf{r}^i d^3\mathbf{r}, \quad (21d) $$

where $\mathbf{r}$ is the spatial position within the pore throat.

Furthermore, for capillary tubes, Algive et al. [25] provided analytical expressions of $\gamma^*, v^*$, and $D^*$ as functions of Pe and PeDa. These functions are deduced from a numerical application of the formulation of Sankarasubramanian and Gill [34] and from the propagation of a particle cloud using a random-walk technique. Thanks to these formulas, macroscopic coefficients can be assigned to each pore throat and pore body of the pore network.

Moreover, since only long-term phenomena are studied, the asymptotic regime can be generalized to the whole pore network. Thus, a unique exponential decrease rate is defined which is common to every element of the pore network. Let $\Upsilon$ be this decrease rate. Then $\bar{c}'$ in a pore throat can be written as

$$ \bar{c}'(z,t) = \bar{X}^t(z) \exp(-\Upsilon t). \quad (22) $$

Introduction of (22) into (20a) yields a second-order ordinary differential equation

$$ \frac{\partial^2 \bar{X}^t}{\partial z'^2} - Pe^t \frac{\partial \bar{X}^t}{\partial z'} - PeDa^t \bar{X}^t = 0 \quad (23) $$

$$ \text{with } Pe^t = \frac{v^*l}{D^*} \quad \text{and} \quad PeDa^t = \frac{(\gamma^* - \Upsilon)l^2}{D^*}, $$

where $l$ is the length of the pore throat and $z' = z/l$ is the dimensionless longitudinal coordinate along the tube. The distribution of the solute within pore throats can be deduced easily from this second-order ordinary differential equation as a function of $\bar{X}^t$ at the edges of the pore throat, denoted $\bar{X}^t(z' = 0)$ and $\bar{X}^t(z' = 1)$ (cf. Fig. 3).

## 2. Pore bodies

The reasoning applied to the resolution of the concentration distribution within pore throats cannot be generalized to pore bodies where the velocity is not calculated. Hence, $X(\rho')$ in pore bodies is limited to only two forms based on dominant diffusion or perfect mixing.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Schematization of a pore throat connecting the two pore bodies i and $I_{ij}$.</figcaption></figure>

Perfect mixing implies a uniform $c'$ in the pore body. For dominant diffusion, $c'$ in a sphere of radius $R$ is controlled by diffusion and reaction; therefore,

$$ X(\rho') = \frac{\omega^2}{3\text{PeDa}} \frac{\sin(\omega\rho')}{\rho \sin(\omega)} \bar{X} \quad (24a) $$

$$ \text{with } \bar{X} = \frac{1}{V} \int_V X d^3x, \quad (24b) $$

where $\rho' = \rho/R$ is the dimensionless radial coordinate. $\omega$ is the first positive solution of

$$ 1 - \frac{\omega}{\tan \omega} = \text{PeDa.} \quad (25) $$

The flux at $\rho' = 1$ is given by

$$ \Phi_w = \left(1 - \frac{\omega}{\tan \omega}\right) \frac{\omega^2}{3\text{PeDa}} \bar{X}. \quad (26) $$

Details of the calculations are provided in Appendix A 2.

## 3. Pore network

Then, the mass balance over each pore body $i$ of the network yields

$$ \text{PeDa}_i^p \bar{X}_i = \sum_{j=1}^{n_i} \Phi_{ij} \quad \text{with} \quad \text{PeDa}^p = \frac{(\gamma^* - \Upsilon)R^2}{D^*}, \quad (27) $$

where the subscripts $i$ and $j$ correspond to the $i$th and $j$th pore bodies; $\Phi_{ij}$ are the solute fluxes at the interface between the $n_j$ connected neighbor pore throats (with $n_j$ also called the coordination number) and the $i$th pore body. The left side of (27) is the sink-source term of the reaction in the pore body. $\Phi_{ij}$ is derived from the analytical solution of (23):

$$ \Phi = \left[ - \frac{\partial \bar{X}^t}{\partial z'} \Big|_{z'=0} + Pe^t \bar{X}^t(0) \right] = \varphi \bar{X}^t(0) + \psi \bar{X}^t(1), \quad (28) $$

where $\varphi$ and $\psi$ are coefficients derived from the resolution of (23).

The neighbor pore body connected through pore throat $j$ to the $i$th pore body is denoted by $I_{ij}$ (cf. Fig. 3). Using this notation, the end conditions of a pore throat become

$$ \bar{X}^t(z' = 0) = \xi_i \bar{X}_i \quad \text{and} \quad \bar{X}^t(z' = 1) = \xi_{I_{ij}} \bar{X}_{I_{ij}}, \quad (29) $$

where $\xi$ is equal to 1 when perfect mixing is assumed or equal to $\frac{\omega^2}{3\text{PeDa}}$ for dominant diffusion [$\rho' = 1$ in (24a)].

Thus, the mass balance (27) becomes

$$ \text{PeDa}_i^p \bar{X}_i = \sum_{j=1}^{n_i} (\varphi_{ij} \xi_i \bar{X}_i + \psi_{ij} \xi_{I_{ij}} \bar{X}_{I_{ij}}). \quad (30) $$
$N$ equations (where $N$ is the number of pore bodies in the network) with $N + 1$ unknowns can be written for the mass balance of each pore body. Indeed, the unknowns are $\bar{X}_i=1,...,N$ for each pore body, to which should be added the overall decrease rate of the normalized concentration $\Upsilon$ (22). Thus, a closure equation is needed and it is provided by the normalization of the dimensionless concentration (4).

Moreover, because $\varphi_{ij}$ and $\psi_{ij}$ are related to $\Upsilon$ through PeDa$^p$, (30) involves some nonlinear terms which are products of unknowns. Thus, the iterative Newton-Raphson method [33] is used to solve the resulting nonlinear system.

In the first iteration, perfect mixing is assumed in all the pore bodies. Then, during the iterative process, the solution of dominant diffusion (24a) is assigned if all the fluxes $\Phi_{ij}$ of a given pore body, computed at the previous iteration, have the same sign. At any iteration, when this minimal criterion is no longer met, the pore body switches back to perfect mixing. Of course, the closer to (26) solute fluxes $\Phi_{ij}$ are, the more accurate the dominant diffusion solution is.

### D. Algorithm description

The evolution of the geometry due to reaction is computed through an iterative process based on porosity modifications as detailed in Fig. 4. For a given initial pore network, the flow field is determined in step 1 for an arbitrary pressure difference between inlet and outlet [8]. Then, the porosity and the permeability of the pore network are calculated as well as the mean interstitial velocity. In step 2, a linear correction is applied to the pressure difference to adjust the velocity field to the imposed Pe. As in Algive et al. [25], the pore-scale transport coefficients $\gamma^*, \nu^*$, and $D^*$ defined in (20a) are determined for each pore body and each pore throat in step 3.

<figure><img src="image_5.png" /><figcaption>FIG. 4. General scheme of the reactive transport resolution by using the PNM.</figcaption></figure>

In step 4, (30) is solved over the whole network and it yields $\bar{X}$. Then, in step 5, the evolution of the geometry is taken into account with (5f). Due to the variations of $c'$, the reaction is not uniformly distributed within an element of the pore network. To be consistent with the PNM formalism, the wall evolution is averaged over each element of the pore network in order to keep the original shape of the element. The wall evolution is adjusted in order to obtain small and controlled evolution of the porosity. Since the geometry of the new pore network is compatible with the PNM, steps 1 to 5 are iterated (step 6) on the updated pore network.

At the end of the simulation, porosity-permeability curves and the macroscopic coefficients $\bar{\gamma}^*, \bar{\nu}^*$, and $\bar{D}^*$ are calculated.

## IV. MODEL COMPARISONS

In order to evaluate the accuracy of the reactive PNM, it is compared with the PSM and observations of Daccord et al. [37] at the end of the first iteration of the general schemes schematized in Fig. 1 for the PSM and Fig. 4 for the PNM. It should be noticed that the pore-throat surface displacement deduced from solute flux at the wall (3) is taken into account differently in the two models. In the PSM, the surface displacement is calculated for each voxel, while in the PNM the pore-throat diameter evolution is deduced from the average solute flux over each surface element. Thus, the evolution of the pore geometry by the PSM would not stay consistent with PNM formalism for long times; it should be noticed that this paper is not focused on this geometrical evolution and only the flux at the wall of the network is calculated.

Two samples are addressed by both numerical codes. The first sample is used to validate the analytical solutions implemented in the PNM. The second one, which contains only six pores, has been designed in order to highlight and understand the various possible reaction patterns observed in the asymptotic regime. Moreover, thanks to this reduced geometry, a new reaction regime has been found.

### A. Validation of the PNM assumptions

Consider the cubic pore network of $4 \times 4 \times 4$ spheres interconnected by capillary tubes, which is displayed in Fig. 5; it is called PN1. Each pore body is connected to six other ones. The mean flow is parallel to one of the main directions of the cubic network.

The pore-body diameters are stochastically generated with a Weibull probability density function. Pore-throat diameters are related to the smallest neighbor pore body by a ratio equal to 4 between pore-body and pore-throat diameters. The distance between pore bodies is equal to 100 µm along all directions.

The study of the $c'$ field is focused on the main flow path and on the secondary flow path which passes through the largest pore body. These two flow paths are crucial for the description and understanding of the field $c'$ (see Fig. 5). Both flow paths are in the plane $y = 100$ µm. The main flow path is the one whose axis is $z = 200$ µm and the center of the largest pore body is at $x = 300$ µm and $z = 300$ µm.
<figure><img src="image_6.png" /><figcaption>FIG. 5. Pore network PN1.</figcaption></figure>

### 1. Flow field validation

In order to minimize discrepancy between the fields $c'$, the computed flow fields by both models have to be compared.

The transverse velocity profile obtained in pore throats by the PSM can be used to validate the Poiseuille profile used in the PNM. It should be emphasized that the pore space is the one shown in Fig. 5 and discretized by the LSM. The velocity profile of the shortest pore throat obtained by using the PSM is compared in Fig. 6(a) to the analytical solution (13) used in PNM. As shown in Fig. 6(b), a good agreement is observed for the mean velocity between two pore-body centers, despite the hydraulic resistance of the pore bodies, which is not taken into account in PNM.

The reactive transport solutions in this network are calculated with local values of Pe and PeDa. For the sake of simplicity, the mean diameter of the network elements, $l_c = \langle d \rangle = 18.98$ µm, is chosen as the characteristic length. The molecular diffusion coefficient $D$ is supposed to be constant and equal to $10^{-9}$ m$^2$s$^{-1}$. When the same Pe is imposed on both models, a comparable local velocity field is obtained, as observed in Fig. 6(c).

### 2. Concentration field validation

The normalized concentration obtained by the PSM in the pore space displayed in Fig. 5 and discretized by the LSM is used to validate the hypotheses implemented in the PNM; along a pore throat, $X$ follows the analytical solution (18a) and does not depend on the longitudinal position; i.e., there is no effect of the boundaries of the pore throats. Figure 7 shows the PSM $X$ in one pore throat for PeDa = 1 and Pe = 1. Error bars in Fig. 7 represent the spreading along the pore-throat axis of the values of $X$ for the same radius. The PSM confirms that $X$ can be considered as established along the pore throats and follows the analytical solution (18a).

In Fig. 7, $X$ is normalized by $\bar{c}'$ in the pore-throat cross section. However, $\bar{c}'$ in PSM calculations depends on the longitudinal position in the pore throat. It matches the solution of the average equation (20a), which appears to be a good

<figure><img src="image_7.png" /><figcaption>FIG. 6. Velocity comparison between PSM results (□) and the PNM solution (—). (a) Comparison between the transverse velocity profiles $v'$ in a throat of PN1; (b) comparison between the longitudinal velocity profiles $\bar{v}'$ along the same pore throat. The mean velocity in each pore throat of PN1 computed by PNM is compared to PSM results in (c). The bisecting line corresponds to a perfect agreement between the two models.</figcaption></figure>

approximation for reactive transport in pore throats that can be used in PNM formalism.
<figure><img src="image_8.png" /><figcaption>FIG. 7. The transverse normalized concentration profile $X$ in the cross section of the pore throat PT1. Comparison between PSM results ($\times$) and the analytical solution provided by (18a) and implemented into the PNM (—). The small error bars relative to PSM results indicate the spreading of the transverse profile along the pore throat.</figcaption></figure>

Then, the local distributions of $c'$ obtained by using the PSM for PeDa = 1 and Pe = 0.1, 1, and 2 in the whole network are compared to PNM results. The values of the dimensionless numbers are chosen in order to induce various solute distributions and to correspond to the regimes described by Daccord et al. [37]. These reaction regimes will be discussed in the rest of this section.

For Pe equal to 0.1, diffusion is found to be dominant in the largest pore body (PB1). Therefore, the local solute distribution in this pore body is assumed to follow Eq. (24a). The PNM solution is compared to PSM calculations in Fig. 8(a). Despite some small convective effects in PB1 which cannot be taken into account by the PNM, a very good agreement between the models is observed. The pore bodies $\bar{X}$ are satisfactorily compared in Fig. 9(a). It is seen that the largest pore body contains the largest amount of solute, which is in agreement with the results of Daccord et al. [37]. Indeed, for low Pe and high PeDa, dominant diffusion implies a dominant reaction in the largest elements of the network and this induces vuggy porosity formation.

For Pe equal to 1, the condition for dominant diffusion in the largest pore body in the PNM is no longer satisfied and perfect mixing is assigned to PB1 in PNM calculations. For clarity, PNM and PSM results are compared only along the secondary flow path in Fig. 8(b). The PSM results, which are displayed in Fig. 9(b), show that convection influences the solute distribution, an effect which cannot be taken into account by the PNM. However, Fig. 9(b) shows that distributions of $\bar{X}$ in all pore bodies of the network are generally in good agreement.

For Pe equal to 2, according to the PSM, the dominant reaction in the asymptotic regime takes place along the main flow path, which should induce, with a modification of the geometry, the wormhole regime observed experimentally by Daccord et al. [37]. Indeed, for high Pe and high PeDa, the reaction occurs preferentially along the main flow path and it induces wormholing dissolution patterns. Variations of $X$

<figure><img src="image_9.png" /><figcaption>FIG. 8. Variations of the normalized concentration $X$ at PeDa = 1 along the secondary flow path for Pe = 0.1 (a) and 1 (b) and along the main flow path for Pe = 2 (c). Data are for PSM calculations (\*-) and the solution (24a) used in PNM pore bodies (—).</figcaption></figure>

along the main flow path are compared in Fig. 8(c). A ratio of roughly 3 is observed between PSM and PNM values, but the same qualitative longitudinal profiles, with little variation, are observed. This is also seen in Fig. 9(c), where the values of $\bar{X}$ in pore bodies computed by PNM and PSM are compared. Underestimation of $\bar{X}$ along the main flow path by the PNM has consequences on the whole network. Thus, by normalization, the largest pore body $\bar{X}$ is overestimated. Anywhere else, the PNM and the PSM are very close.

Therefore, when diffusion is dominant, both the local variations of $c'$ and the macroscopic behavior of the solute ($\bar{X}$ within pore bodies) obtained by the PNM and the PSM are in perfect agreement. Moreover, despite the discrepancies observed between models when convection increases, the
<figure><img src="image_10.png" /><figcaption>FIG. 9. Comparison of the average normalized concentration $\bar{X}$ between the two models (PNM and PSM) within pore bodies of the network PN1 (Fig. 5) for PeDa = 1. Data are for Pe = 0.1 (a), 1 (b), and 2 (c). The bisecting line corresponds to a perfect agreement between the two models.</figcaption></figure>

global distributions of the solute in the network are in good agreement.

### 3. Validation of core scale coefficients

The macroscopic behavior of the PSM and the PNM are compared by using the two macroscopic parameters $\bar{\gamma}^*$ and $\bar{v}'$. In the asymptotic regime, the apparent reaction rate $\bar{\gamma}^*$ is computed by using the PSM as the exponential decrease rate of $c'$ and by using the PNM as the average of $\gamma^*$ in each pore

<figure><img src="image_11.png" /><figcaption>FIG. 10. The apparent reaction rate $\bar{\gamma}^*$ (a) and the average solute velocity $\bar{v}'$ (b) in the synthetic network PN1 as functions of the Péclet number Pe. PeDa = 1. Data are for the PSM (■) and the PNM (■-●).</figcaption></figure>

body and pore throat (21a) weighted by $c'$. Therefore,

$$ \bar{\gamma}_{\text{PSM}}^* = \frac{2 \langle c'_{it} \rangle - \langle c'_{it+1} \rangle}{\delta t \langle c'_{it} \rangle + \langle c'_{it+1} \rangle}, \quad (31a) $$

$$ \bar{\gamma}_{\text{PNM}}^* = \frac{\langle \bar{c}' V \gamma^* \rangle}{\langle \bar{c}' V \rangle}, \quad (31b) $$

where $\langle \cdot \rangle$ is the volume integral over the void space of the porous medium. In the PSM, $\delta t$ is the time step and $c_{it}$ is the local field of $c'$ at iteration it. $\bar{c}'$ is determined by using the PNM [cf. (20b)]; $V$ is the volume of each element of the network. These global reaction rates are satisfactorily compared in Fig. 10(a). The observed discontinuity of $\bar{\gamma}_{\text{PNM}}^*$ at Pe ≈ 0.8 is due to the shift from dominant diffusion to perfect mixing in the pore bodies.

Let $\bar{v}'$ be the average fluid velocity within the pore space weighted by $\bar{c}'$:

$$ \bar{v}' = \frac{\langle \bar{c}' \bar{v}_x \rangle}{\langle \bar{c}' \rangle \langle \bar{v}_x \rangle}, \quad (32) $$

where $\langle \bar{c}' \rangle$ and $\langle \bar{v}_x \rangle$ are the volume averages of $c'$ and of the velocity in each element of the network determined by both the PNM and the PSM, respectively. $\langle \bar{v}_x \rangle$ is also equal to the interstitial velocity of the fluid. $\bar{v}'$ is used to evaluate the velocity of the solute on the core scale by using both the PNM and the PSM. The variations of this quantity with Pe are displayed for both models in Fig. 10(b). For the apparent reaction rate, a discontinuity of $\bar{v}_{\text{PNM}}^*$ is observed for Pe ≈ 0.8. For Pe = 2, the observed discrepancy is mostly induced by the low value of $c'$ obtained by the PNM along the main flow
<figure><img src="image_12.png" /><figcaption>FIG. 11. The normalized concentration $c'$ obtained by the PSM in PN2 for PeDa = 1. Data are for Pe = 0.1 (a), 1 (b), and 50 (c). (d) Gray scale used to display the local normalized concentration.</figcaption></figure>

path where the dominant reaction takes place. Despite the assumptions made in the PNM, a good quantitative agreement between models is observed for the macroscopic parameters.

## B. Validation of reaction regimes

The second proposed geometry, reduced to six pore bodies in a single plane and called PN2, has been devised for a better understanding of the various reaction patterns observed in the asymptotic regime. A cross section of the chosen geometry is shown in Fig. 11. Three facts are important to notice about this geometry: first, the main flow path is composed of the pore bodies PB1 and PB2; second, the largest pore body PB1 belongs to the main flow path; third, the largest pore body out of the main flow path, PB3, belongs to the secondary flow path.

Three local fields $c'$ are illustrated in Fig. 11; the cross sections of the fields using the PSM are drawn for Pe = 0.1, 1, and 50 to illustrate the three reaction regimes observed at PeDa = 1. Indeed, Fig. 11(a) illustrates the first reaction regime identified for PN1, i.e., reaction occurring within the largest pore body for low Pe. Moreover, in Fig. 11(c), a dominant reaction is observed in the main flow path for very high Pe. However, for an intermediate Pe (Pe = 1), a particular reaction regime is observed [Fig. 11(b)]. In this transitional regime, reaction occurs preferentially within the largest pore body out of the main flow path.

Of course, the governing regime has important consequences on the value of the macroscopic parameters and on the porosity-permeability evolution. The transition regime [Fig. 11(b)] would induce a lower macroscopic velocity because the fluid velocity in the main flow path is faster than in the secondary flow path where the reaction is dominant.

The variations of $\bar{X}$ in the three main pore bodies (PB1, PB2, and PB3) as a function of Pe are displayed in Fig. 12 for the PNM and the PSM.

A perfect agreement between models is observed for Pe = 0.1 in each pore body. The reaction occurs preferentially for the two models in the largest pore body, PB1.

For Pe = 1 and Pe = 10, both PNM and PSM calculations show that the reaction occurs preferentially in PB3 (Fig. 12). However, small discrepancies of concentrations appear which are due to the zero velocity assumed by the PNM in pore bodies. Nevertheless, the observed discrepancy between models remains small since the reaction mostly occurs in PB3 where diffusion is still dominant for Pe ≪ 10.

For Pe = 50, PSM calculations show that reaction occurs preferentially along the main flow path and that it induces wormholing dissolution patterns [Fig. 11(c)]; indeed, $\bar{X}_{PB1}$ and $\bar{X}_{PB2}$ are larger than $\bar{X}_{PB3}$, as observed in Fig. 12. Moreover, according to the PNM, the wormholing reaction regime is not fully established for Pe = 50 since the network still undergoes a dominant reaction in PB3. However, the larger Pe becomes, the more dominant the main flow path becomes until $\bar{X}_{PB1}$ and $\bar{X}_{PB2}$ overtake $\bar{X}_{PB3}$ for Pe ≈ 200. The observed discrepancy is due to the PNM, where the pore-body hydraulic resistance is neglected. As a consequence, higher Pe is needed to observe a dominant reaction along the main flow path in the PNM.

Therefore, it is important to notice that the same three reaction regimes are obtained by both models.

## C. Discussion

### 1. Model complementarity

The proposed study shows a good agreement between the two models despite their different application scales. The same qualitative behavior and the same regime changes are observed by using both models for the studied networks. Local resolution of flow and concentration allows a better understanding of the limitations of PNM assumptions. For low Pe, a very good quantitative agreement between models is observed. Thus, in diffusive conditions, the simplifying assumptions of the PNM are justified and compatible with the exact PSM calculations.

For high Pe, the strong assumptions of the PNM relative to pore bodies can lead to significant discrepancies in $c'$ when compared to the PSM. One should notice that, generally in the
<figure><img src="image_13.png" /></figure>

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

<figure><img src="image_15.png" /><figcaption>FIG. 12. The average normalized concentration $\bar{X}$ as a function of Pe in (a) the largest pore body (PB1), (b) the second pore body along the main flow path (PB2), and (c) the largest pore body out of the main flow path (PB3). Data are for the PSM (■) and the PNM (—).</figcaption></figure>

PNM, the porous medium is not modeled by a regular lattice aligned with the pressure gradient ($\theta = 0^\circ$). This kind of pore-network construction is not able to take into account the totuosity of a real porous medium. A classical way to avoid this problem is to consider a regular lattice rotated by $\theta = 45^\circ$ with respect to the applied pressure gradient, or to consider an irregular lattice. For these porous media, the comparison between PSM and PNM results should be even better since the perfect mixing assumption in pore bodies would be more justified.

## 2. Reaction regimes

According to Daccord et al. [37], three main regimes occur during a surface reaction depending on Pe and PeDa numbers: a uniform reaction for low PeDa, a flow path reaction for high PeDa and high Pe (inducing wormholing), and a compact

reaction for high PeDa and low Pe (inducing vuggy porosity formation). Thus, for high PeDa, this classification predicts a change of reaction regime with Pe, from compact reaction to wormholing. This change is confirmed by both PNM and PSM wall flux calculations performed on PN1 (Sec. IV A). The evolution of the reaction regime with Pe results from the competition, in terms of $c'$, between the largest pore body and the main flow path. However, in PN2, where the largest pore body is on the main flow path, an intermediary reaction regime appears between compact and wormholing reaction regimes. To the best of our knowledge, this specific regime was not observed in the literature.

For a better understanding of conditions for this new intermediary regime to appear, the main and secondary flow paths of PN2, denoted $FP_M$ and $FP_S$, respectively, in Fig. 13, are independently studied by using the PNM. The flow path with the slowest $\bar{\gamma}^*$ is likely to dominate the solute transport. In Fig. 13(c), the reaction rates $\bar{\gamma}^*$ of $FP_M$ and $FP_S$ are plotted as functions of Pe (cf. Fig. 12); two crossing points are noticed ($Pe \approx 0.25$ and $Pe \approx 300$). Thus, for a given Pe, the local velocities in $FP_M$ and $FP_S$ are the same as in the corresponding flow paths in PN2. For low Pe ($Pe < 0.25$), $c'$ in PN2 is localized essentially in the main flow path since the reaction rate of $FP_M$ is slower than that of $FP_S$ [Fig. 13(c)]; the same behavior is observed in Fig. 12. For moderate $Pe$ ($0.25 < Pe < 300$), $c'$ in PN2 is localized essentially in $FP_S$ because its reaction rate is slower than the $FP_M$ one [Fig. 13(c)]; the same behavior is also observed in Fig. 12. For high Pe ($Pe > 300$), $c'$ in PN2 is localized along $FP_M$ because its rate is again slower than the $FP_S$ one [Fig. 13(c)], as observed in Fig. 12.

By considering only the calculations performed either on $FP_M$ or on $FP_S$, one can notice the increase of the apparent reaction rate $\bar{\gamma}^*$ with Pe, which corresponds to an acceleration of the chemical reaction. For a better understanding of this phenomenon, one should notice that, for high PeDa, the mean reaction rate of each element of the network (pore body or pore throat) is inversely related to the square of its radius. Indeed, for high PeDa, $\omega$, which is the solution of (A7) [or (A12) depending on the shape of the element], tends to a constant. Thus, $\bar{\gamma}^*$, which is equal to $\lambda$ in (A3), is related to the radius of an element through (A6) (see Appendix A for more details). Therefore, the larger an element, the slower the reaction rate (because the reaction is limited by diffusion). Nevertheless, for low Pe, solute exchanges between elements are limited; each element reacts according to its own rate $\bar{\gamma}^*$. Since the element (or group of elements) with the slowest reaction rate dominates the solute transport in the asymptotic regime, the reaction occurs in the larger pore body for low Pe. Thus, the calculated $\bar{\gamma}^*$ is slower than for any higher value of Pe.

Moreover, the exchanges of solute along $FP_M$ and $FP_S$ increase with Pe. Thus, $c'$ in the largest pore body of each flow path is convected out and consumed by the nearby pore throats. Since these elements are necessarily smaller and their reaction rate is faster, a speed up of the apparent reaction rate of the flow path is induced; for instance, this acceleration is observed in Fig. 13(c) around $Pe = 1$ and $Pe = 100$ for $FP_M$ and $FP_S$, respectively. Therefore, $\bar{\gamma}^*$ is an increasing function of Pe.

Finally, when convection is dominant over diffusion within $FP_M$ or $FP_S$ (for very high Pe), the distribution of $c'$ is uniform
<figure><img src="image_16.png" /><figcaption>FIG. 13. The main flow path (a) and the secondary flow path (b) of PN2. (c) The apparent reaction rate $\bar{\gamma}^*$ as a function of Pe for the main (■) and secondary (—) flow paths.</figcaption></figure>

along the flow path. The obtained reaction rate, denoted $\bar{\gamma}_c^*$, is the quickest reaction rate that can be reached by the considered geometry. Since $c'$, in this case, is uniform, $\bar{\gamma}_c^*$ can be evaluated by using (31b) as the volume average of $\gamma^*$ along the considered flow path (FP):

$$ \bar{\gamma}_c^* = \sum_{FP} V\gamma^* / \sum_{FP} V. \qquad (33) $$

Thus, the intermediary reaction regime, highlighted by the study of PN2, appears when two conditions are satisfied. First, the largest pore body of the network must belong to the main flow path. Second, at least one element of the pore network out of the main flow path must have a reaction rate slower than the critical reaction rate $\bar{\gamma}_c^*$ of the main flow path.

Since $\gamma^*$ is inversely related to the square of the diameter D of the element for high PeDa, the above conditions can be formulated by using element diameters. A critical diameter $D_c$ is defined as the diameter of a pore body where the solute reacts at the same rate as on the main flow path of the network. The second condition necessitates that at least one pore body diameter is greater than $D_c$:

$$ \forall m \notin FP_M, \quad D_m < \max_{FP_M} D, \qquad (34a) $$

$$ \exists m \notin FP_M, \quad D_m > D_c = \sqrt{\sum_{FP_M} V / \sum_{FP_M} \frac{V}{D^2}}. \qquad (34b) $$

Of course, these two conditions are verified in PN2. The largest pore body is on the main flow path with a diameter of $D_{max} = 65$ µm; the critical diameter is $D_c \approx 55$ µm and $D_{PB3} = 60$ µm is the diameter of the largest pore body out of the main flow path. Thus, this intermediate regime disappears when the diameter of PB3 is smaller than 55 µm. Figure 14 shows the mean $c'$ calculated by using the PNM in the modified PN2 (PB3 = 55 µm instead of 60 µm). Whatever the value of Pe, $c'$ in PB1 is always higher than in the new PB3 and there is no longer an intermediate reaction regime. On the one hand, a compact reaction is observed at low Pe; diffusion is dominant in PB1 and its mean $c'$ is high. On the other hand, for high Pe,

the network undergoes a wormholing regime; perfect mixing is assumed in PB1 and PB3 and their $c'$ values are around 1.

## V. APPLICATION OF THE PNM TO A PORE NETWORK EXTRACTED FROM A REAL ROCK SAMPLE

In this section, the reactive flow modeling using the PNM is applied to an irregular 3D pore network that is extracted from microtomography images (Fig. 15). Due to the complexity of real porous media, the mixing in pore bodies is improved and the local tortuosity of pore throats reduces the relative contribution of pore bodies to permeability. Therefore, in such a geometry, the PNM assumptions are more realistic. The imaged sample is used to calculate the macroscopic coefficients of the reactive transport.

The studied network is extracted from a sample of Fontainebleau sandstone of (3000 µm)$^3$ with a porosity of 14.1%, an absolute permeability of 765 mD, and a formation factor equal to 32. The selected sample is homogeneous with an average value of pore-body and of pore-throat diameters equal to 90 and 25 µm, respectively (Fig. 16). Only the pore structure of the Fontainebleau sandstone is taken into account. For simplicity, the chemical interaction at the fluid-rock interface is modeled by a constant and homogeneous $\kappa$.

<figure><img src="image_17.png" /><figcaption>FIG. 14. The average normalized concentration $\bar{X}$ calculated by using the PNM in the modified PN2 as a function of Pe. In contrast to the original PN2, the PB3 diameter is reduced in order not to satisfy (34b). Data are for PB1 (■) and PB3 (—).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIG. 15. Geometry of a pore network extracted from a real rock sample.</figcaption></figure>

By using the PNM, the reactive transport is solved and the apparent reaction rate and the effective solute velocity are computed from (31b) and (B1). These macroscopic coefficients are presented in Figs. 17(a) and 17(b). $\tilde{\nu}'$ is the apparent reaction rate on the core scale normalized by the reaction rate in perfect mixing conditions, whereas $\bar{v}^*/\langle v \rangle$ is the normalization of the effective solute velocity by the mean interstitial velocity of the fluid (see Appendix B).

At low PeDa, the solute is uniformly distributed; i.e., the solute is not influenced by the structure of the rock ($\tilde{\nu}'$ and $\bar{v}'$ equal one). Under these conditions, the reaction is so slow that the solute has time to spread by diffusion over the whole porous medium and a uniform reaction is induced.

The influence of the structure of the rock on the macroscopic coefficients is an increasing function of PeDa. It induces the decrease of the apparent reaction rate and velocity. A discontinuity of $\bar{v}^*$ is observed in Fig. 17(b) for Pe around 2. This collapse of $\bar{v}^*$ is due to the appearance of the intermediate reaction regime highlighted in Sec. IV B.

A remark can be made on this case about the used computation resources. Because of the distribution of diameters in this network (Fig. 16) and of the large volume, a cubic discretization for the PSM would imply a lattice of roughly $10^9$ cells. Such a volume cannot be addressed yet by the PSM.

Using the PNM, roughly 30 000 elements are needed to define the network geometry. Nevertheless, around 3 computation days and over 2.4 gigabytes of RAM are necessary to obtain the results displayed in Figs. 17(a) and 17(b).

<figure><img src="image_19.png" /><figcaption>FIG. 16. The probability density functions of pore-body and pore-throat diameters in the network displayed in Fig. 15.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>FIG. 17. The dimensionless effective reaction rate $\tilde{\nu}'$ (a) and the effective solute velocity $\bar{v}^*/\langle v \rangle$ (b) of an irregular network as functions of the dimensionless numbers Pe and PeDa.</figcaption></figure>

## VI. CONCLUDING REMARKS

Two methods are used to solve reactive transport in a porous medium, the PSM, which solves on the microscale, and the PNM, which reformulates this problem on a larger scale. When the simulations for both methods are made in identical media under the same conditions, the two models provide analogous solutions. The present investigation validates the hypotheses and the formulations used in the PNM in order to describe flow, surface reaction, and transport in each element of the network. PSM calculations show that $v$ and $X$ can be considered as established along the pore throats and follow the analytical solutions for an infinite capillary. Despite the PNM assumptions in pore bodies—zero velocity and perfect mixing or dominant diffusion for $c'$—a very good agreement between the models is observed.

From a computing point of view, for PN1, the PSM requires 3 gigabytes during 2 days while the PNM only needs 100 megabytes and a few seconds on the same computer. Because of the use of the analytical solutions in the network elements,
the PNM can be applied to large velocities and reaction rates. However, this approach cannot represent exactly the geometry and the local evolutions of the solute. The PNM only provides an average description of solute distribution in pore bodies and pore throats. Because of the simplifications used in the PNM, a large volume can be addressed with affordable computation resources.

A new reaction regime, not mentioned in the literature, is brought into light which depends on the localization and the diameter of the largest elements of the porous medium with respect to the main flow path. This new reaction regime can have consequences on the macroscopic behavior of the solute which can change the distribution pattern of the concentration in reservoir simulators [38].

## APPENDIX A: ANALYTICAL RESOLUTION OF THE TRANSVERSE PROFILE

This Appendix is devoted to the resolution of reactive transport in infinite capillary tubes and in closed spheres. As a reminder, the reactive transport of a solute is governed by

$$ \mathrm{Pe} \nabla' c' \cdot \mathbf{v}' - \nabla'^2 c' = - \frac{L^2}{DT} \frac{\partial c'}{\partial t} \quad \text{in } \Omega_F, \quad (A1) $$

$$ \mathbf{n} \cdot \nabla' c' = -\mathrm{PeDa} c' \quad \text{on } \Gamma. \quad (A2) $$

### 1. Infinite capillary tube

In the asymptotic regime, (A1) becomes a well-known second-order ordinary differential equation because $c'(\rho,t)$ is decomposed as

$$ c'(\rho,t) = X(\rho) \exp(-\lambda t). \quad (A3) $$

For an infinite cylinder, the local problem is governed by

$$ \frac{1}{\rho'} \frac{d}{d\rho'} \left( \rho' \frac{dX}{d\rho'} \right) + \frac{\lambda r^2}{D} X(\rho') = 0, \quad (A4) $$

whose solutions are linear combinations of Bessel functions of first and second kinds (commonly denoted $J_0$ and $Y_0$, respectively), where $r$ is the radius of the cylinder and $\rho'$ is the dimensionless radial coordinate ($\rho' = \rho/r$). However, $Y_0$ is discontinuous at zero, which is not compatible with the continuity of the concentration profile. Thus, the normalized concentration profile is

$$ X(\rho') = a J_0(\omega \rho'), \quad (A5) $$

where $\omega$ is a dimensionless factor related to the exponential decrease rate $\lambda$; $\omega$ depends on the surface reaction rate. Indeed, by substitution of (A5) into (A4), the following relation between $\omega$ and $\lambda$ is obtained:

$$ \lambda = \frac{D\omega^2}{r^2}. \quad (A6) $$

Moreover, by substituting the local solution (A5) into the boundary condition equation (A2), a relation can be deduced between $\omega$ and the local PeDa of the cylinder:

$$ \frac{\omega J_1(\omega)}{J_0(\omega)} = \frac{\kappa r}{D} = \mathrm{PeDa}. \quad (A7) $$

Since a general analytical solution of this equation does not exist, $\omega$ must be numerically determined for a given PeDa. In order to deal with a dimensionless profile $X(\rho')$, its mean value over a section of the cylinder is fixed to 1. This leads to the final expression of $X$ for a cylinder:

$$ X(\rho') = \frac{\omega^2}{2\mathrm{PeDa}} \frac{J_0(\omega\rho')}{J_0(\omega)}. \quad (A8) $$

### 2. Closed sphere

Within a closed sphere, the local problem can be written as

$$ \frac{1}{\rho'^2} \frac{d}{d\rho'} \left( \rho'^2 \frac{dX}{d\rho'} \right) + \frac{\lambda R^2}{D} X(\rho') = 0, \quad (A9) $$

where $R$ is the radius of the sphere and $\rho'$ is the dimensionless radial coordinate ($\rho' = \rho/R$). The physical solution of this equation is

$$ X(\rho') = a \frac{\sin(\omega\rho')}{\omega\rho'}. \quad (A10) $$

By substituting Eq. (A10) into (A9), an analogous relation to (A6) is obtained:

$$ \lambda = \frac{D\omega^2}{R^2}. \quad (A11) $$

Since the same reasoning can be used for the boundary condition, the following relation is deduced:

$$ 1 - \frac{\omega}{\tan \omega} = \frac{\kappa R}{D} = \mathrm{PeDa}. \quad (A12) $$

Similarly to the cylinder geometry, the determination of $\omega$ requires a numerical evaluation for a given PeDa. The final expression of $X$ within a closed sphere is obtained by normalizing the mean concentration over the whole sphere:

$$ X(\rho') = \frac{\omega^2 \sin(\omega\rho')}{3\mathrm{PeDa} \rho' \sin(\omega)}. \quad (A13) $$

## APPENDIX B: MACROSCOPIC SOLUTE VELOCITY

According to Shapiro and Brenner [39], the macroscopic velocity of the solute $\bar{v}^*$ can be expressed as

$$ \bar{v}^* = \frac{\langle X_0 Y_0 \mathbf{v} + D(Y_0 \nabla X_0 - X_0 \nabla Y_0) \rangle}{\langle X_0 Y_0 \rangle}, \quad (B1) $$

where $X_0$ and $Y_0$ are the solutions of the local reactive transport problem and its adjoint problem, respectively. This expression is adapted to the PNM as proposed by Varloteaux et al. [38]. The macroscopic velocity of the solute is used to represent the advective contribution of the reactive transport at a macroscopic scale. Of course, this velocity is not equal to $\bar{v}'$ defined in Sec. IV A.
[1] O. Izgec, B. Demiral, H. Bertin, and S. Akin, Transp. Porous Media 72, 1 (2008).

[2] O. Izgec, B. Demiral, H. Bertin, and S. Akin, Transp. Porous Media 73, 57 (2007).

[3] S. Bekri, J.-F. Thovert, and P. Adler, Chem. Eng. Sci. 50, 2765 (1995).

[4] C. L. Phillips, MIT Undergrad. J. Math. 1, 155 (1999).

[5] J. A. Sethian and P. Smereka, Annu. Rev. Fluid Mech. 35, 341 (2003).

[6] A. C. Payatakes and M. M. Dias, Rev. Chem. Eng. 2, 85 (1984).

[7] B. W. Lindquist, A. Venkatarangan, J. Dunsmuir, and T.-F. Wong, J. Geophys. Res. 105, 21509 (2000).

[8] C. Laroche and O. Vizika, Transp. Porous Media 61, 77 (2005).

[9] A. Al-Kharusi and M. Blunt, J. Petrol. Sci. Eng. 56, 219 (2007).

[10] A. Raoof and S. Hassanizadeh, Transp. Porous Media 81, 391 (2009).

[11] Q. Kang, P. C. Lichtner, H. S. Viswanathan, and A. I. Abdel-Fattah, Transp. Porous Media 82, 197 (2009).

[12] I. Fatt, Trans. AIME 207, 144 (1956).

[13] M. Blunt, D. Zhou, and D. Fenwick, Transp. Porous Media 20, 77 (1995).

[14] M. I. J. van Dijke and K. S. Sorbie, Phys. Rev. E 66, 046302 (2002).

[15] P. Øren, J. Petrol. Sci. Eng. 39, 177 (2003).

[16] C. Lu and Y. C. Yortsos, Ind. Eng. Chem. Res. 43, 3008 (2004).

[17] J. Ochi, Transp. Porous Media 37, 303 (1999).

[18] M. Sahimi, A. R. Mehrabi, N. Mirzaee, and H. Rassamdana, Transp. Porous Media 41, 325 (2000).

[19] S. Bhat, in SPE Western Regional Meeting (Bakersfield, USA, 1998), p. 46210.

[20] C. Lu and Y. C. Yortsos, Phys. Rev. E 72, 036201 (2005).

[21] A. Raoof, S. M. Hassanizadeh, and A. Leijnse, Vadose Zone J. 9, 624 (2010).

[22] L. Li, C. A. Peters, and M. Celia, Adv. Water Resour. 29, 1351 (2006).

[23] D. Kim, Scale-up of reactive flow through network flow modelling (State University of New York at Stony Brook, New York, USA, 2009).

[24] L. Algive, S. Bekri, M. Robin, and O. Vizika, in *Proceeding of the International Symposium of the Society of Core Analysts*, 10 (SCA, Calgary, 2007).

[25] L. Algive, S. Bekri, and O. Vizika, SPE J. 15, 618 (2010).

[26] J. A. Sethian, Proc. Nat. Acad. Sci. USA 93, 1591 (1996).

[27] C. Lu and Y. C. Yortsos, AIChE J. 51, 1279 (2005).

[28] S. Bekri and O. Vizika, in *Proceeding of the International Symposium of the Society of Core Analysts*, 22 (SCA, Trondheim, 2006).

[29] S. Youssef, E. Rosenberg, N. Gland, S. Bekri, and O. Vizika, in *Proceeding of the International Symposium of the Society of Core Analysts*, 17 (SCA, Calgary, 2007).

[30] P. E. Øren and S. Bakke, Transp. Porous Media 46, 311 (2002).

[31] R. Sok, M. Knackstedt, A. Sheppard, W. Pinczewski, W. Lindquist, A. Venkatarangan, and L. Paterson, Transp. Porous Media 46, 345 (2002).

[32] J.-F. Thovert, J. Salles, and P. M. Adler, J. Microsc. 170, 65 (1993).

[33] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, *Numerical Recipes in C: The Art of Scientific Computing*, 2nd ed. (Cambridge University Press, Cambridge, England, 1992).

[34] R. Sankarasubramanian and W. N. Gill, Proc. R. Soc. London, Ser. A 333, 115 (1973).

[35] M. Shapiro and H. Brenner, Chem. Eng. Sci. 41, 1417 (1986).

[36] H. Brenner, Philos. Trans. R. Soc. London A 297, 81 (1980).

[37] G. Daccord, O. Lietard, and R. Lenormand, Chem. Eng. Sci. 48, 179 (1993).

[38] C. Varloteaux, S. Békri, and P. M. Adler, Adv. Water Resour. 53, 87 (2013).

[39] M. Shapiro and H. Brenner, Chem. Eng. Sci. 43, 551 (1988).