© The Author(s), published by EDP Sciences, 2024

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 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 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 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 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 plane-polarized 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₂-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 entire sample. The effective porosity is determined using the three-weight 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 CO2-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₂ 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 α [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 α = 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 Darcy-Brinkman 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 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 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:(5)where is the Darcy velocity, μ is the fluid dynamic viscosity, the pressure gradient, and is the second-order 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 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).

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.

Figure 5 Pictures of the sample after alteration (40-mm diameter, 100-mm 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 wormhole-type channel.

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.

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

Figure 7 Velocity field, using limited scale (Lattice-Boltzmann units) in order to enhance the contrast between low and high-velocity zones, kmicro = 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 cross-section of the sample. In order to better understand the influence of the 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.

Figure 8 [Left]: Initial image of selected sub-volume (SV7 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 SV7, [Bottom Right]: Flow fields in the alterated images.

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

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.

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

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.

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 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 the percolation threshold and a fluid flow dominated by macropores above this threshold. Permeability values of the numerically altered sub-volume SV₇ clearly remain in the regime of the matrix-controlled 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 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.

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

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

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 post-collision 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 anti-symmetric components. Symmetric components are defined as while the anti-symmetric components are defined as for . For q = 0, and .

The two relaxation parameters are defined as follows: λ ⁻ for all anti-symmetric non-equilibrium components and λ ⁺ for all symmetric non-equilibrium components . Equilibrium components for the zero velocity are given by: and .

Equation (A1) is then expressed as [22]:(A3)with . The LBM-TRT 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 Darcy-Brinkman equation are presented in the following.

Darcy-Brinkman

Darcy-Brinkman equation, in the absence of gravity, is given by(A6)or by using and ν = μ/ρ₀ as(A7)where ν_eff is the effective Brinkman viscosity, ν the kinematic viscosity of the fluid and k the permeability. ρ₀ is an average density fixed to ρ₀ = 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₂q₉ 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 kmicro = 0.1.

All Tables

All Figures

	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.
In the text

	Figure 2 Schematic diagram of coreflood device showing upstream CO2-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]: 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.
In the text

	Figure 5 Pictures of the sample after alteration (40-mm diameter, 100-mm 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 microscanner-imaged cross-sections, [Middle]: Microscanner images before alteration, [Right]: Microscanner images after alteration.
In the text

	Figure 7 Velocity field, using limited scale (Lattice-Boltzmann units) in order to enhance the contrast between low and high-velocity zones, kmicro = 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 sub-volume (SV7 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 SV7, [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 sub-volumes. Each sub-volume is characterized by a specific fraction of macropores and microporous matrix.
In the text

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

In the text

	Figure A.1 Flow field in the model porous medium for kmicro = 0.1.
In the text