# Towards a better comprehension of reactive transport coupling experimental and numerical approaches

Daniela Bauer, Théo Briolet, Mathilde Adelinet, Maxime Moreaud, Olivier

Sissmann, Maxime Pelerin, Jérôme Fortin, Elisabeth Bemer

► To cite this version:

Daniela Bauer, Théo Briolet, Mathilde Adelinet, Maxime Moreaud, Olivier Sissmann, et al.. Towards a better comprehension of reactive transport coupling experimental and numerical approaches. Science and Technology for Energy Transition, 2024, 79, pp.22. 10.2516/stet/2024010. hal-04532718

HAL Id: hal-04532718

https://ifp.hal.science/hal-04532718

Submitted on 4 Apr 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
# Towards a better comprehension of reactive transport coupling experimental and numerical approaches

Daniela Bauer, Théo Briolet, Mathilde Adelinet, Maxime Moreaud, Olivier

Sissmann, Maxime Pelerin, Jerome Fortin, Elisabeth Bemer

► To cite this version:

Daniela Bauer, Théo Briolet, Mathilde Adelinet, Maxime Moreaud, Olivier Sissmann, et al. Towards a better comprehension of reactive transport coupling experimental and numerical approaches. Science and Technology for Energy Transition, 2024, 79, pp.22. 10.2516/stet/2024010. hal-04532718

HAL Id: hal-04532718

https://ifp.hal.science/hal-04532718v1

Submitted on 4 Apr 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
REGULAR ARTICLE

OPEN ACCESS

Towards a better comprehension of reactive transport coupling experimental and numerical approaches

Daniela Bauer<fnref n="1" />,*, Théo Briolet<fnref n="1" />, Mathilde Adelinet<fnref n="1" />, Maxime Moreaud<fnref n="2" /> ID, Olivier Sissmann<fnref n="1" />, Maxime Pelerin<fnref n="1" />, Jérôme Fortin<fnref n="3" />, and Elisabeth Bemer<fnref n="1" /> ID

Received: 7 September 2023 / Accepted: 31 January 2024

**Abstract.** In this work we focus on further understanding reactive transport in carbonate rocks, in particular limestones characterized by a bimodal pore size distribution. To this end, we performed injection experiments with CO₂-saturated water on a sample of Euville limestone and monitored the experiments with a medical CT scanner. Microscanner imaging was performed before and after alteration. Experiments showed that permeability increased by nearly two decades due to the alteration process. This increase could be attributed to the formation of a preferential dissolution path visualized on the CT images. Microscanner images show that preferential dissolution areas are characterized by the presence of numerous enlarged macropores. The preferential dissolution path created therefore retains a porous structure and does not correspond to a wormhole-type channel. To provide further knowledge of the small-scale physics of reactive transport, we performed Lattice-Boltzmann simulations of flow in a numerically generated model 2D porous medium having geometrical and topological features designed to approach Euville limestone. We showed that the fluid velocity increased in nearly percolating paths of macropores. Considering the experiments, this means that the CO₂-saturated water starts to enter high-velocity zones earlier than low-velocity zones, inducing an earlier onset of the alteration process and a more pronounced local dissolution. However, numerical results showed that the alteration of non-connected macropores leads to an increase of permeability much smaller than the experimentally observed one. To explain this fact we used effective medium modelling that permits predicting the variation in permeability as a function of the fraction of macropores and consequently as a function of alteration. It proved that as long as there is no alteration-induced percolating path consisting of macropores, the increase in permeability is relatively low as shown by the Lattice-Boltzmann simulations. An increase in permeability of several orders of magnitude is only observed when the macroporosity is close to the percolation threshold. This fact is in accordance with the experimentally observed results.

**Keywords:** Reactive transport, Lattice-Boltzmann simulation, Permeability, Homogenization.

# 1 Introduction

As reaffirmed in the IPCC 2022 report [1], Carbon Capture and Storage (CCS) will play a key role in reducing total net carbon emissions and thus mitigating global climate change. Deep saline aquifers, including carbonate aquifers, offer the greatest potential storage capacity. However, the injection of large amounts of CO₂ into carbonate formations can result in various chemical and physical interactions inducing changes in the host rock properties that need to be assessed to ensure the storage integrity [2, 3].

The injected CO₂ will notably dissolve into the *in situ* brine to form a weak acid likely to react with carbonate minerals through dissolution and/or precipitation mechanisms modifying the porous structure of the host formation [4]. These processes involve strong couplings between geochemistry and fluid transport.

Carbonate rocks are characterized by complex microstructure resulting from the initial sedimentary processes and further diagenetic modifications [5]. They often present a bimodal porous structure including macropores and micropores. Micropores are often too small to be imaged by microscanner analysis [6], carbonate samples then appear as consisting of large macropores, a microporous matrix, and a dense solid phase. The topology of

<footnote n="1">¹IFP Energies Nouvelles, 1 et 4 Av. Bois Préau, 92852 Rueil Malmaison, France</footnote>

<footnote n="2">²IFP Energies Nouvelles, Rond-point de l'échangeur de Solaize, 69360 Solaize, France</footnote>

<footnote n="3">³Laboratoire de Géologie, Ecole Normale Supérieure - PSL Research University, CNRS, UMR8538, Paris, France</footnote>

<footnote>* Corresponding author: daniela.bauer@ifpen.fr</footnote>
REGULAR ARTICLE

OPEN ACCESS

# Towards a better comprehension of reactive transport coupling experimental and numerical approaches

Daniela Bauer<fnref n="1" />,*, Théo Briolet<fnref n="1" />, Mathilde Adelinet<fnref n="1" />, Maxime Moreaud<fnref n="2" /> ID, Olivier Sissmann<fnref n="1" />, Maxime Pelerin<fnref n="1" />, Jérôme Fortin<fnref n="3" />, and Elisabeth Bemer<fnref n="1" /> ID

Received: 7 September 2023 / Accepted: 31 January 2024

**Abstract.** In this work we focus on further understanding reactive transport in carbonate rocks, in particular limestones characterized by a bimodal pore size distribution. To this end, we performed injection experiments with CO₂-saturated water on a sample of Euville limestone and monitored the experiments with a medical CT scanner. Microscanner imaging was performed before and after alteration. Experiments showed that permeability increased by nearly two decades due to the alteration process. This increase could be attributed to the formation of a preferential dissolution path visualized on the CT images. Microscanner images show that preferential dissolution areas are characterized by the presence of numerous enlarged macropores. The preferential dissolution path created therefore retains a porous structure and does not correspond to a wormhole-type channel. To provide further knowledge of the small-scale physics of reactive transport, we performed Lattice-Boltzmann simulations of flow in a numerically generated model 2D porous medium having geometrical and topological features designed to approach Euville limestone. We showed that the fluid velocity increased in nearly percolating paths of macropores. Considering the experiments, this means that the CO₂-saturated water starts to enter high-velocity zones earlier than low-velocity zones, inducing an earlier onset of the alteration process and a more pronounced local dissolution. However, numerical results showed that the alteration of non-connected macropores leads to an increase of permeability much smaller than the experimentally observed one. To explain this fact we used effective medium modelling that permits predicting the variation in permeability as a function of the fraction of macropores and consequently as a function of alteration. It proved that as long as there is no alteration-induced percolating path consisting of macropores, the increase in permeability is relatively low as shown by the Lattice-Boltzmann simulations. An increase in permeability of several orders of magnitude is only observed when the macroporosity is close to the percolation threshold. This fact is in accordance with the experimentally observed results.

**Keywords:** Reactive transport, Lattice-Boltzmann simulation, Permeability, Homogenization.

# 1 Introduction

As reaffirmed in the IPCC 2022 report [1], Carbon Capture and Storage (CCS) will play a key role in reducing total net carbon emissions and thus mitigating global climate change. Deep saline aquifers, including carbonate aquifers, offer the greatest potential storage capacity. However, the injection of large amounts of CO₂ into carbonate formations can result in various chemical and physical interactions inducing changes in the host rock properties that need to be assessed to ensure the storage integrity [2, 3].

The injected CO₂ will notably dissolve into the *in situ* brine to form a weak acid likely to react with carbonate minerals through dissolution and/or precipitation mechanisms modifying the porous structure of the host formation [4]. These processes involve strong couplings between geochemistry and fluid transport.

Carbonate rocks are characterized by complex microstructure resulting from the initial sedimentary processes and further diagenetic modifications [5]. They often present a bimodal porous structure including macropores and micropores. Micropores are often too small to be imaged by microscanner analysis [6], carbonate samples then appear as consisting of large macropores, a microporous matrix, and a dense solid phase. The topology of

<footnote n="1">¹IFP Energies Nouvelles, 1 et 4 Av. Bois Préau, 92852 Rueil Malmaison, France</footnote>

<footnote n="2">²IFP Energies Nouvelles, Rond-point de l'échangeur de Solaize, 69360 Solaize, France</footnote>

<footnote n="3">³Laboratoire de Géologie, Ecole Normale Supérieure - PSL Research University, CNRS, UMR8538, Paris, France</footnote>

<footnote>* Corresponding author: daniela.bauer@ifpen.fr</footnote>
the macropores can vary from one sample to another. Indeed, on one hand, macropores might be linked or partially connected forming a more or less percolating network. On the other hand, macropores can be isolated and embedded within the matrix, which then forms the percolating network including the macropores. The specific structure of each sample strongly influences the flow and transport patterns within the porous medium and hence the alteration processes.

In the laboratory, the impact of CO₂ injection in carbonate formations is generally studied through coreflight experiments that can be coupled with scanner imaging to derive 3D dissolution patterns [7–10]. Coreflight tests performed under various experimental conditions show the development of different dissolution patterns that can be classified into three main categories: 1) compact dissolution characterized by extensive dissolution close to the injection face, 2) wormhole development characterized by the development of localized dissolution channels, and 3) uniform dissolution characterized by a roughly homogeneous dissolution within the sample [8]. These distinct dissolution profiles are naturally associated with different evolutions of porosity and permeability. However, the impact of the rock microstructure on the development of alteration is not yet fully understood [11–13].

In order to provide a further understanding of the reactive transport processes and the resulting alteration of the rock, we combine experimental techniques with direct numerical simulations and a homogenization approach. Indeed, numerical simulation allows the determination of the velocity field in the porous structure, providing valuable information on the physics at the pore scale. Strictly speaking, in order to obtain the complete velocity field of a bimodal porous medium, it is necessary to solve the Stokes equation in the entire domain. Nevertheless, beyond the associated computational cost, this approach can be defeated if the micropores are too small to be described through imaging techniques. We overcome this deficiency by solving a specific equation, the Darcy-Brinkman equation [14, 15], allowing the description of flow in a bimodal porous medium. This equation can be seen as an interpolation between the Stokes equation, valid in the free fluid, and the Darcy equation, applied in the microporous phase. The major advantage of the Darcy-Brinkman equation lies in its simplicity of implementation as there is only one domain, meaning that the same equation is solved at any point of the domain. Simulations are performed in 2D images allowing the visualization of the flow field in the porous structure. The 2D simulations allow a refined representation of the flow, but their computational cost is too high to consider a systematic calculation of permeability, and the 3D nature of this property should also be taken into account. To this end, we propose to use a three-dimensional homogenization approach that will allow us to derive the initial permeability and its evolution under the effect of weathering.

The objective of this article is threefold: In the first step, we perform reactive transport experiments and analyze the resulting dissolution patterns through the scanner and

microscanner imaging. Then, in order to gain a further understanding of these experiments we conduct small-scale simulations to provide information on the local velocity field. Finally, we use a homogenization method to predict the variation in effective permeability due to a modification of the porous structure induced by alteration. The article is structured as follows: In the first step, we present the experimental techniques, the equation to be solved the numerical method to perform small-scale simulations, and the homogenization method to compute the permeability. Then, we discuss the corresponding results. Further details on the numerical scheme can be found in the Appendix A.

# 2 Methods

## 2.1 Alteration experiments

### 2.1.1 Sample description

In order to focus on the impact of the microstructure, we choose to study a limestone composed almost entirely of calcite: Euville limestone [16]. This limestone is Oxfordian in age and comes from the quarry of Euville located in the eastern part of the Paris Basin (Lorraine region). It consists of a crinoidal grainstone (Fig. 1) mostly made of crinoid ossicles (more than 95%), a few debris of echinoids, and rare lamellibranch bioclasts. The main cement phase corresponds to an important development of syntaxial overgrowths around crinoid debris. This cement partially occluded the primary interparticle macroporosity. A secondary intra-particle patchy microporosity is also clearly observed in the crinoid ossicles. Finally, the micritic sediment may present a microporosity (possibly intercrystalline). Measurements carried out on samples from different blocks show that Euville limestone has a porosity between 12% and 20% and a permeability between 10 and 150 mD. Its bimodal pore size distribution is confirmed by both mercury porosimetry and nuclear magnetic resonance, which lead to similar volume fractions for macroporosity and microporosity of about 65% and 35% (±5%) respectively [16].

### 2.1.2 Experimental workflow

The studied sample has a diameter of 40 mm and a length of 100 mm. The same sequence of analyses is performed before and after alteration to precisely characterize the impact of the injection of CO₂-saturated water: scanner and microscanner imaging and measurement of porosity and permeability. Medical CT scanner imaging provides a 3D image of the whole sample with a resolution of 0.625 mm. A porosity profile along the sample length can moreover be inferred from a correlation between the radiological densities measured on successive 0.625-mm-thick slices and data corresponding to pure calcite samples. Microscanner imaging allows a significantly higher resolution (20.28 µm for a 40-mm diameter sample) but requires four series of acquisition with an overlapping to image the
the macropores can vary from one sample to another. Indeed, on one hand, macropores might be linked or partially connected forming a more or less percolating network. On the other hand, macropores can be isolated and embedded within the matrix, which then forms the percolating network including the macropores. The specific structure of each sample strongly influences the flow and transport patterns within the porous medium and hence the alteration processes.

In the laboratory, the impact of CO₂ injection in carbonate formations is generally studied through coreflight experiments that can be coupled with scanner imaging to derive 3D dissolution patterns [7–10]. Coreflight tests performed under various experimental conditions show the development of different dissolution patterns that can be classified into three main categories: 1) compact dissolution characterized by extensive dissolution close to the injection face, 2) wormhole development characterized by the development of localized dissolution channels, and 3) uniform dissolution characterized by a roughly homogeneous dissolution within the sample [8]. These distinct dissolution profiles are naturally associated with different evolutions of porosity and permeability. However, the impact of the rock microstructure on the development of alteration is not yet fully understood [11–13].

In order to provide a further understanding of the reactive transport processes and the resulting alteration of the rock, we combine experimental techniques with direct numerical simulations and a homogenization approach. Indeed, numerical simulation allows the determination of the velocity field in the porous structure, providing valuable information on the physics at the pore scale. Strictly speaking, in order to obtain the complete velocity field of a bimodal porous medium, it is necessary to solve the Stokes equation in the entire domain. Nevertheless, beyond the associated computational cost, this approach can be defeated if the micropores are too small to be described through imaging techniques. We overcome this deficiency by solving a specific equation, the Darcy-Brinkman equation [14, 15], allowing the description of flow in a bimodal porous medium. This equation can be seen as an interpolation between the Stokes equation, valid in the free fluid, and the Darcy equation, applied in the microporous phase. The major advantage of the Darcy-Brinkman equation lies in its simplicity of implementation as there is only one domain, meaning that the same equation is solved at any point of the domain. Simulations are performed in 2D images allowing the visualization of the flow field in the porous structure. The 2D simulations allow a refined representation of the flow, but their computational cost is too high to consider a systematic calculation of permeability, and the 3D nature of this property should also be taken into account. To this end, we propose to use a three-dimensional homogenization approach that will allow us to derive the initial permeability and its evolution under the effect of weathering.

The objective of this article is threefold: In the first step, we perform reactive transport experiments and analyze the resulting dissolution patterns through the scanner and

microscanner imaging. Then, in order to gain a further understanding of these experiments we conduct small-scale simulations to provide information on the local velocity field. Finally, we use a homogenization method to predict the variation in effective permeability due to a modification of the porous structure induced by alteration. The article is structured as follows: In the first step, we present the experimental techniques, the equation to be solved the numerical method to perform small-scale simulations, and the homogenization method to compute the permeability. Then, we discuss the corresponding results. Further details on the numerical scheme can be found in the Appendix A.

# 2 Methods

## 2.1 Alteration experiments

### 2.1.1 Sample description

In order to focus on the impact of the microstructure, we choose to study a limestone composed almost entirely of calcite: Euville limestone [16]. This limestone is Oxfordian in age and comes from the quarry of Euville located in the eastern part of the Paris Basin (Lorraine region). It consists of a crinoidal grainstone (Fig. 1) mostly made of crinoid ossicles (more than 95%), a few debris of echinoids, and rare lamellibranch bioclasts. The main cement phase corresponds to an important development of syntaxial overgrowths around crinoid debris. This cement partially occluded the primary interparticle macroporosity. A secondary intra-particle patchy microporosity is also clearly observed in the crinoid ossicles. Finally, the micritic sediment may present a microporosity (possibly intercrystalline). Measurements carried out on samples from different blocks show that Euville limestone has a porosity between 12% and 20% and a permeability between 10 and 150 mD. Its bimodal pore size distribution is confirmed by both mercury porosimetry and nuclear magnetic resonance, which lead to similar volume fractions for macroporosity and microporosity of about 65% and 35% (±5%) respectively [16].

### 2.1.2 Experimental workflow

The studied sample has a diameter of 40 mm and a length of 100 mm. The same sequence of analyses is performed before and after alteration to precisely characterize the impact of the injection of CO₂-saturated water: scanner and microscanner imaging and measurement of porosity and permeability. Medical CT scanner imaging provides a 3D image of the whole sample with a resolution of 0.625 mm. A porosity profile along the sample length can moreover be inferred from a correlation between the radiological densities measured on successive 0.625-mm-thick slices and data corresponding to pure calcite samples. Microscanner imaging allows a significantly higher resolution (20.28 µm for a 40-mm diameter sample) but requires four series of acquisition with an overlapping to image the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Petrographic overview of Euville limestone (transmitted light plane-polarized photomicrographs). Porosity appears in blue. [Left]: Cr: Crinoid ossicle, Ec: Echinoid, La: Lamellibranch bioclast, [Right]: Mi: Micropore, Ma: Macropore.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Schematic diagram of coreflight device showing upstream CO₂-saturated water injection system and downstream effluent sampling system.</figcaption></figure>

entire sample. The effective porosity is determined using the three-weight method:

$$ \phi_e = \frac{m_{\text{sat}} - m_{\text{dry}}}{m_{\text{sat}} - m_{\text{imm}}} \quad (1) $$

where $m_{\text{dry}}$ is the mass of the dry sample, $m_{\text{sat}}$ is the mass of the sample saturated with a 25 g/l NaCl brine, and $m_{\text{imm}}$ is the mass of the sample immersed in the fluid in which it is saturated.

Water permeability is directly measured in the coreflight device used for alteration (Fig. 2). A constant confining pressure of 30 bar is applied during the whole experiment which is performed at room temperature. A back pressure regulator ensures a constant outlet pore pressure of 2 bar. During permeability measurements, several flow rates are applied at the inlet. Each flow rate $Q$ is maintained until

equilibration of the differential pressure $\Delta P$. The sample intrinsic permeability $K$ is then deduced from the slope of the $Q - \Delta P$ curve using Darcy's law:

$$ K = \mu \frac{hQ}{S\Delta P} \quad (2) $$

where $\mu$ is the fluid dynamic viscosity, $h$ is the sample length, and $S$ is its section. At least four flow rates are applied for each permeability measurement.

The reactive fluid used for the alteration phase is pure water initially degassed and then kept in equilibrium with CO₂ gas at a pressure of 2 bar. The resulting solution of pH between 3.8 and 4 is injected into the sample at a flow rate of 140 cm³/h for 5 days. These experimental conditions were chosen to generate significant dissolution. The monitoring of the differential pressure allows us to follow the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Petrographic overview of Euville limestone (transmitted light plane-polarized photomicrographs). Porosity appears in blue. [Left]: Cr: Crinoid ossicle, Ec: Echinoid, La: Lamellibranch bioclast, [Right]: Mi: Micropore, Ma: Macropore.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Schematic diagram of coreflight device showing upstream CO₂-saturated water injection system and downstream effluent sampling system.</figcaption></figure>

entire sample. The effective porosity is determined using the three-weight method:

$$ \phi_e = \frac{m_{\text{sat}} - m_{\text{dry}}}{m_{\text{sat}} - m_{\text{imm}}} \qquad (1) $$

where $m_{\text{dry}}$ is the mass of the dry sample, $m_{\text{sat}}$ is the mass of the sample saturated with a 25 g/l NaCl brine, and $m_{\text{imm}}$ is the mass of the sample immersed in the fluid in which it is saturated.

Water permeability is directly measured in the coreflight device used for alteration (Fig. 2). A constant confining pressure of 30 bar is applied during the whole experiment which is performed at room temperature. A back pressure regulator ensures a constant outlet pore pressure of 2 bar. During permeability measurements, several flow rates are applied at the inlet. Each flow rate $Q$ is maintained until

equilibration of the differential pressure $\Delta P$. The sample intrinsic permeability $K$ is then deduced from the slope of the $Q - \Delta P$ curve using Darcy's law:

$$ K = \mu \frac{hQ}{S\Delta P} \qquad (2) $$

where $\mu$ is the fluid dynamic viscosity, $h$ is the sample length, and $S$ is its section. At least four flow rates are applied for each permeability measurement.

The reactive fluid used for the alteration phase is pure water initially degassed and then kept in equilibrium with CO₂ gas at a pressure of 2 bar. The resulting solution of pH between 3.8 and 4 is injected into the sample at a flow rate of 140 cm³/h for 5 days. These experimental conditions were chosen to generate significant dissolution. The monitoring of the differential pressure allows us to follow the
evolution of the permeability during the alteration. Fluid samples are taken at the outlet to monitor dissolved species.

## 2.2 Direct numerical simulations in 2D model porous structure

### 2.2.1 2D porous medium generation: Porosity-controlled aggregation (PlugIm)

The objective here is to create a 2D model porous structure representative of limestones with bimodal porosity consisting of solid grains, macropores, and a microporous matrix. The generation of a realistic porous rock structure is based on the aggregation model from [17]. Rock grains are aggregated one by one according to the position probabilities of their center calculated using morphological operators [18]. The possible positions of the center of the next grain are obtained by morphological dilation. Concave positions are extracted from morphological openings. The adjustment of these positions through a compacity parameter $\alpha$ [17] allows to set a final targeted porosity. The model proposed here, presenting relatively large areas free of any grains, is modified after each expansion step with a hole-filling operation [18] performed before extracting the admissible positions. The different shapes and sizes of the rock grains are drawn from a library of 2500 grains imaged by tomography from a real rock sample. The matrix phase is computed by the geodesic dilation of points from the skeleton of the macroporous network where the pore radius is superior to 4 pixels, the size of the geodesic dilation being equal to the initial porous radius. Simulations of size 1000 × 1000 pixels are generated with 1000 grains, with compacity parameter $\alpha$ = 1.0. Further details for the model and its parameters can be found in [17]. Implementation of the model is available from [19].

Figure 3 shows the model porous media on which we perform simulations, composed of 48.1% solid grains, 19.6% macropores, and 32.3% microporous matrix. Grains are enclosed in the microporous matrix, whereas macropores are in contact with grains and matrix. The microporous phase forms a percolating network in contrast to macropores which might be large but not percolating. However, in order to best represent the heterogeneity of natural porous media at this scale, macropores are not uniformly distributed. Here, their density is higher in the bottom and top of the image and lower in the center.

### 2.2.2 Flow equation

Neglecting the gravity forces, the 2D velocity field in the model porous structure built with PlugIm is obtained by solving the Darcy-Brinkman equation [14, 15] given by

$$-\mu_{\text{eff}}\Delta\vec{u} + \frac{\mu}{k}\vec{u} = -\nabla P \quad (3)$$

where $\mu_{\text{eff}}$ is the effective Brinkman viscosity, $\mu$ the dynamic viscosity of the fluid, and $k$ the local permeability assumed isotropic. $\nabla P$ is the driving pressure gradient and $\vec{u}$ the local velocity. In the matrix $k$ takes the value of the permeability of the microporous phase, whereas it becomes extremely large in the macropores.

<figure><img src="image_4.png" /><figcaption>Figure 3. Model porous medium generated with PlugIm. Macropores: black, solid grains: white and microporous matrix: grey.</figcaption></figure>

From the velocity field the global permeability $K_i$, with $i = x, y$, of the 2D domain can be computed as:

$$K_i = \frac{\mu u_{i,D} L_i}{\Delta P_i} \quad (4)$$

where $u_{i,D}$ is the Darcy velocity, $L_i$ is the domain length and $\Delta P_i$ is the pressure difference in direction $i$.

### 2.2.3 Lattice-Boltzmann method for the resolution of the Darcy-Brinkman equation

In this work, flow is simulated by means of a Lattice-Boltzmann Two-Relaxation-Time scheme [20–22]. It employs two relaxation parameters characterizing the relaxation of particle distributions towards an equilibrium state. The first relaxation rate is related to a physical variable whereas the second one corresponds to a numerical parameter ensuring numerical stability. The approach combines the simplicity of implementation of the single relaxation rate (BGK) method with the numerical advantages (accuracy and stability) of the Multi Relaxation Times (MRT) scheme. After the specification of the geometry of the porous medium, the flow field of the carrier fluid is computed by means of a specific Lattice Boltzmann scheme for the resolution of the Darcy-Brinkman equation [23, 24]. A detailed description of the numerical scheme is given in the Appendix A.

## 2.3 Homogenization method

Effective medium modeling allows the calculation of the homogenized permeability of a 3D porous medium without solving the flow in the 3D pore network. The calculations are performed on a Representative Elementary Volume (REV). The macroscopic scale considers the REV as an
evolution of the permeability during the alteration. Fluid samples are taken at the outlet to monitor dissolved species.

## 2.2 Direct numerical simulations in 2D model porous structure

### 2.2.1 2D porous medium generation: Porosity-controlled aggregation (PlugIm)

The objective here is to create a 2D model porous structure representative of limestones with bimodal porosity consisting of solid grains, macropores, and a microporous matrix. The generation of a realistic porous rock structure is based on the aggregation model from [17]. Rock grains are aggregated one by one according to the position probabilities of their center calculated using morphological operators [18]. The possible positions of the center of the next grain are obtained by morphological dilation. Concave positions are extracted from morphological openings. The adjustment of these positions through a compacity parameter $\alpha$ [17] allows to set a final targeted porosity. The model proposed here, presenting relatively large areas free of any grains, is modified after each expansion step with a hole-filling operation [18] performed before extracting the admissible positions. The different shapes and sizes of the rock grains are drawn from a library of 2500 grains imaged by tomography from a real rock sample. The matrix phase is computed by the geodesic dilation of points from the skeleton of the macroporous network where the pore radius is superior to 4 pixels, the size of the geodesic dilation being equal to the initial porous radius. Simulations of size 1000 × 1000 pixels are generated with 1000 grains, with compacity parameter $\alpha$ = 1.0. Further details for the model and its parameters can be found in [17]. Implementation of the model is available from [19].

Figure 3 shows the model porous media on which we perform simulations, composed of 48.1% solid grains, 19.6% macropores, and 32.3% microporous matrix. Grains are enclosed in the microporous matrix, whereas macropores are in contact with grains and matrix. The microporous phase forms a percolating network in contrast to macropores which might be large but not percolating. However, in order to best represent the heterogeneity of natural porous media at this scale, macropores are not uniformly distributed. Here, their density is higher in the bottom and top of the image and lower in the center.

### 2.2.2 Flow equation

Neglecting the gravity forces, the 2D velocity field in the model porous structure built with PlugIm is obtained by solving the Darcy-Brinkman equation [14, 15] given by

$$-\mu_{\text{eff}}\Delta\vec{u} + \frac{\mu}{k}\vec{u} = -\nabla P \quad (3)$$

where $\mu_{\text{eff}}$ is the effective Brinkman viscosity, $\mu$ the dynamic viscosity of the fluid, and $k$ the local permeability assumed isotropic. $\nabla P$ is the driving pressure gradient and $\vec{u}$ the local velocity. In the matrix $k$ takes the value of the permeability of the microporous phase, whereas it becomes extremely large in the macropores.

<figure><img src="image_4.png" /><figcaption>Figure 3. Model porous medium generated with PlugIm. Macropores: black, solid grains: white and microporous matrix: grey.</figcaption></figure>

From the velocity field the global permeability $K_i$, with $i = x, y$, of the 2D domain can be computed as:

$$K_i = \frac{\mu u_{i,D} L_i}{\Delta P_i} \quad (4)$$

where $u_{i,D}$ is the Darcy velocity, $L_i$ is the domain length and $\Delta P_i$ is the pressure difference in direction $i$.

### 2.2.3 Lattice-Boltzmann method for the resolution of the Darcy-Brinkman equation

In this work, flow is simulated by means of a Lattice-Boltzmann Two-Relaxation-Time scheme [20–22]. It employs two relaxation parameters characterizing the relaxation of particle distributions towards an equilibrium state. The first relaxation rate is related to a physical variable whereas the second one corresponds to a numerical parameter ensuring numerical stability. The approach combines the simplicity of implementation of the single relaxation rate (BGK) method with the numerical advantages (accuracy and stability) of the Multi Relaxation Times (MRT) scheme. After the specification of the geometry of the porous medium, the flow field of the carrier fluid is computed by means of a specific Lattice Boltzmann scheme for the resolution of the Darcy-Brinkman equation [23, 24]. A detailed description of the numerical scheme is given in the Appendix A.

## 2.3 Homogenization method

Effective medium modeling allows the calculation of the homogenized permeability of a 3D porous medium without solving the flow in the 3D pore network. The calculations are performed on a Representative Elementary Volume (REV). The macroscopic scale considers the REV as an
elementary particle with effective properties and the micro- scopic scale refers to the different phases composing the REV (solid grains, matrix phase, and macropores here). An exhaustive theoretical development of effective medium calculation for permeability is presented in [25].

In the framework of periodic homogenization, the steady-state, laminar, single-phase, and incompressible flow of a Newtonian viscous fluid in a porous medium obeys to Darcy's law at the macroscopic scale:

$$
\vec{u}_D = - \frac{\mathbb{K}^{\text{hom}}}{\mu} \vec{\nabla} P \qquad (5)
$$

where $\vec{u}_D$ is the Darcy velocity, $\mu$ is the fluid dynamic viscosity, $\vec{\nabla} P$ the pressure gradient, and $\mathbb{K}^{\text{hom}}$ is the second-order tensor of permeability that depends only on the porous medium geometry. The objective of the homogenization process is to estimate $\mathbb{K}^{\text{hom}}$. To this end, the different components of the REV (solid grains, micro- porous matrix, macropores) are considered fictitious Darcy media [26], characterized by negligible permeability for the solid phase and very high permeability for the macropores. The same type of equation, namely Darcy's law, then applies at both the microscopic and macroscopic scales. The first step is to define the properties of each phase. As a first approach, we introduce the solid grains, the microporous matrix, and the macropores as isotropi- cally oriented ellipsoids with a high aspect ratio (fixed to 0.8 for all three phases) and we assume that the perme- ability of each phase is isotropic. The macroscopic perme- ability tensor will thus also be isotropic. The second step is to assemble the different phases with respect to their volume fractions and compute the homogenized perme- ability through the self-consistent scheme based on Eshelby's auxiliary problem. In this approach, the interac- tion between the different phases is accounted for by considering each phase as an inclusion embedded in the homogenized medium. The self-consistent scheme is particularly well-suited for the disordered microstructures typical of carbonate rocks and has the further advantage of highlighting percolation thresholds in the case of high permeability contrasts between the phases (in this case the microporous phase and the macropores) [27].

3 Results and discussion

3.1 Experimental results

The studied sample has an initial porosity of 12.9% and an initial permeability of 40 mD. After five days of alteration within the coreflood cell, the permeability has drastically increased up to about 3400 mD, and the porosity measured after alteration is 15.2%. Analysis of the fluid at the coreflood outlet confirms the dissolution of calcite. The differential pressure decreases sharply from the beginning of alteration, leading to a rapid increase in permeability. Figure 4 shows a 3D representation of the increase in poros- ity across the sample length obtained through image differ- ence between the altered and intact state together with

porosity profiles showing the average porosity per section before and after alteration. The sample face on the injection side shows various localized dissolution areas (see Fig. 5A). Dissolution then develops globally homogeneously in the core of the sample, with a preferential dissolution path localized close to the periphery which spreads at around one-third of the sample length before relocating near the periphery. The preferential dissolution path is directly visi- ble on the sample surface (Fig. 5).

To investigate more finely a microstructural changes during dissolution, microscanner has been used to image the studied sample in intact and altered states. Image comparison shows that preferential dissolution mainly propagates through more porous areas (Fig. 6), either asso- ciated with large macropores and/or highly microporous areas. After alteration, areas of preferential dissolution show numerous significantly enlarged macropores. The preferential dissolution path created therefore retains a porous structure and does not correspond to a wormhole- type channel.

3.2 Flow field

In this section we present two-dimensional simulation results of flow in the model porous medium shown in Figure 3. To this goal we solve the Darcy-Brinkman equa- tion to obtain the flow field.

Figure 7 shows the velocity field in the model porous medium. A sensitivity analysis has been performed to study the impact of the permeability of the microporous matrix on the flow field. The results showed (see Fig. A.1 in the Appendix A for the flow field obtained with $k_{\text{micro}} = 0.1$), that the influence of the value of $k_{\text{micro}}$ lies in the intensity of the fluid velocity within the matrix phase. The overall behavior and particularly the location of the percolating paths are controlled by the geometric repartition of the different phases and are not modified by a change in $k_{\text{micro}}$. As the objective of the simulations is to provide further comprehension of the structure-dependent flow field on reactive transport and not to match experimental data at this stage, the permeability of the matrix is arbitrarily set to $k_{\text{micro}} = 1$. In Figure 7 the scale of the velocity field was slightly limited in order to enhance the contrast between low and high velocity zones. A nearly percolating path of higher velocities can be observed in the bottom part of the image although macropores do not form a percolating network. In this region, percolation of high-velocity path- ways is favored by the vicinity of the macropores and the presence of the microporous matrix. Indeed, velocity increases strongly in the matrix situated between two macropores. In the upper part of the image, there are high-velocity zones that are interrupted by a zone of lower porosity in the center left of the image. In the center of the porous medium, velocity becomes low due to the small number of macropores.

Consequently, when injected into the model bimodal porous medium, a passive or reactive fluid thus propagates more rapidly along paths with a higher macropore fraction, but nevertheless invades the entire cross-section of the sample. In order to better understand the influence of the
elementary particle with effective properties and the microscopic scale refers to the different phases composing the REV (solid grains, matrix phase, and macropores here). An exhaustive theoretical development of effective medium calculation for permeability is presented in [25].

In the framework of periodic homogenization, the steady-state, laminar, single-phase, and incompressible flow of a Newtonian viscous fluid in a porous medium obeys to Darcy's law at the macroscopic scale:

$$ \vec{u}_D = -\frac{\mathbb{K}^{\text{hom}}}{\mu} \vec{\nabla} P \qquad (5) $$

where $\vec{u}_D$ is the Darcy velocity, $\mu$ is the fluid dynamic viscosity, $\vec{\nabla} P$ the pressure gradient, and $\mathbb{K}^{\text{hom}}$ is the second-order tensor of permeability that depends only on the porous medium geometry. The objective of the homogenization process is to estimate $\mathbb{K}^{\text{hom}}$. To this end, the different components of the REV (solid grains, microporous matrix, macropores) are considered fictitious Darcy media [26], characterized by negligible permeability for the solid phase and very high permeability for the macropores. The same type of equation, namely Darcy's law, then applies at both the microscopic and macroscopic scales. The first step is to define the properties of each phase. As a first approach, we introduce the solid grains, the microporous matrix, and the macropores as isotropically oriented ellipsoids with a high aspect ratio (fixed to 0.8 for all three phases) and we assume that the permeability of each phase is isotropic. The macroscopic permeability tensor will thus also be isotropic. The second step is to assemble the different phases with respect to their volume fractions and compute the homogenized permeability through the self-consistent scheme based on Eshelby's auxiliary problem. In this approach, the interaction between the different phases is accounted for by considering each phase as an inclusion embedded in the homogenized medium. The self-consistent scheme is particularly well-suited for the disordered microstructures typical of carbonate rocks and has the further advantage of highlighting percolation thresholds in the case of high permeability contrasts between the phases (in this case the microporous phase and the macropores) [27].

# 3 Results and discussion

## 3.1 Experimental results

The studied sample has an initial porosity of 12.9% and an initial permeability of 40 mD. After five days of alteration within the coreflood cell, the permeability has drastically increased up to about 3400 mD, and the porosity measured after alteration is 15.2%. Analysis of the fluid at the coreflood outlet confirms the dissolution of calcite. The differential pressure decreases sharply from the beginning of alteration, leading to a rapid increase in permeability. Figure 4 shows a 3D representation of the increase in porosity across the sample length obtained through image difference between the altered and intact state together with

porosity profiles showing the average porosity per section before and after alteration. The sample face on the injection side shows various localized dissolution areas (see Fig. 5A). Dissolution then develops globally homogeneously in the core of the sample, with a preferential dissolution path localized close to the periphery which spreads at around one-third of the sample length before relocating near the periphery. The preferential dissolution path is directly visible on the sample surface (Fig. 5).

To investigate more finely a microstructural changes during dissolution, microscanner has been used to image the studied sample in intact and altered states. Image comparison shows that preferential dissolution mainly propagates through more porous areas (Fig. 6), either associated with large macropores and/or highly microporous areas. After alteration, areas of preferential dissolution show numerous significantly enlarged macropores. The preferential dissolution path created therefore retains a porous structure and does not correspond to a wormhole-type channel.

## 3.2 Flow field

In this section we present two-dimensional simulation results of flow in the model porous medium shown in Figure 3. To this goal we solve the Darcy-Brinkman equation to obtain the flow field.

Figure 7 shows the velocity field in the model porous medium. A sensitivity analysis has been performed to study the impact of the permeability of the microporous matrix on the flow field. The results showed (see Fig. A.1 in the Appendix A for the flow field obtained with $k_{\text{micro}} = 0.1$), that the influence of the value of $k_{\text{micro}}$ lies in the intensity of the fluid velocity within the matrix phase. The overall behavior and particularly the location of the percolating paths are controlled by the geometric repartition of the different phases and are not modified by a change in $k_{\text{micro}}$. As the objective of the simulations is to provide further comprehension of the structure-dependent flow field on reactive transport and not to match experimental data at this stage, the permeability of the matrix is arbitrarily set to $k_{\text{micro}} = 1$. In Figure 7 the scale of the velocity field was slightly limited in order to enhance the contrast between low and high velocity zones. A nearly percolating path of higher velocities can be observed in the bottom part of the image although macropores do not form a percolating network. In this region, percolation of high-velocity pathways is favored by the vicinity of the macropores and the presence of the microporous matrix. Indeed, velocity increases strongly in the matrix situated between two macropores. In the upper part of the image, there are high-velocity zones that are interrupted by a zone of lower porosity in the center left of the image. In the center of the porous medium, velocity becomes low due to the small number of macropores.

Consequently, when injected into the model bimodal porous medium, a passive or reactive fluid thus propagates more rapidly along paths with a higher macropore fraction, but nevertheless invades the entire cross-section of the sample. In order to better understand the influence of the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** [Top]: Alteration-induced porosity variations deduced from CT-scanner image difference, longitudinal sections at different rotation angles (0°, 45°, 90°, 135°, 180°), [Bottom]: CT-scanner porosity profile before (solid line) and after (dashed line) alteration.</figcaption></figure>

alteration on the flow field, we selected a sub-volume along the highest velocity path observed in the bottom part of the model (SV₇, see Fig. 9), with solid, matrix, and macropore fractions as close as possible to those of the Euville sample on which the alteration experiment was performed. As observed on microscanner images, preferentially altered areas show enlarged macropores and dissolved microporous material. Thus, we applied a dilation algorithm to macropores of the selected sub-volume and dilated the latter by 3, 6, and 9 pixels. Only the matrix was modified, supposing that the solid phase is not altered (the different images can be seen in Fig. 8 [Top]). Table 1 [Top] shows the fractions of the solid, matrix, and macropore phase before and after alteration of 3, 6, and 9 pixels. As can be seen, the solid phase remains constant, whereas the macropore fraction increases and the matrix fraction decreases.

We then computed the velocity field for each image (see Fig. 8). By increasing the size of the altered macropores a

path of higher velocity appears, which would further facilitate transport of the reactive fluid.

## 3.3 Permeability

In this section, we present the computation results of the effective permeability obtained by direct 2D simulations of the velocity field using Lattice Boltzmann (LB) and by applying 3D Effective Medium modeling (EM). Effective permeability was calculated for the total model porous medium and different sub-volumes (see Fig. 9).

The model displayed in Figure 3 is considered a macroscopic REV. Then, we split this model into three horizontal bands to separate the central lower velocity zone and the top and bottom higher velocity zones. The top and bottom bands are further divided into several sub-volumes showing rather homogeneous phase distribution but different fractions of macropores and microporous matrix (Tab. 2).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** [Top]: Alteration-induced porosity variations deduced from CT-scanner image difference, longitudinal sections at different rotation angles (0°, 45°, 90°, 135°, 180°), [Bottom]: CT-scanner porosity profile before (solid line) and after (dashed line) alteration.</figcaption></figure>

alteration on the flow field, we selected a sub-volume along the highest velocity path observed in the bottom part of the model (SV₇, see Fig. 9), with solid, matrix, and macropore fractions as close as possible to those of the Euville sample on which the alteration experiment was performed. As observed on microscanner images, preferentially altered areas show enlarged macropores and dissolved microporous material. Thus, we applied a dilation algorithm to macropores of the selected sub-volume and dilated the latter by 3, 6, and 9 pixels. Only the matrix was modified, supposing that the solid phase is not altered (the different images can be seen in Fig. 8 [Top]). Table 1 [Top] shows the fractions of the solid, matrix, and macropore phase before and after alteration of 3, 6, and 9 pixels. As can be seen, the solid phase remains constant, whereas the macropore fraction increases and the matrix fraction decreases.

We then computed the velocity field for each image (see Fig. 8). By increasing the size of the altered macropores a

path of higher velocity appears, which would further facilitate transport of the reactive fluid.

## 3.3 Permeability

In this section, we present the computation results of the effective permeability obtained by direct 2D simulations of the velocity field using Lattice Boltzmann (LB) and by applying 3D Effective Medium modeling (EM). Effective permeability was calculated for the total model porous medium and different sub-volumes (see Fig. 9).

The model displayed in Figure 3 is considered a macroscopic REV. Then, we split this model into three horizontal bands to separate the central lower velocity zone and the top and bottom higher velocity zones. The top and bottom bands are further divided into several sub-volumes showing rather homogeneous phase distribution but different fractions of macropores and microporous matrix (Tab. 2).
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Pictures of the sample after alteration (40-mm diameter, 100-mm length). A: Injection face, B: Outlet face, C: Longitudinal view.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Velocity field, using limited scale (Lattice-Boltzmann units) in order to enhance the contrast between low and high-velocity zones, $k_{\text{micro}} = 1$ (no flow on top and bottom faces, application of a pressure gradient between left and right faces).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison of microscanner images before and after alteration. [Left]: Position of the microscanner-imaged cross-sections, [Middle]: Microscanner images before alteration, [Right]: Microscanner images after alteration.</figcaption></figure>

Note that sub-volume $SV_7$ in Figure 9 corresponds to the volume altered in Section 3.2.

LB simulations in x- and y-direction can be performed on the images of the sub-volumes, thus providing 2D directional permeability estimates. The homogenization method is applied to each sub-volume considered here as an

independent REV to provide an effective 3D isotropic permeability according to the fractions of solid grains, macropores, and microporous matrix. Figure 10 presents the effective permeability of each sub-volume normalized by the permeability of the global REV (total image). Following the same rationale as for the LB simulations, the permeability of the microporous phase is set to 1 for homogenization computations. For both methods, permeability increases with increasing macropore fraction. However, for large macropore fractions, normalized permeability values obtained by EM are higher than those obtained by LB. This can be explained by the fact that the homogenization algorithm supposes a 3D structure based on the solid, matrix, and macropore fractions of the 2D images used for the LB velocity field computations. And, the percolation threshold of the macropore network in 3D is located at lower macropore fractions than for 2D structures [28, 29]. In general, close to the percolation threshold permeability increases more rapidly. In our case, even though the macropore fraction of certain sub-volumes is relatively high it is still below the percolation threshold of both LB and EM. Nevertheless, the fact that certain sub-volumes are close to the percolation threshold explains the increase in permeability and the difference observed between LB simulations and homogenization.

Table 1 [Bottom] shows the effective permeability obtained by LB simulations for the altered images shown in Figure 8 [Top]. As can be seen, alteration induces an increase in permeability. However, this increase is much lower than the experimentally observed one. This is due to the fact that no percolating path consisting of only macropores yet been created. This approach is too costly to be
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Pictures of the sample after alteration (40-mm diameter, 100-mm length). A: Injection face, B: Outlet face, C: Longitudinal view.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison of microscanner images before and after alteration. [Left]: Position of the microscanner-imaged cross-sections, [Middle]: Microscanner images before alteration, [Right]: Microscanner images after alteration.</figcaption></figure>

Note that sub-volume SV₇ in Figure 9 corresponds to the volume altered in Section 3.2.

LB simulations in x- and y-direction can be performed on the images of the sub-volumes, thus providing 2D directional permeability estimates. The homogenization method is applied to each sub-volume considered here as an

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Velocity field, using limited scale (Lattice-Boltzmann units) in order to enhance the contrast between low and high-velocity zones, $k_{\text{micro}} = 1$ (no flow on top and bottom faces, application of a pressure gradient between left and right faces).</figcaption></figure>

independent REV to provide an effective 3D isotropic permeability according to the fractions of solid grains, macropores, and microporous matrix. Figure 10 presents the effective permeability of each sub-volume normalized by the permeability of the global REV (total image). Following the same rationale as for the LB simulations, the permeability of the microporous phase is set to 1 for homogenization computations. For both methods, permeability increases with increasing macropore fraction. However, for large macropore fractions, normalized permeability values obtained by EM are higher than those obtained by LB. This can be explained by the fact that the homogenization algorithm supposes a 3D structure based on the solid, matrix, and macropore fractions of the 2D images used for the LB velocity field computations. And, the percolation threshold of the macropore network in 3D is located at lower macropore fractions than for 2D structures [28, 29]. In general, close to the percolation threshold permeability increases more rapidly. In our case, even though the macropore fraction of certain sub-volumes is relatively high it is still below the percolation threshold of both LB and EM. Nevertheless, the fact that certain sub-volumes are close to the percolation threshold explains the increase in permeability and the difference observed between LB simulations and homogenization.

Table 1 [Bottom] shows the effective permeability obtained by LB simulations for the altered images shown in Figure 8 [Top]. As can be seen, alteration induces an increase in permeability. However, this increase is much lower than the experimentally observed one. This is due to the fact that no percolating path consisting of only macropores yet been created. This approach is too costly to be
<figure><img src="image_9.png" /><figcaption>**Figure 8.** [Left]: Initial image of selected sub-volume ($SV_7$ in Fig. 9), [Right]: Altered images, geodesic dilation of macropores with 3 pixels, 6 pixels, 9 pixels. Macropores: black, solid grains: white and microporous matrix: grey. [Bottom Left]: Flow field in the initial image of selected sub-volume $SV_7$, [Bottom Right]: Flow fields in the altered images.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Different volumes considered for effective medium modeling and Lattice Boltzmann simulations. [Left]: Initial representative elementary volume (REV) with corresponding fractions of the different phases. [Right]: REV image divided into 9 sub-volumes. Each sub-volume is characterized by a specific fraction of macropores and microporous matrix.</figcaption></figure>

implemented on the whole REV. In contrast to that, effective medium modeling allows a straightforward prediction of the permeability of altered porous media as a function of the degree of alteration. The straight line in Figure 11 presents the permeability evolution predicted using effective medium modeling. A constant fraction of solid grains is used and fixed to one of the REV models (48.1%). Then, the macropore fraction is increased from 0% to 51.9%. Consequently, the microporous matrix fraction decreases from 51.9% to 0%. The red point marks the exact phase proportions of the global REV (total image). The EM effective permeability computed for each sub-volume is also reported. As their solid grain fractions are close to the one of the total volume, the corresponding blue points follow the overall permeability trend. The percolation threshold is obtained for a macropore fraction of about 33% (corresponding to a matrix fraction around 19%). Two flow regimes are highlighted: a matrix-controlled fluid flow for macropore fractions below

<table><caption>**Table 1.** [Top] Fractions of the solid grains, microporous matrix, and macropores before and after alteration of 3, 6 and 9 pixels. [Bottom] Normalized permeability obtained by LB of the altered sub-volume $SV_7$. $K_0$ corresponds to the effective permeability of sub-volume $SV_7$ before alteration.</caption><thead><tr><th>Dilation (in pixels)</th><th>0</th><th>3</th><th>6</th><th>9</th></tr></thead><tbody><tr><th>Macro</th><td>0.12</td><td>0.15</td><td>0.17</td><td>0.19</td></tr><tr><th>Matrix</th><td>0.36</td><td>0.33</td><td>0.31</td><td>0.29</td></tr><tr><th>Solid</th><td>0.51</td><td>0.51</td><td>0.51</td><td>0.51</td></tr><tr><th>K/K<sub>0</sub></th><td>1</td><td>1.47</td><td>1.66</td><td>1.96</td></tr></tbody></table>

<table><caption>**Table 2.** Fractions of macropores, microporous matrix, and solid in the different sub-volumes.</caption><thead><tr><th>(%)</th><th>Total (REV)</th><th>SV<sub>1</sub></th><th>SV<sub>2</sub></th><th>SV<sub>3</sub></th><th>SV<sub>4</sub></th></tr></thead><tbody><tr><th>Macro</th><td>19.6</td><td>27.7</td><td>11.3</td><td>27.9</td><td>0.5</td></tr><tr><th>Matrix</th><td>32.3</td><td>21.9</td><td>38.9</td><td>27.7</td><td>54.2</td></tr><tr><th>Solid</th><td>48.1</td><td>50.4</td><td>49.8</td><td>44.4</td><td>45.3</td></tr><tr><th>(%)</th><th>SV<sub>5</sub></th><th>SV<sub>6</sub></th><th>SV<sub>7</sub></th><th>SV<sub>8</sub></th><th>SV<sub>9</sub></th></tr><tr><th>Macro</th><td>31.2</td><td>30.8</td><td>12.2</td><td>13.1</td><td>21.8</td></tr><tr><th>Matrix</th><td>25.8</td><td>25.6</td><td>37.6</td><td>42.1</td><td>33.9</td></tr><tr><th>Solid</th><td>43.0</td><td>43.6</td><td>50.2</td><td>44.8</td><td>44.3</td></tr></tbody></table>

the percolation threshold and a fluid flow dominated by macropores above this threshold. Permeability values of the numerically altered sub-volume $SV_7$ clearly remain in the regime of the matrix-controlled fluid flow. If we now
<figure><img src="image_9.png" /><figcaption>**Figure 8.** [Left]: Initial image of selected sub-volume ($SV_7$ in Fig. 9), [Right]: Altered images, geodesic dilation of macropores with 3 pixels, 6 pixels, 9 pixels. Macropores: black, solid grains: white and microporous matrix: grey. [Bottom Left]: Flow field in the initial image of selected sub-volume $SV_7$, [Bottom Right]: Flow fields in the altered images.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Different volumes considered for effective medium modeling and Lattice Boltzmann simulations. [Left]: Initial representative elementary volume (REV) with corresponding fractions of the different phases. [Right]: REV image divided into 9 sub-volumes. Each sub-volume is characterized by a specific fraction of macropores and microporous matrix.</figcaption></figure>

implemented on the whole REV. In contrast to that, effective medium modeling allows a straightforward prediction of the permeability of altered porous media as a function of the degree of alteration. The straight line in Figure 11 presents the permeability evolution predicted using effective medium modeling. A constant fraction of solid grains is used and fixed to one of the REV models (48.1%). Then, the macropore fraction is increased from 0% to 51.9%. Consequently, the microporous matrix fraction decreases from 51.9% to 0%. The red point marks the exact phase proportions of the global REV (total image). The EM effective permeability computed for each sub-volume is also reported. As their solid grain fractions are close to the one of the total volume, the corresponding blue points follow the overall permeability trend. The percolation threshold is obtained for a macropore fraction of about 33% (corresponding to a matrix fraction around 19%). Two flow regimes are highlighted: a matrix-controlled fluid flow for macropore fractions below

<table><caption>**Table 1.** [Top] Fractions of the solid grains, microporous matrix, and macropores before and after alteration of 3, 6 and 9 pixels. [Bottom] Normalized permeability obtained by LB of the altered sub-volume $SV_7$. $K_0$ corresponds to the effective permeability of sub-volume $SV_7$ before alteration.</caption><thead><tr><th>Dilation (in pixels)</th><th>0</th><th>3</th><th>6</th><th>9</th></tr></thead><tbody><tr><th>Macro</th><td>0.12</td><td>0.15</td><td>0.17</td><td>0.19</td></tr><tr><th>Matrix</th><td>0.36</td><td>0.33</td><td>0.31</td><td>0.29</td></tr><tr><th>Solid</th><td>0.51</td><td>0.51</td><td>0.51</td><td>0.51</td></tr><tr><th>K/K<sub>0</sub></th><td>1</td><td>1.47</td><td>1.66</td><td>1.96</td></tr></tbody></table>

<table><caption>**Table 2.** Fractions of macropores, microporous matrix, and solid in the different sub-volumes.</caption><thead><tr><th>(%)</th><th>Total (REV)</th><th>SV<sub>1</sub></th><th>SV<sub>2</sub></th><th>SV<sub>3</sub></th><th>SV<sub>4</sub></th></tr></thead><tbody><tr><th>Macro</th><td>19.6</td><td>27.7</td><td>11.3</td><td>27.9</td><td>0.5</td></tr><tr><th>Matrix</th><td>32.3</td><td>21.9</td><td>38.9</td><td>27.7</td><td>54.2</td></tr><tr><th>Solid</th><td>48.1</td><td>50.4</td><td>49.8</td><td>44.4</td><td>45.3</td></tr><tr><th>(%)</th><th>SV<sub>5</sub></th><th>SV<sub>6</sub></th><th>SV<sub>7</sub></th><th>SV<sub>8</sub></th><th>SV<sub>9</sub></th></tr><tr><th>Macro</th><td>31.2</td><td>30.8</td><td>12.2</td><td>13.1</td><td>21.8</td></tr><tr><th>Matrix</th><td>25.8</td><td>25.6</td><td>37.6</td><td>42.1</td><td>33.9</td></tr><tr><th>Solid</th><td>43.0</td><td>43.6</td><td>50.2</td><td>44.8</td><td>44.3</td></tr></tbody></table>

the percolation threshold and a fluid flow dominated by macropores above this threshold. Permeability values of the numerically altered sub-volume $SV_7$ clearly remain in the regime of the matrix-controlled fluid flow. If we now
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Effective permeability calculated for the different sub-volumes using Effective Medium modeling (EM) and LB. Permeability is normalized by the effective permeability of the REV (total image). For LB simulations permeability was calculated in the x- and y-direction.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Normalized permeability as a function of the macropore fraction. The straight line was obtained by means of the effective medium approach. Blue dots correspond to the permeability computed by effective medium modeling of the sub-volumes and green triangles indicate permeability changes due to alteration obtained by LB. All values are normalized by the permeability value of the total volume ($K_{\text{total}}$).</figcaption></figure>

return to the experimental data, the total increase in porosity induced by alteration is equal to 2.3% (porosity unit) and thus remains limited. According to the modeling approach, the associated 85-fold increase in permeability is only compatible with macropore fractions approaching the percolation threshold for consecutive sub-volumes of the sample. And this is indeed what can be seen on the microscanner images after alteration in Figure 6.

# 4 Conclusion

In this article we presented different approaches to provide further understanding of the physics of reactive

transport within carbonate rocks. We first conducted injection experiments with CO₂-saturated water in a sample of Euville limestone combined with medical CT scanner and microscanner imaging before and after alteration. This limestone has a bimodal pore size distribution. Experimental results showed that a preferential dissolution path characterized by a porous structure with enlarged macropores was created resulting in an increase in permeability of nearly two decades. To further investigate small scale flow processes involved in reactive transport, we performed Lattice Boltzmann simulations in a numerically-generated model 2D porous medium. To approach the structure of the Euville limestone the generated porous medium consisted of solid grains, microporous phase, and macropores. In some parts of this artificial porous structure, macropores were nearly connected, whereas in other parts they were isolated, but connected through the matrix. We showed that the fluid velocity increased in the nearly percolating paths of macropores. Velocity was also enhanced in the microporous matrix situated between close macropores. In high-velocity zones CO₂-saturated water input starts earlier than in low-velocity zones leading probably to an earlier onset of dissolution that results in a stronger local alteration. When numerically altering the size of the macropores we could show that the velocity increased particularly in the altered zones, but not as much as it was experimentally observed. We therefore used a homogenization approach allowing a straightforward prediction of the permeability evolution of altered porous media as a function of the degree of alteration. Effective medium modeling showed that as long as the alteration does not create percolating macroporosity, permeability variation remains low as has been observed with the local Lattice-Boltzmann simulations. However, once macropores approach percolation, the increase in permeability reaches several orders of magnitude. In conclusion, by applying different experimental, Lattice-Boltzmann, and homogenization approaches we were able to provide further understanding of the process of dissolution of porous media due to CO₂-saturated water injection. In future work, it might be interesting to apply the modeling approaches to more representative microstructure models derived from imaging techniques and studying other limestones with different microstructures, in particular ones displaying wormhole development.

## Acknowledgments

The authors would like to acknowledge Leo Agelas, Gabriel Farag, and Isa- belle Faille for implementing the Darcy-Brinkmann equation in the open- source Lattice-Boltzmann platform Walberla. They would also like to thank Souhail Youssef for fruitful discussions on the subject.

## Authors contribution statement

Daniela Bauer, Elisabeth Bemer, Théo Briolet, Mathilde Adelinet, and Maxime Moreaud, wrote the article. Elisabeth Bemer, Théo Briolet, Olivier Sissmann, and Jérôme Fortin, designed the experimental workflow and analyzed the data. Théo Briolet carried out the alteration experiment. Maxime Pelerin performed the scanner and micro-scanner analyses. Maxime Moreaud created the 2D model porous structure. Daniela Bauer performed the Lattice-Boltzmann simulations and Mathilde Adelinet the homogenization modeling.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Effective permeability calculated for the different sub-volumes using Effective Medium modeling (EM) and LB. Permeability is normalized by the effective permeability of the REV (total image). For LB simulations permeability was calculated in the x- and y-direction.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Normalized permeability as a function of the macropore fraction. The straight line was obtained by means of the effective medium approach. Blue dots correspond to the permeability computed by effective medium modeling of the sub-volumes and green triangles indicate permeability changes due to alteration obtained by LB. All values are normalized by the permeability value of the total volume ($K_{\text{total}}$).</figcaption></figure>

return to the experimental data, the total increase in porosity induced by alteration is equal to 2.3% (porosity unit) and thus remains limited. According to the modeling approach, the associated 85-fold increase in permeability is only compatible with macropore fractions approaching the percolation threshold for consecutive sub-volumes of the sample. And this is indeed what can be seen on the microscanner images after alteration in Figure 6.

## 4 Conclusion

In this article we presented different approaches to provide further understanding of the physics of reactive

transport within carbonate rocks. We first conducted injection experiments with CO₂-saturated water in a sample of Euville limestone combined with medical CT scanner and microscanner imaging before and after alteration. This limestone has a bimodal pore size distribution. Experimental results showed that a preferential dissolution path characterized by a porous structure with enlarged macropores was created resulting in an increase in permeability of nearly two decades. To further investigate small scale flow processes involved in reactive transport, we performed Lattice Boltzmann simulations in a numerically-generated model 2D porous medium. To approach the structure of the Euville limestone the generated porous medium consisted of solid grains, microporous phase, and macropores. In some parts of this artificial porous structure, macropores were nearly connected, whereas in other parts they were isolated, but connected through the matrix. We showed that the fluid velocity increased in the nearly percolating paths of macropores. Velocity was also enhanced in the microporous matrix situated between close macropores. In high-velocity zones CO₂-saturated water input starts earlier than in low-velocity zones leading probably to an earlier onset of dissolution that results in a stronger local alteration. When numerically altering the size of the macropores we could show that the velocity increased particularly in the altered zones, but not as much as it was experimentally observed. We therefore used a homogenization approach allowing a straightforward prediction of the permeability evolution of altered porous media as a function of the degree of alteration. Effective medium modeling showed that as long as the alteration does not create percolating macroporosity, permeability variation remains low as has been observed with the local Lattice-Boltzmann simulations. However, once macropores approach percolation, the increase in permeability reaches several orders of magnitude. In conclusion, by applying different experimental, Lattice-Boltzmann, and homogenization approaches we were able to provide further understanding of the process of dissolution of porous media due to CO₂-saturated water injection. In future work, it might be interesting to apply the modeling approaches to more representative microstructure models derived from imaging techniques and studying other limestones with different microstructures, in particular ones displaying wormhole development.

## Acknowledgments

The authors would like to acknowledge Leo Agelas, Gabriel Farag, and Isa- belle Faille for implementing the Darcy-Brinkmann equation in the open- source Lattice-Boltzmann platform Walberla. They would also like to thank Souhail Youssef for fruitful discussions on the subject.

## Authors contribution statement

Daniela Bauer, Elisabeth Bemer, Théo Briolet, Mathilde Adelinet, and Maxime Moreaud, wrote the article. Elisabeth Bemer, Théo Briolet, Olivier Sissmann, and Jérôme Fortin, designed the experimental workflow and analyzed the data. Théo Briolet carried out the alteration experiment. Maxime Pelerin performed the scanner and micro-scanner analyses. Maxime Moreaud created the 2D model porous structure. Daniela Bauer performed the Lattice-Boltzmann simulations and Mathilde Adelinet the homogenization modeling.
## References

1 Pörtner H., Roberts D., Tignor M., Poloczanska E., Mintenbeck K., Alegria A., Craig M., Langsdorf S., Löschke S., Möller V., Okem A., Rama B. (2022) *IPCC 2022: Climate change 2022: Impacts, adaptation, and vulnerability, contribution of working group II to the sixth assessment report of the intergovernmental panel on climate change*, Cambridge University Press.

2 Lombard J.M., Azaroual M., Pironon J., Broseta D., Egermann P., Munier G., Mouronval G. (2010) CO₂ injectivity in geological storages: An overview of program and results of the GeoCarbone-Injectivity project, *Oil Gas Sci. Technol.* 65, 4, 533-539.

3 Siqueira T.A., Iglesias R.S., Ketzer J.M. (2017) Carbon dioxide injection in carbonate reservoirs – a review of CO₂-water-rock interaction studies, *Greenh. Gases: Sci. Technol.* 7, 5, 802–816.

4 André L., Audigane P., Azaroual M., Menjoz A. (2007) Numerical modeling of fluid-rock chemical interactions at the supercritical CO₂-liquid interface during CO₂ injection into a carbonate reservoir, the Dogger aquifer (Paris Basin, France), *Energy Convers. Manag.* 48, 6, 1782–1797.

5 Lønøy A. (2006) Making sense of carbonate pore systems, *AAPG Bull.* 90, 9, 1381–1405.

6 Bauer D., Youssef S., Han M., Bekri S., Rosenberg E., Fleury M., Vizika O. (2011) From computed microtomography images to resistivity index calculations of heterogeneous carbonates using a dual-porosity pore-network approach: Influence of percolation on the electrical transport properties, *Phys. Rev. E* 84, 1, 011133.

7 Bazin B. (2001) From matrix acidizing to acid fracturing: a laboratory evaluation of acid/rock interactions, *SPE Prod. Facil.* 16, 01, 22-29.

8 Golfier F., Zarcone C., Bazin B., Lenormand R., Lasseux D., Quintard M. (2002) On the ability of a darcy-scale model to capture wormhole formation during the dissolution of a porous medium, *J. Fluid Mech.* 457, 213–254.

9 Vialle S., Contraires S., Zinzsner B., Clavaud J.-B., Mahiouz K., Zuddas P., Zamora M. (2014) Percolation of CO₂-rich fluids in a limestone sample: Evolution of hydraulic, electrical, chemical, and structural properties, *J. Geophys. Res. Solid Earth* 119, 4, 2828–2847.

10 Menke H.P., Bijeljic B., Andrew M.G., Blunt M.J. (2015) Dynamic three-dimensional pore-scale imaging of reaction in a carbonate at reservoir conditions, *Environ. Sci. Technol.* 49, 7, 4407–4414.

11 Vialle S., Dvorkin J., Mavko G. (2013) Implications of pore microgeometry heterogeneity for the movement and chemical reactivity of CO₂ in carbonates, *Geophysics* 78, 5, L69–L86.

12 Yang Y., Li Y., Yao J., Iglauer S., Luquot L., Zhang K., Sun H., Zhang L., Song W., Wang Z. (2020) Dynamic pore-scale dissolution by CO₂-saturated brine in carbonates: Impact of homogeneous versus fractured versus vuggy pore structure, *Water Resour. Res.* 56, 4, e2019WR026112.

13 Peter A., Yang D., Eshiet K.I.I.I., Sheng Y. (2022) A review of the studies on CO₂-brine-rock interaction in geological storage process, *Geosciences* 12, 4, 168.

14 Brinkman H.C. (1949) A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles, *Flow, Turbul. Combust.* 1, 1, 27.

15 Brinkman H.C. (1949) On the permeability of media consisting of closely packed porous particles, *Flow Turbul. Combust.* 1, 1, 81.

16 Bemer E., Nguyen M.T., Dautriet J., Adelinet M., Fleury M., Youssef S. (2016) Impact of chemical alteration on the poromechanical properties of carbonate rocks, *Geophys. Prospect.* 64, 4-Advances in Rock Physics, 810–827.

17 Ferri G., Humbert S., Digne M., Schweitzer J.-M., Moreaud M. (2021) Simulation of large aggregate particles systems with a new morphological model, *Image Anal. Stereol.* 40, 71–84.

18 Serra J. (1988) *Image analysis and mathematical morphology, part II: Theoretical advances.*

19 “plugim!” (2018) *An open access and customizable software for signal and image processing*. https://www.plugim.fr.

20 Ginzburg I. (2008) Consistent Lattice Boltzmann schemes for the brinkman model of porous flow and infinite Chapman-Enskog expansion, *Phys. Rev. E* 77, (6), 066704.

21 Ginzburg I. (2007) Lattice Boltzmann modeling with discontinuous collision components: Hydrodynamic and advection-diffusion equations, *J. Stat. Phys.* 126, 157–206.

22 Ginzburg I., d’Humières D., Kuzmin A. (2010) Optimal stability of advection-diffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium, *J. Stat. Phys.* 139, 6, 1090–1143.

23 Ginzburg I., Silva G., Talon L. (2015) Analysis and improvement of Brinkman lattice Boltzmann schemes: Bulk, boundary, interface. Similarity and distinctness with finite elements in heterogeneous porous media, *Phys. Rev. E* 91, 2, 023307.

24 Silva G., Ginzburg I. (2016) Stokes–Brinkman–Darcy solutions of bimodal porous flow across periodic array of permeable cylindrical inclusions: Cell model, lubrication theory and LBM/FEM numerical simulations, *Transp. Porous Media* 111, 3, 795–825.

25 Barthélémy J.-F. (2009) Effective permeability of media with a dense network of long and micro fractures, *Transp. Porous Media* 76, 1, 153–178.

26 Fokker P.A. (2001) General anisotropic effective medium theory for the effective permeability of heterogeneous reservoirs, *Transp. Porous Media* 44, 2, 205–218.

27 Dormieux L., Kondo D. (2004) Approche micromécanique du couplage perméabilité-endommagement, *CR Mécanique* 332, 2, 135–140.

28 Stauffer D., Aharony A. (1992) *Introduction to percolation theory*, Taylor and Francis, London.

29 Masihi M., Gago P.A., King P.R. (2016) Estimation of the effective permeability of heterogeneous porous media by using percolation concepts, *Transp. Porous Media* 114, 169–199.

30 Talon L., Bauer D., Gland N., Youssef S., Auradou H., Ginzburg I. (2012) Assessment of the two relaxation time Lattice-Boltzmann scheme to simulate Stokes flow in porous media, *Water Resour. Res.* 48, 4.

# Appendix A

## Lattice Boltzmann schemes

The basis of the Lattice Boltzmann Method (LBM) lies in the specific discretization of the Boltzmann equation. Particle velocities and space are discretized on a regular grid, where nodes are related by the velocity vectors $\vec{v}_q$, $q \in [0, q_m]$, $q_m = 8$. Particles are only displaced with velocities $\vec{v}_q$. Populations $f_q(\vec{r}, t)$ are defined as the density of particles at position $\vec{r}$ with velocity $\vec{v}_q$. When implementing the Lattice Boltzmann model, the computation of the temporal evolution of the particle density is split in two distinct steps, the propagation and collision step, that are processed individually.

**Collision step:** In this step, distribution functions $f_q$ evolve according to the collision operator $\Omega$, specific for the equation to be solved. Distribution functions $f_q$ are updated as follows:

$$ \tilde{f}_q(\vec{r}, t) = f_q(\vec{r}, t) + \Omega[f(\vec{r}, t)]_q \quad (A1) $$

where $\tilde{f}_q(\vec{r}, t)$ stands for the post-collision distributions.

**Propagation step:** Here, distribution functions $f_q$ are exchanged between neighboring nodes with respect to the velocity set:

$$ f_q(\vec{r} + \vec{v}_q \Delta t, t + \Delta t) = \tilde{f}_q(\vec{r}, t). \quad (A2) $$

In the Two Relaxation Time Lattice Boltzmann scheme (TRT) the collision operator contains two different relaxation rates, one for the symmetric and one for the antisymmetric components. Symmetric components are defined
## References

1 Pörtner H., Roberts D., Tignor M., Polozanska E., Mintenbeck K., Alegria A., Craig M., Langsdorf S., Löschke S., Möller V., Okem A., Rama B. (2022) *IPCC 2022: Climate change 2022: Impacts, adaptation, and vulnerability, contribution of working group II to the sixth assessment report of the intergovernmental panel on climate change*, Cambridge University Press.

2 Lombard J.M., Azaroual M., Pironon J., Broseta D., Egermann P., Munier G., Mouronval G. (2010) CO₂ injectivity in geological storages: An overview of program and results of the GeoCarbone-Injectivity project, *Oil Gas Sci. Technol.* 65, 4, 533–539.

3 Siqueira T.A., Iglesias R.S., Ketzer J.M. (2017) Carbon dioxide injection in carbonate reservoirs – a review of CO₂-water-rock interaction studies, *Greenh. Gases: Sci. Technol.* 7, 5, 802–816.

4 André L., Audigane P., Azaroual M., Menjoz A. (2007) Numerical modeling of fluid-rock chemical interactions at the supercritical CO₂-liquid interface during CO₂ injection into a carbonate reservoir, the Dogger aquifer (Paris Basin, France), *Energy Convers. Manag.* 48, 6, 1782–1797.

5 Lønøy A. (2006) Making sense of carbonate pore systems, *AAPG Bull.* 90, 9, 1381–1405.

6 Bauer D., Youssef S., Han M., Bekri S., Rosenberg E., Fleury M., Vizika O. (2011) From computed microtomography images to resistivity index calculations of heterogeneous carbonates using a dual-porosity pore-network approach: Influence of percolation on the electrical transport properties, *Phys. Rev. E* 84, 1, 011133.

7 Bazin B. (2001) From matrix acidizing to acid fracturing: a laboratory evaluation of acid/rock interactions, *SPE Prod. Facil.* 16, 01, 22–29.

8 Golfier F., Zarcone C., Bazin B., Lenormand R., Lasseux D., Quintard M. (2002) On the ability of a darcy-scale model to capture wormhole formation during the dissolution of a porous medium, *J. Fluid Mech.* 457, 213–254.

9 Vialle S., Contraires S., Zinzsner B., Clavaud J.-B., Mahiouz K., Zuddas P., Zamora M. (2014) Percolation of CO₂-rich fluids in a limestone sample: Evolution of hydraulic, electrical, chemical, and structural properties, *J. Geophys. Res. Solid Earth* 119, 4, 2828–2847.

10 Menke H.P., Bijeljic B., Andrew M.G., Blunt M.J. (2015) Dynamic three-dimensional pore-scale imaging of reaction in a carbonate at reservoir conditions, *Environ. Sci. Technol.* 49, 7, 4407–4414.

11 Vialle S., Dvorkin J., Mavko G. (2013) Implications of pore microgeometry heterogeneity for the movement and chemical reactivity of CO₂ in carbonates, *Geophysics* 78, 5, L69–L86.

12 Yang Y., Li Y., Yao J., Iglauer S., Luquot L., Zhang K., Sun H., Zhang L., Song W., Wang Z. (2020) Dynamic pore-scale dissolution by CO₂-saturated brine in carbonates: Impact of homogeneous versus fractured versus vuggy pore structure, *Water Resour. Res.* 56, 4, e2019WR026112.

13 Peter A., Yang D., Eshiet K.I.I.I., Sheng Y. (2022) A review of the studies on CO₂-brine-rock interaction in geological storage process, *Geosciences* 12, 4, 168.

14 Brinkman H.C. (1949) A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles, *Flow, Turbul. Combust.* 1, 1, 27.

15 Brinkman H.C. (1949) On the permeability of media consisting of closely packed porous particles, *Flow Turbul. Combust.* 1, 1, 81.

16 Bemer E., Nguyen M.T., Dautriet J., Adelinet M., Fleury M., Youssef S. (2016) Impact of chemical alteration on the poromechanical properties of carbonate rocks, *Geophys. Prospect.* 64, 4-Advances in Rock Physics, 810–827.

17 Ferri G., Humbert S., Digne M., Schweitzer J.-M., Moreaud M. (2021) Simulation of large aggregate particles systems with a new morphological model, *Image Anal. Stereol.* 40, 71–84.

18 Serra J. (1988) *Image analysis and mathematical morphology, part II: Theoretical advances.*

19 “plugim!” (2018) *An open access and customizable software for signal and image processing*. https://www.plugim.fr.

20 Ginzburg I. (2008) Consistent Lattice Boltzmann schemes for the brinkman model of porous flow and infinite Chapman-Enskog expansion, *Phys. Rev. E* 77, (6), 066704.

21 Ginzburg I. (2007) Lattice Boltzmann modeling with discontinuous collision components: Hydrodynamic and advection-diffusion equations, *J. Stat. Phys.* 126, 157–206.

22 Ginzburg I., d’Humières D., Kuzmin A. (2010) Optimal stability of advection-diffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium, *J. Stat. Phys.* 139, 6, 1090–1143.

23 Ginzburg I., Silva G., Talon L. (2015) Analysis and improvement of Brinkman lattice Boltzmann schemes: Bulk, boundary, interface. Similarity and distinctness with finite elements in heterogeneous porous media, *Phys. Rev. E* 91, 2, 023307.

24 Silva G., Ginzburg I. (2016) Stokes–Brinkman–Darcy solutions of bimodal porous flow across periodic array of permeable cylindrical inclusions: Cell model, lubrication theory and LBM/FEM numerical simulations, *Transp. Porous Media* 111, 3, 795–825.

25 Barthélémy J.-F. (2009) Effective permeability of media with a dense network of long and micro fractures, *Transp. Porous Media* 76, 1, 153–178.

26 Fokker P.A. (2001) General anisotropic effective medium theory for the effective permeability of heterogeneous reservoirs, *Transp. Porous Media* 44, 2, 205–218.

27 Dormieux L., Kondo D. (2004) Approche micromécanique du couplage perméabilité-endommagement, *CR Mécanique* 332, 2, 135–140.

28 Stauffer D., Aharony A. (1992) *Introduction to percolation theory*, Taylor and Francis, London.

29 Masihi M., Gago P.A., King P.R. (2016) Estimation of the effective permeability of heterogeneous porous media by using percolation concepts, *Transp. Porous Media* 114, 169–199.

30 Talon L., Bauer D., Gland N., Youssef S., Auradou H., Ginzburg I. (2012) Assessment of the two relaxation time Lattice-Boltzmann scheme to simulate Stokes flow in porous media, *Water Resour. Res.* 48, 4.

## Appendix A

### Lattice Boltzmann schemes

The basis of the Lattice Boltzmann Method (LBM) lies in the specific discretization of the Boltzmann equation. Particle velocities and space are discretized on a regular grid, where nodes are related by the velocity vectors $\vec{v}_q$, $q \in [0, q_m]$, $q_m = 8$. Particles are only displaced with velocities $\vec{v}_q$. Populations $f_q(\vec{r}, t)$ are defined as the density of particles at position $\vec{r}$ with velocity $\vec{v}_q$. When implementing the Lattice Boltzmann model, the computation of the temporal evolution of the particle density is split in two distinct steps, the propagation and collision step, that are processed individually.

**Collision step:** In this step, distribution functions $f_q$ evolve according to the collision operator $\Omega$, specific for the equation to be solved. Distribution functions $f_q$ are updated as follows:

$$ \tilde{f}_q(\vec{r}, t) = f_q(\vec{r}, t) + \Omega[f(\vec{r}, t)]_q \quad (A1) $$

where $\tilde{f}_q(\vec{r}, t)$ stands for the post-collision distributions.

**Propagation step:** Here, distribution functions $f_q$ are exchanged between neighboring nodes with respect to the velocity set:

$$ f_q(\vec{r} + \vec{v}_q \Delta t, t + \Delta t) = \tilde{f}_q(\vec{r}, t). \quad (A2) $$

In the Two Relaxation Time Lattice Boltzmann scheme (TRT) the collision operator contains two different relaxation rates, one for the symmetric and one for the antisymmetric components. Symmetric components are defined
as $f_q^+ = \frac{1}{2}(f_q + f_{\bar{q}})$ while the anti-symmetric components are defined as $f_q^- = \frac{1}{2}(f_q - f_{\bar{q}})$ for $q \in \{1, ..., q_m/2\}$. For $q = 0$, $f_0^+ = f_0$ and $f_0^- = 0$.

The two relaxation parameters are defined as follows: $\lambda^-$ for all anti-symmetric non-equilibrium components $n_q^- = f_q^- - e_q^-$ and $\lambda^+$ for all symmetric non-equilibrium components $n_q^+ = f_q^+ - e_q^+$. Equilibrium components for the zero velocity are given by: $e_0^+ = e_0$ and $e_0^- = 0$.

Equation (A1) is then expressed as [22]:

$$ \begin{aligned} \tilde{f}_q(\vec{r}, t) &= f_q(\vec{r}, t) + \lambda^+ n_q^+ + \lambda^- n_q^- \\ \tilde{f}_{\bar{q}}(\vec{r}, t) &= f_{\bar{q}}(\vec{r}, t) + \lambda^+ n_q^+ - \lambda^- n_q^- \\ \tilde{f}_0(\vec{r}, t) &= f_0(\vec{r}, t)(1 + \lambda^+) - \lambda^+ e_0 \end{aligned} \quad (A3) $$

with $q \in \{1, ..., q_m/2\}$. The LBM-TRT requires the definition of a set of numerical parameters $\Lambda^\pm$ and $\Lambda$ to ensure the stability of the algorithm and to reduce the numerical error. These parameters are related to the relaxation parameters $\lambda^\pm$ as follows [22]:

$$ \Lambda = \Lambda^+ \Lambda^- \quad (A4) $$

with

$$ \Lambda^{\pm} = -(1/2 + 1/\lambda^{\pm}) \quad (A5) $$

for $-2 < \lambda^\pm < 0$. Values of the latter parameters depend on the equation to be solved. Specificities of the resolution of Darcy-Brinkman equation are presented in the following.

### Darcy-Brinkman

Darcy-Brinkman equation, in the absence of gravity, is given by

$$ -\mu_{\text{eff}} \Delta \vec{u} + \frac{\mu}{k} \vec{u} = -\nabla P \quad (A6) $$

or by using $\vec{u} = \vec{j}/\rho_0$ and $v = \mu/\rho_0$ as

$$ -v_{\text{eff}} \Delta \vec{j} + \frac{v}{k} \vec{j} = -\nabla P \quad (A7) $$

where $v_{\text{eff}}$ is the effective Brinkman viscosity, $v$ the kinematic viscosity of the fluid and $k$ the permeability. $\rho_0$ is an average density fixed to $\rho_0 = 1$. Equation (A7) can be expressed as:

$$ \nabla P = \vec{F} + v_{\text{eff}} \Delta \vec{j} \quad (A8) $$

corresponding to the stationary Stokes equation where the force term $\vec{F}$ has the following form:

$$ \vec{F} = -\frac{v}{k} \vec{j} \quad (A9) $$

As can be seen $\vec{F}$ depends on the macroscopic momentum $\vec{j}$. Before each propagation step, $j$ is computed from [23]

$$ \vec{j} = \frac{\vec{J}}{1 + \frac{v}{2k}} \quad (A10) $$

<figure><img src="image_13.png" /><figcaption>Figure A.1. Flow field in the model porous medium for $k_{\text{micro}} = 0.1$.</figcaption></figure>

with $\vec{J} = \sum_{q=1}^{q_m} f_q \vec{v}_q$.

Equation (A10) is then used to compute the equilibrium functions defined as [21, 30]:

$$ \begin{aligned} e_q^+(\vec{r}, t) &= c_s^2 \rho(\vec{r}, t) t_q^* \\ e_q^-(\vec{r}, t) &= t_q^*(\vec{j}(\vec{r}, t) \cdot \vec{v}_q) \\ e_0^+(\vec{r}, t) &= \rho(\vec{r}, t) - \sum_{q=1}^{q_m} e_q^+(\vec{r}, t) \\ e_0^-(\vec{r}, t) &= \rho(\vec{r}, t) - \sum_{q=1}^{q_m} e_q^-(\vec{r}, t) \end{aligned} \quad (A11) $$

where $c_s$ is the sound velocity ($c_s^2 = 1/3$) and $t_q^*$ correspond to isotropic physical weights. $t_q^*$ take the value $t_q^* = \left\{ \frac{1}{3}; \frac{1}{12} \right\}$ for the first (coordinate) and the second (diagonal) neighbor link in the $d_2 q_9$ scheme [21].

The calculation of the equilibrium functions requires the computation of the fluid density $\rho$ which is expressed as [21]: $\rho(\vec{r}, t) = \sum_{q=0}^{q_m} f_q(\vec{r}, t)$. The pressure is defined as $P = c_s^2 \rho$ and the kinematic viscosity as $\nu = \frac{\Lambda^+}{3} = -\frac{1}{3} \left( \frac{1}{2} + \frac{1}{\lambda^+} \right)$. In our case, numerical accuracy is obtained by using $\Lambda = 3/16$ [20] and prescribing $\nu = 1/12$. $\lambda^+$ and $\lambda^-$ can then be computed from equations (A4) and (A5). At the solid/fluide interface (when $\vec{v}_q$ is pointing towards the wall) the bounce back rule is applied: $f_{\bar{q}}(\vec{r}, t + \Delta t) = \vec{f}_q(\vec{r}, t)$.

More details on this method can be found in [21, 23, 30].

### Sensitivity analysis

Figure A.1 shows the flow field in the model porous medium for $k_{\text{micro}} = 0.1$. As can be seen the overall behavior is identical to the one of $k_{\text{micro}} = 1.0$.
as $f_q^+ = \frac{1}{2}(f_q + f_{\bar{q}})$ while the anti-symmetric components are defined as $f_q^- = \frac{1}{2}(f_q - f_{\bar{q}})$ for $q \in \{1, ..., q_m/2\}$. For $q = 0$, $f_0^+ = f_0$ and $f_0^- = 0$.

The two relaxation parameters are defined as follows: $\lambda^-$ for all anti-symmetric non-equilibrium components $n_q^- = f_q^- - e_q^-$ and $\lambda^+$ for all symmetric non-equilibrium components $n_q^+ = f_q^+ - e_q^+$. Equilibrium components for the zero velocity are given by: $e_0^+ = e_0$ and $e_0^- = 0$.

Equation (A1) is then expressed as [22]:

$$ \begin{aligned} \tilde{f}_q(\vec{r}, t) &= f_q(\vec{r}, t) + \lambda^+ n_q^+ + \lambda^- n_q^- \\ \tilde{f}_{\bar{q}}(\vec{r}, t) &= f_{\bar{q}}(\vec{r}, t) + \lambda^+ n_q^+ - \lambda^- n_q^- \\ \tilde{f}_0(\vec{r}, t) &= f_0(\vec{r}, t)(1 + \lambda^+) - \lambda^+ e_0 \end{aligned} \quad (A3) $$

with $q \in \{1, ..., q_m/2\}$. The LBM-TRT requires the definition of a set of numerical parameters $\Lambda^\pm$ and $\Lambda$ to ensure the stability of the algorithm and to reduce the numerical error. These parameters are related to the relaxation parameters $\lambda^\pm$ as follows [22]:

$$ \Lambda = \Lambda^+ \Lambda^- \quad (A4) $$

with

$$ \Lambda^{\pm} = -(1/2 + 1/\lambda^{\pm}) \quad (A5) $$

for $-2 < \lambda^\pm < 0$. Values of the latter parameters depend on the equation to be solved. Specificities of the resolution of Darcy-Brinkman equation are presented in the following.

### Darcy-Brinkman

Darcy-Brinkman equation, in the absence of gravity, is given by

$$ -\mu_{\text{eff}} \Delta \vec{u} + \frac{\mu}{k} \vec{u} = -\nabla P \quad (A6) $$

or by using $\vec{u} = \vec{j}/\rho_0$ and $v = \mu/\rho_0$ as

$$ -v_{\text{eff}} \Delta \vec{j} + \frac{v}{k} \vec{j} = -\nabla P \quad (A7) $$

where $v_{\text{eff}}$ is the effective Brinkman viscosity, $v$ the kinematic viscosity of the fluid and $k$ the permeability. $\rho_0$ is an average density fixed to $\rho_0 = 1$. Equation (A7) can be expressed as:

$$ \nabla P = \vec{F} + v_{\text{eff}} \Delta \vec{j} \quad (A8) $$

corresponding to the stationary Stokes equation where the force term $\vec{F}$ has the following form:

$$ \vec{F} = -\frac{v}{k} \vec{j} \quad (A9) $$

As can be seen $\vec{F}$ depends on the macroscopic momentum $\vec{j}$. Before each propagation step, $j$ is computed from [23]

$$ \vec{j} = \frac{\vec{J}}{1 + \frac{v}{2k}} \quad (A10) $$

<figure><img src="image_13.png" /><figcaption>Figure A.1. Flow field in the model porous medium for $k_{\text{micro}} = 0.1$.</figcaption></figure>

with $\vec{J} = \sum_{q=1}^{q_m} f_q \vec{v}_q$.

Equation (A10) is then used to compute the equilibrium functions defined as [21, 30]:

$$ \begin{aligned} e_q^+(\vec{r}, t) &= c_s^2 \rho(\vec{r}, t) t_q^* \\ e_q^-(\vec{r}, t) &= t_q^*(\vec{j}(\vec{r}, t) \cdot \vec{v}_q) \\ e_0^+(\vec{r}, t) &= \rho(\vec{r}, t) - \sum_{q=1}^{q_m} e_q^+(\vec{r}, t) \\ e_0^-(\vec{r}, t) &= \rho(\vec{r}, t) - \sum_{q=1}^{q_m} e_q^-(\vec{r}, t) \end{aligned} \quad (A11) $$

where $c_s$ is the sound velocity ($c_s^2 = 1/3$) and $t_q^*$ correspond to isotropic physical weights. $t_q^*$ take the value $t_q^* = \left\{ \frac{1}{3}; \frac{1}{12} \right\}$ for the first (coordinate) and the second (diagonal) neighbor link in the $d_2 q_9$ scheme [21].

The calculation of the equilibrium functions requires the computation of the fluid density $\rho$ which is expressed as [21]: $\rho(\vec{r}, t) = \sum_{q=0}^{q_m} f_q(\vec{r}, t)$. The pressure is defined as $P = c_s^2 \rho$ and the kinematic viscosity as $\nu = \frac{\Lambda^+}{3} = -\frac{1}{3} \left( \frac{1}{2} + \frac{1}{\lambda^+} \right)$. In our case, numerical accuracy is obtained by using $\Lambda = 3/16$ [20] and prescribing $\nu = 1/12$. $\lambda^+$ and $\lambda^-$ can then be computed from equations (A4) and (A5). At the solid/fluide interface (when $\vec{v}_q$ is pointing towards the wall) the bounce back rule is applied: $f_{\bar{q}}(\vec{r}, t + \Delta t) = \vec{f}_q(\vec{r}, t)$.

More details on this method can be found in [21, 23, 30].

### Sensitivity analysis

Figure A.1 shows the flow field in the model porous medium for $k_{\text{micro}} = 0.1$. As can be seen the overall behavior is identical to the one of $k_{\text{micro}} = 1.0$.