Issue 
Sci. Tech. Energ. Transition
Volume 79, 2024
Characterization and Modeling of the Subsurface in the Context of Ecological Transition



Article Number  22  
Number of page(s)  11  
DOI  https://doi.org/10.2516/stet/2024010  
Published online  26 March 2024 
Regular Article
Towards a better comprehension of reactive transport coupling experimental and numerical approaches
^{1}
IFP Energies Nouvelles, 1 et 4 Av. Bois Préau, 92852 Rueil Malmaison, France
^{2}
IFP Energies Nouvelles, Rondpoint de l’échangeur de Solaize, 69360 Solaize, France
^{3}
Laboratoire de Géologie, Ecole Normale Supérieure – PSL Research University, CNRS, UMR8538, Paris, France
^{*} Corresponding author: daniela.bauer@ifpen.fr
Received:
7
September
2023
Accepted:
31
January
2024
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_{2}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 wormholetype channel. To provide further knowledge of the smallscale physics of reactive transport, we performed LatticeBoltzmann 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_{2}saturated water starts to enter highvelocity zones earlier than lowvelocity zones, inducing an earlier onset of the alteration process and a more pronounced local dissolution. However, numerical results showed that the alteration of nonconnected 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 alterationinduced percolating path consisting of macropores, the increase in permeability is relatively low as shown by the LatticeBoltzmann 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.
Key words: Reactive transport / LatticeBoltzmann simulation / Permeability / Homogenization
© The Author(s), published by EDP Sciences, 2024
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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_{2} 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_{2} 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 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_{2} injection in carbonate formations is generally studied through coreflood experiments that can be coupled with scanner imaging to derive 3D dissolution patterns [7–10]. Coreflood 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 DarcyBrinkman 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 DarcyBrinkman 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 threedimensional 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 smallscale 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 smallscale 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 intraparticle 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].
Figure 1 Petrographic overview of Euville limestone (transmitted light planepolarized photomicrographs). Porosity appears in blue. [Left]: Cr: Crinoid ossicle, Ec: Echinoid, La: Lamellibranch bioclast, [Right]: Mi: Micropore, Ma: Macropore. 
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_{2}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.625mmthick slices and data corresponding to pure calcite samples. Microscanner imaging allows a significantly higher resolution (20.28 μm for a 40mm diameter sample) but requires four series of acquisition with an overlapping to image the entire sample. The effective porosity is determined using the threeweight method:(1)where m_{dry} is the mass of the dry sample, m_{sat} is the mass of the sample saturated with a 25 g/l NaCl brine, and m_{imm} is the mass of the sample immersed in the fluid in which it is saturated.
Water permeability is directly measured in the coreflood 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 ΔP. The sample intrinsic permeability K is then deduced from the slope of the Q − ΔP curve using Darcy’s law:(2)where μ 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.
Figure 2 Schematic diagram of coreflood device showing upstream CO_{2}saturated water injection system and downstream effluent sampling system. 
The reactive fluid used for the alteration phase is pure water initially degassed and then kept in equilibrium with CO_{2} 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^{3}/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: Porositycontrolled 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 α [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 holefilling 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 α = 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.
Figure 3 Model porous medium generated with PlugIm. Macropores: black, solid grains: white and microporous matrix: grey. 
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 DarcyBrinkman equation [14, 15] given by(3)where μ_{eff} is the effective Brinkman viscosity, μ the dynamic viscosity of the fluid, and k the local permeability assumed isotropic. is the driving pressure gradient and 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.
From the velocity field the global permeability K_{i}, with i = x, y, of the 2D domain can be computed as:(4)where u_{i,D} is the Darcy velocity, L_{i} is the domain length and ΔP_{i} is the pressure difference in direction i.
2.2.3 LatticeBoltzmann method for the resolution of the DarcyBrinkman equation
In this work, flow is simulated by means of a LatticeBoltzmann TwoRelaxationTime 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 DarcyBrinkman 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 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 steadystate, laminar, singlephase, and incompressible flow of a Newtonian viscous fluid in a porous medium obeys to Darcy’s law at the macroscopic scale:(5)where is the Darcy velocity, μ is the fluid dynamic viscosity, the pressure gradient, and is the secondorder tensor of permeability that depends only on the porous medium geometry. The objective of the homogenization process is to estimate . 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 selfconsistent 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 selfconsistent scheme is particularly wellsuited 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 onethird of the sample length before relocating near the periphery. The preferential dissolution path is directly visible on the sample surface (Fig. 5).
Figure 4 [Top]: Alterationinduced porosity variations deduced from CTscanner image difference, longitudinal sections at different rotation angles (0°, 45°, 90°, 135°, 180°), [Bottom]: CTscanner porosity profile before (solid line) and after (dashed line) alteration. 
Figure 5 Pictures of the sample after alteration (40mm diameter, 100mm length). A: Injection face, B: Outlet face, C: Longitudinal view. 
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 wormholetype channel.
Figure 6 Comparison of microscanner images before and after alteration. [Left]: Position of the microscannerimaged crosssections, [Middle]: Microscanner images before alteration, [Right]: Microscanner images after alteration. 
3.2 Flow field
In this section we present twodimensional simulation results of flow in the model porous medium shown in Figure 3. To this goal we solve the DarcyBrinkman 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_{micro} = 0.1), that the influence of the value of k_{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_{micro}. As the objective of the simulations is to provide further comprehension of the structuredependent flow field on reactive transport and not to match experimental data at this stage, the permeability of the matrix is arbitrarily set to k_{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 highvelocity 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 highvelocity 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.
Figure 7 Velocity field, using limited scale (LatticeBoltzmann units) in order to enhance the contrast between low and highvelocity zones, k_{micro} = 1 (no flow on top and bottom faces, application of a pressure gradient between left and right faces). 
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 crosssection of the sample. In order to better understand the influence of the alteration on the flow field, we selected a subvolume along the highest velocity path observed in the bottom part of the model (SV_{7}, 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 subvolume 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.
Figure 8 [Left]: Initial image of selected subvolume (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 subvolume SV_{7}, [Bottom Right]: Flow fields in the alterated images. 
[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 subvolume SV_{7}. K_{0} corresponds to the effective permeability of subvolume SV_{7} before alteration.
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 subvolumes (see Fig. 9).
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 subvolumes. Each subvolume is characterized by a specific fraction of macropores and microporous matrix. 
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 subvolumes showing rather homogeneous phase distribution but different fractions of macropores and microporous matrix (Tab. 2). Note that subvolume SV_{7} in Figure 9 corresponds to the volume altered in Section 3.2.
Fractions of macropores, microporous matrix, and solid in the different subvolumes.
LB simulations in x and ydirection can be performed on the images of the subvolumes, thus providing 2D directional permeability estimates. The homogenization method is applied to each subvolume 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 subvolume 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 subvolumes is relatively high it is still below the percolation threshold of both LB and EM. Nevertheless, the fact that certain subvolumes are close to the percolation threshold explains the increase in permeability and the difference observed between LB simulations and homogenization.
Figure 10 Effective permeability calculated for the different subvolumes 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 ydirection. 
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 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 subvolume 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 matrixcontrolled fluid flow for macropore fractions below the percolation threshold and a fluid flow dominated by macropores above this threshold. Permeability values of the numerically altered subvolume SV_{7} clearly remain in the regime of the matrixcontrolled fluid flow. If we now 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 85fold increase in permeability is only compatible with macropore fractions approaching the percolation threshold for consecutive subvolumes of the sample. And this is indeed what can be seen on the microscanner images after alteration in Figure 6.
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 subvolumes 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_{total}). 
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_{2}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 numericallygenerated 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 highvelocity zones CO_{2}saturated water input starts earlier than in lowvelocity 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 LatticeBoltzmann simulations. However, once macropores approach percolation, the increase in permeability reaches several orders of magnitude. In conclusion, by applying different experimental, LatticeBoltzmann, and homogenization approaches we were able to provide further understanding of the process of dissolution of porous media due to CO_{2}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 Isabelle Faille for implementing the DarcyBrinkmann equation in the opensource LatticeBoltzmann 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 microscanner analyses. Maxime Moreaud created the 2D model porous structure. Daniela Bauer performed the LatticeBoltzmann simulations and Mathilde Adelinet the homogenization modeling.
References
 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. [Google Scholar]
 Lombard J.M., Azaroual M., Pironon J., Broseta D., Egermann P., Munier G., Mouronval G. (2010) CO_{2} injectivity in geological storages: An overview of program and results of the GeoCarboneInjectivity project, Oil Gas Sci. Technol. 65, 4, 533–539. [CrossRef] [Google Scholar]
 Siqueira T.A., Iglesias R.S., Ketzer J.M. (2017) Carbon dioxide injection in carbonate reservoirs – a review of CO_{2}waterrock interaction studies, Greenh. Gases: Sci. Technol. 7, 5, 802–816. [CrossRef] [Google Scholar]
 André L., Audigane P., Azaroual M., Menjoz A. (2007) Numerical modeling of fluid–rock chemical interactions at the supercritical CO_{2}–liquid interface during CO_{2} injection into a carbonate reservoir, the Dogger aquifer (Paris Basin, France), Energy Convers. Manag. 48, 6, 1782–1797. [CrossRef] [Google Scholar]
 Lønøy A. (2006) Making sense of carbonate pore systems, AAPG Bull. 90, 9, 1381–1405. [CrossRef] [Google Scholar]
 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 dualporosity porenetwork approach: Influence of percolation on the electrical transport properties, Phys. Rev. E 84, 1, 011133. [CrossRef] [PubMed] [Google Scholar]
 Bazin B. (2001) From matrix acidizing to acid fracturing: a laboratory evaluation of acid/rock interactions, SPE Prod. Facil. 16, 1, 22–29. [Google Scholar]
 Golfier F., Zarcone C., Bazin B., Lenormand R., Lasseux D., Quintard M. (2002) On the ability of a darcyscale model to capture wormhole formation during the dissolution of a porous medium, J. Fluid Mech. 457, 213–254. [CrossRef] [Google Scholar]
 Vialle S., Contraires S., Zinzsner B., Clavaud J.B., Mahiouz K., Zuddas P., Zamora M. (2014) Percolation of CO_{2}rich fluids in a limestone sample: Evolution of hydraulic, electrical, chemical, and structural properties, J. Geophys. Res. Solid Earth 119, 4, 2828–2847. [CrossRef] [Google Scholar]
 Menke H.P., Bijeljic B., Andrew M.G., Blunt M.J. (2015) Dynamic threedimensional porescale imaging of reaction in a carbonate at reservoir conditions, Environ. Sci. Technol. 49, 7, 4407–4414. [Google Scholar]
 Vialle S., Dvorkin J., Mavko G. (2013) Implications of pore microgeometry heterogeneity for the movement and chemical reactivity of CO_{2} in carbonates, Geophysics 78, 5, L69–L86. [CrossRef] [Google Scholar]
 Yang Y., Li Y., Yao J., Iglauer S., Luquot L., Zhang K., Sun H., Zhang L., Song W., Wang Z. (2020) Dynamic porescale dissolution by CO_{2}saturated brine in carbonates: Impact of homogeneous versus fractured versus vuggy pore structure, Water Resour. Res. 56, 4, e2019WR026112. [CrossRef] [Google Scholar]
 Peter A., Yang D., Eshiet K.I.I.I., Sheng Y. (2022) A review of the studies on CO_{2}–brine–rock interaction in geological storage process, Geosciences 12, 4, 168. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Brinkman H.C. (1949) On the permeability of media consisting of closely packed porous particles, Flow Turbul. Combust. 1, 1, 81. [CrossRef] [Google Scholar]
 Bemer E., Nguyen M.T., Dautriat J., Adelinet M., Fleury M., Youssef S. (2016) Impact of chemical alteration on the poromechanical properties of carbonate rocks, Geophys. Prospect. 64, 4Advances in Rock Physics, 810–827. [CrossRef] [Google Scholar]
 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. [CrossRef] [MathSciNet] [Google Scholar]
 Serra J. (1988) Image analysis and mathematical morphology, part II: Theoretical advances. Academic Press [Google Scholar]
 “plugim!” (2018) An open access and customizable software for signal and image processing. https://www.plugim.fr. [Google Scholar]
 Ginzburg I. (2008) Consistent Lattice Boltzmann schemes for the brinkman model of porous flow and infinite ChapmanEnskog expansion, Phys. Rev. E 77, (6), 066704. [CrossRef] [PubMed] [Google Scholar]
 Ginzburg I. (2007) Lattice Boltzmann modeling with discontinuous collision components: Hydrodynamic and advectiondiffusion equations, J. Stat. Phys. 126, 157–206. [CrossRef] [MathSciNet] [Google Scholar]
 Ginzburg I., d’Humières D., Kuzmin A. (2010) Optimal stability of advectiondiffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium, J. Stat. Phys. 139, 6, 1090–1143. [CrossRef] [MathSciNet] [Google Scholar]
 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. [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 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. [CrossRef] [MathSciNet] [Google Scholar]
 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. [CrossRef] [MathSciNet] [Google Scholar]
 Fokker P.A. (2001) General anisotropic effective medium theory for the effective permeability of heterogeneous reservoirs, Transp. Porous Media 44, 2, 205–218. [CrossRef] [MathSciNet] [Google Scholar]
 Dormieux L., Kondo D. (2004) Approche micromécanique du couplage perméabilité–endommagement, CR Mécanique 332, 2, 135–140. [CrossRef] [Google Scholar]
 Stauffer D., Aharony A. (1992) Introduction to percolation theory, Taylor and Francis, London. [Google Scholar]
 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. [Google Scholar]
 Talon L., Bauer D., Gland N., Youssef S., Auradou H., Ginzburg I. (2012) Assessment of the two relaxation time LatticeBoltzmann scheme to simulate Stokes flow in porous media, Water Resour. Res. 48, 4. [CrossRef] [Google Scholar]
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 , q ∈ [0,q_{m}], q_{m} = 8. Particles are only displaced with velocities . Populations are defined as the density of particles at position with velocity . 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 Ω, specific for the equation to be solved. Distribution functions f_{q} are updated as follows:(A1)where stands for the postcollision distributions.
Propagation step: Here, distribution functions f_{q} are exchanged between neighboring nodes with respect to the velocity set:(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 while the antisymmetric components are defined as for . For q = 0, and .
The two relaxation parameters are defined as follows: λ ^{−} for all antisymmetric nonequilibrium components and λ ^{+} for all symmetric nonequilibrium components . Equilibrium components for the zero velocity are given by: and .
Equation (A1) is then expressed as [22]:(A3)with . The LBMTRT requires the definition of a set of numerical parameters Λ^{±} and Λ to ensure the stability of the algorithm and to reduce the numerical error. These parameters are related to the relaxation parameters λ ^{±} as follows [22]:
(A4)with(A5)for . Values of the latter parameters depend on the equation to be solved. Specificities of the resolution of DarcyBrinkman equation are presented in the following.
DarcyBrinkman
DarcyBrinkman equation, in the absence of gravity, is given by(A6)or by using and ν = μ/ρ_{0} as(A7)where ν_{eff} is the effective Brinkman viscosity, ν the kinematic viscosity of the fluid and k the permeability. ρ_{0} is an average density fixed to ρ_{0} = 1. Equation (A7) can be expressed as:(A8)corresponding to the stationary Stokes equation where the force term has the following form:(A9)
As can be seen depends on the macroscopic momentum . Before each propagation step, j is computed from [23](A10)with .
Equation (A10) is then used to compute the equilibrium functions defined as [21, 30]:(A11)where c_{s} is the sound velocity () and correspond to isotropic physical weights. take the value 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 ρ which is expressed as [21]: . The pressure is defined as and the kinematic viscosity as . In our case, numerical accuracy is obtained by using Λ = 3/16 [20] and prescribing ν = 1/12. λ ^{+} and λ ^{−} can then be computed from equations (A4) and (A5). At the solid/fluide interface (when is pointing towards the wall) the bounce back rule is applied: .
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_{micro} = 0.1. As can be seen the overall behavior is identical to the one of k_{micro} = 1.0.
Figure A.1 Flow field in the model porous medium for k_{micro} = 0.1. 
All Tables
[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 subvolume SV_{7}. K_{0} corresponds to the effective permeability of subvolume SV_{7} before alteration.
Fractions of macropores, microporous matrix, and solid in the different subvolumes.
All Figures
Figure 1 Petrographic overview of Euville limestone (transmitted light planepolarized photomicrographs). Porosity appears in blue. [Left]: Cr: Crinoid ossicle, Ec: Echinoid, La: Lamellibranch bioclast, [Right]: Mi: Micropore, Ma: Macropore. 

In the text 
Figure 2 Schematic diagram of coreflood device showing upstream CO_{2}saturated water injection system and downstream effluent sampling system. 

In the text 
Figure 3 Model porous medium generated with PlugIm. Macropores: black, solid grains: white and microporous matrix: grey. 

In the text 
Figure 4 [Top]: Alterationinduced porosity variations deduced from CTscanner image difference, longitudinal sections at different rotation angles (0°, 45°, 90°, 135°, 180°), [Bottom]: CTscanner porosity profile before (solid line) and after (dashed line) alteration. 

In the text 
Figure 5 Pictures of the sample after alteration (40mm diameter, 100mm length). A: Injection face, B: Outlet face, C: Longitudinal view. 

In the text 
Figure 6 Comparison of microscanner images before and after alteration. [Left]: Position of the microscannerimaged crosssections, [Middle]: Microscanner images before alteration, [Right]: Microscanner images after alteration. 

In the text 
Figure 7 Velocity field, using limited scale (LatticeBoltzmann units) in order to enhance the contrast between low and highvelocity zones, k_{micro} = 1 (no flow on top and bottom faces, application of a pressure gradient between left and right faces). 

In the text 
Figure 8 [Left]: Initial image of selected subvolume (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 subvolume SV_{7}, [Bottom Right]: Flow fields in the alterated images. 

In the text 
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 subvolumes. Each subvolume is characterized by a specific fraction of macropores and microporous matrix. 

In the text 
Figure 10 Effective permeability calculated for the different subvolumes 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 ydirection. 

In the text 
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 subvolumes 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_{total}). 

In the text 
Figure A.1 Flow field in the model porous medium for k_{micro} = 0.1. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.