The governing equations consider flow, transport and reaction at the fluid-solid interface. The domain Ω is partitioned into fluid Ω _f and solid Ω _s . Under isothermal conditions and in the absence of gravitational effects, fluid motion in Ω _f is governed by the incompressible Navier-Stokes equations ∇ ⋅ u = 0 , (1) ∂ u ∂ t + ∇ ⋅ u ⊗ u = − ∇ p + ν ∇ 2 u , (2) with the continuity condition at the fluid-solid interface Γ, ρ u − w s ⋅ n s = − ρ s w s ⋅ n s a t Γ , (3) where u (m/s) is the velocity, p (m²/s²) is the kinematic pressure, ν (m²/s) is the kinematic viscosity, ρ (kg/m³) is the fluid density, ρ _s (kg/m³) is the solid density, n _s is the normal vector to the fluid-solid interface pointing toward the solid phase, and w _s (m/s) is the velocity of the fluid-solid interface, which is controlled by the surface reaction rate R (kmol/m²/s) such that w s = M w s ρ s R n s , (4) where M _ws is the molecular weight of the solid. The concentration c (kmol/m³) of a species in the system satisfies an advection-diffusion equation ∂ c ∂ t + ∇ ⋅ c u = ∇ ⋅ D ∇ c , (5) where D (m²/s) is the diffusion coefficient. The chemical reaction occurs at the fluid-solid interface Γ, such that c u − w s − D ∇ c ⋅ n s = ζ R a t Γ , (6) where ζ is the stoichiometric coefficient of the species in the reaction. In this work, we assume that the surface reaction rate depends only on the concentration of one reactant species, following R = k c c , (7) where k _c (m/s) is the reaction constant.

In the micro-continuum approach, the entire domain Ω is considered, i. e fluid Ω _f and solid Ω _s , and the fluid-solid interface is tracked in terms of V _f and V _s , the volume of fluid and solid phase in each control volume V, and their volume fraction ɛ = V _f /V and ɛ _s = 1 − ɛ. The flow, transport and chemical reaction are solved in term of the volume-averaged velocity u ̄ = 1 V ∫ V f u d V , (8) and the phase-averaged pressure and reactant concentration. p ̄ = 1 V f ∫ V f p d V , (9) c ̄ = 1 V f ∫ V f c d V . (10) The averaging process results in an extension of the Darcy-Brinkman-Stokes equation (Soulaine et al., 2017) 1 ε ∂ u ̄ ∂ t + ∇ ⋅ u ̄ ⊗ u ̄ ε = − ∇ p ̄ + ν ε ∇ 2 u ̄ − ν k − 1 u ̄ , (11) where k (m²) is the permeability of the cell. ν k − 1 u ̄ represents the momentum exchange between the fluid and the solid phase, i.e., the Darcy resistance. This term is dominant in the solid phase and vanishes in the fluid phase. To model this, the local permeability field k is assumed to be a function of the local porosity ɛ, following a Kozeny-Carman relationship k = k 0 ε 3 1 − ε 2 , (12) where k ₀ (m²) is the Kozeny-Carman constant. For the acid transport, the mass-balance equation averaged over the control volume gives ∂ ε c ̄ ∂ t + ∇ ⋅ c ̄ u ̄ − ∇ ⋅ ε D ∗ ∇ c ̄ + R ̄ = 0 , (13) where ɛD* (m²/s) is the effective diffusion coefficient and R ̄ (kmol/m³/s) is the volume-averaged surface reaction rate. The effective diffusion coefficient takes into account the reduction of the total diffusion due to the presence of solid phase. In this paper, we take D∗ = D. The volume-averaged surface reaction rate is defined as R ̄ = 1 V ∫ A k c c d S , (14) where A = V ∩Γ is the reactive surface area in the control volume. The specific surface area a _s (m⁻¹) in a control volume is defined as a s = 1 V ∫ A d S . (15) During the dissolution process, the pore-space evolves with the chemical reaction at the pore surfaces. The number of moles of solid consummed by the chemical reaction is equal to the number of moles of reactant consummed. Therefore, the mass balance equation for the solid phase writes ∂ ρ s ε ∂ t = R ̄ M w s . (16) and the volume averaged velocity satisfies ∇ ⋅ u ̄ = R ̄ M w s 1 ρ s − 1 ρ . (17)

The method presented here is analogue to the calculation of the mass transfer across a multiphase interface presented in Maes and Soulaine (2020), for which the mass transfer is calculated as the scalar product between a diffusive flux and the gradient of the phase indicator function. To calculate the volume-averaged surface reaction rate, the improved Volume-of-Solid (iVoS) therefore introduces the reactive flux, Φ _R (kmol/m²/s), defined as Φ R = k c c n s , (18) and the volume-averaged surface reaction rate can be rewritten as R ̄ = 1 V ∫ A Φ R ⋅ n s d S . (19) Assuming that the concentration of the reactant on the reactive surface can be approximated by its volume-averaged on the control volume, and that the normal vector to the interface can be approximated by n ̄ s = − ∇ ε ‖ ∇ ε ‖ , (20) the reactive flux can be approximated by Φ ̄ R = k c c ̄ n ̄ s . (21) Moreover, the average surface normal in a control volume can be calculated as (Quintard and Whitaker, 1994) 1 V ∫ A n s d S = − ∇ ε . (22) Therefore, the volume-averaged surface reaction rate can be calculated as R ̄ = − Φ ̄ R ⋅ ∇ ε . (23) To avoid problems related to the calculation of the gradient of ɛ (see Section 2.4), the divergence theorem is used to recast R ̄ as R ̄ = ε ∇ ⋅ Φ ̄ R − ∇ ⋅ ε Φ ̄ R , (24) With this formulation, the reactive rate is the sum of two terms, an overall mass transfer term ( ∑ R ̄ f = ∑ ε ∇ ⋅ Φ ̄ R ) and a conservative term ( − ∑ ∇ ⋅ ε Φ ̄ R = 0 ) , which balances the local reaction rate between two adjacent control volumes. This means that the reaction can consume reactant in 1 cell and use it to dissolve solid in a neighbor cell. By using a second-order high resolution difference scheme (van Leer, 1974), the reaction rate is balanced toward the reactive surface and the diffusion of the solid interface is limited. This is an accurate representation of a reaction at an interface between 2 cells with ɛ ≈ 0 and ɛ = 1, where the reaction rate is calculated using the concentration in the fluid where ɛ = 1 and the reaction dissolves the solid where ɛ ≈ 0. This formulation is labelled iVoS.

As an alternative to Eq. 24, the volume-averaged surface reaction rate can be calculated as R ̄ = k c c ̄ a s . (25) The specific surface area can be directly calculated as = ‖∇ɛ‖. However, this can lead to a diffuse interface that spreads across a large number of layers in the computational grid. To enforce localization of the dissolution front on the fluid-solid interface, the VoS-psi method introduces a diffuse interface localization function ψ (Soulaine et al., 2017) so that R ̄ = k c c ̄ ψ ‖ ∇ ε ‖ . (26) The main advantage of VoS-ψ compared to iVoS is that it provides a direct calculation of the reactive surface area in a control volume. Therefore, the reaction rate can be calculated by a dedicated geochemical solver, such as Phreeqc (Parkhurst and Wissmeier, 2015) or Reaktoro (Walsh et al.), as it is often done for standard Reactive Transport Modelling dedicated to multi-scale applications (Soulaine et al., 2021b).

While the VoS-ψ method has proved to be a fast and flexible method to match experimental results (Molins et al., 2020), this formulation has one main limitation. It is strongly dependent on the choice of the localization function ψ. Several functions have been proposed by Luo et al. (2012), but their accuracy depends on the case considered and on the discretization scheme used for the gradient. For example, using a centered difference scheme requires ψ(1.0) = 0 for stability. However, using a decentered scheme (in the direction of ɛ = 0 to avoid instabilities) will result in a higher and more diffuse reaction rate, due to a higher reactant concentration away from the surface and a larger specific surface area when ɛ = 0. Centered difference schemes are less diffuse, but they result in incomplete dissolution, since a cell with ɛ < 1 but ɛ = 1 for all its neighbors will have a zero reaction rate. Currently, there is no consensus on the ideal combination of localization function and discretization scheme to use for every scenario, as this will depend on the geometry and flow conditions.

In this work, we use a centered difference scheme for the gradient and ψ = λɛ(1 − ɛ), which is the most accurate combination of discretization scheme and localization function proposed by Luo et al. (2012) for the case of dissolution of a calcite post by acid injection (Soulaine et al., 2017). Typically, λ = 4, but since ψ ≤ 1.0, this will lead to a reduction in interfacial area. For this reason, the VoS-ψ method gives an overall lower dissolution rate that the iVoS method. Instead, λ can be calculated as λ ε = ∫ Ω ‖ ∇ ε ‖ ∫ Ω ‖ ∇ ε ‖ ε 1 − ε . (27) In Equ. (27), λ _ɛ is a parameter that depends on the porosity field ɛ in the whole domain and chosen so that the surface area is conserved globally, i. e ∫ Ω ‖ ∇ ε ‖ λ ε ε 1 − ε = ∫ Ω ‖ ∇ ε ‖ . (28) In this work, we label VoS-ψ the formulation using ψ = 4 ε 1 − ε and VoS-ψ′ the formulation using ψ = λ ε ε 1 − ε . Using the VoS-ψ′ formulation, the total surface area is conserved globally. However, it is not conserved locally. At an interface between 2 cells with ɛ ≈ 0 and ɛ = 1, the interface area will remain close to 0 in both cells even after the correction.

In order to investigate the capabilities of our numerical model to calculate upscaled properties for macro-scale simulations, the flow and reaction in the whole domain are characterised by the total porosity ϕ, the Darcy velocity U _D (m/s), U D = Q i A i , (29) and the permeability K (m²), K = ν U D L D Δ P . (30) where Q _i (m³/s) is the inlet flow rate, A _i (m²) is the inlet area, L _D (m) is the distance between the inlet and outlet, ΔP (m²/s²) is the kinematic pressure drop. In addition, the chemical reaction is characterised by the total specific surface area a (m⁻¹) a = 1 V Ω ∫ Γ d S , (31) and the correction factor α, that represents the reduction of the reactive surface area accessible to reactant, and is defined as α = 1 a ∫ Γ c d S ∫ Ω f c d V . (32)

Change in effective upscaling parameters K, a and α are often modelled as a function of porosity with power law functions (Noiriel et al., 2009; Wen and Li, 2017; Seigneur et al., 2019). The parameters of these power-law functions depend strongly on the flow, transport and reaction conditions, which are characterized by the Péclet number P e = U L D , (33) which quantifies the relative importance of advective and diffusive transport, and the Damköhler number D a = k c U , (34) which quantifies the relative importance of chemical reaction and advective transport. Here U and L are the reference velocity and length. The product of the Damköhler and Péclet numbers is also a relevant quantity called the Kinetic number, defined as K i = D a P e = k c L D . (35) In addition, the reactant strength is characterized by β = c i M w s ζ ρ s , (36) where c _i is the concentration of reactant at the inlet.

In benchmark 1, we dissolve a calcite post with all four numerical methods (ALE, VoS-ψ, VoS-ψ′ and iVoS) and compare the results with the experimental data from Soulaine et al. (2017). In the experiment, an octagonal-shaped calcite post is placed at the center of a straight microchannel and is dissolved by an acidic solution that is flowing past. The experiment and methods are described in detail in Soulaine et al. (2017).

Simulations are performed using the VoS methods, and using the ALE method for reference. The domain is a straight microchannel of size 2.67 mm × 1.5 mm × 0.2 mm with a calcite post of height 0.2 mm in its center. Images of the post (in stl and h5 formats) are given in the supplementary material. In each case, a Cartesian mesh of resolution Δx = 20 μm is generated, which is snapped onto the calcite solid surfaces for the ALE method (C). Increased resolutions of Δx = 10 μm and Δx = 5 μm at the interface are obtained using LMR or AMR (C). The initial meshes have respectively 95,900, 109,476 and 160,390 cells for the ALE simulations and 105,000, 114,430 and 174,280 for both of the micro-continuum simulations.

At t = 0, a solution of hydrochloride acid is injected from the left boundary at constant flow rate, extrapolated from a zero-gradient pressure velocity field (OpenCFD, 2016), and constant concentration c _i = 0.0126 kmol/m³. The acid reacts with the calcite surface to produce Ca ²⁺ and H C O 3 − . However, due to the very low pH, H C O 3 − reacts instantaneously with H ⁺ to give H 2 C O 3 * . The two reactions can be added and modelled as the single reaction described below. C a C O 3 + H + ⇋ C a 2 + + H C O 3 − H C O 3 − + H + ⇋ H 2 C O 3 C a C O 3 + 2 H + ⇋ C a 2 + + H 2 C O 3 ∗ (37) The simulation parameters are summarized in Table 1. Each simulation is run until t = 12,000 s or until all of the solid has been dissolved, whichever happens first.

Figure 1 shows the concentration map at different times for the different methods, with c = 0 in the solid phase. We see that for the ALE and iVoS methods, the concentration maps are very similar and the calcite post has fully disappeared at t = 12,000 s. However, for the VoS-ψ method, the dissolution is delayed and there is still a small but significant volume of calcite at t = 12,000 s. The VoS-ψ’ method corrects most of this error and the concentration maps are similar to the one obtained with the ALE and iVoS methods.

The results of numerical simulations at different mesh resolutions are compared with experimental results from Soulaine et al. (2017) in Figure 2. We observe that the ALE method converges toward a solution close to the experimental results. The difference between the experimental results and the numerical simulation at mesh resolution Δx = 5 μm in grain volume and grain area are less than 1% of the initial values. Although the iVoS method gives a significantly lower grain volume and surface area than the experiment at a resolution Δx = 20 μm, the results are very similar to the ALE results for Δx = 10 μm and Δx = 5 μm. In addition, with the iVoS method as with the ALE method, the difference between the experimental results and the numerical simulation at mesh resolution Δx = 5 μm in grain volume and grain area are less than 1% of the initial values.

Although the VoS-ψ method matches the trend of the experiment, it overestimates the grain volume and surface area for all resolutions. The results do not improve as the mesh resolution increases and converge toward a solution with an error of 5% in the grain volume and 10% in the grain area. This error is corrected by the surface area correction provided by the VoS-ψ′ method. At all resolution, the VoS-ψ’ method gives errors in grain volume and grain area that are less than 1% of the initial values. During the simulation, λ _ɛ is a value between 8 and 16 that changes as ɛ changes. This shows that the VoS-ψ underestimates the overall surface area by a factor between 2 and 4.

Table 2 shows the CPU times for all simulations. The simulations are significantly faster using the VoS methods, with the iVoS being approximately 2.5x faster than the ALE for each simulation. For Δx = 20 μm and Δx = 10 μm, the VoS-ψ method is slightly faster than the iVoS and VoS-ψ′ methods, which is due to the larger time steps (average of 2.1 s for VoS-ψ compared to 1.7 s for iVoS and VoS-ψ′ for the simulation with Δx = 10 μm) and a faster convergence of the steady-state solver (average of six iterations per time-step for VoS-ψ compared to seven iterations for iVoS and VoS-ψ’ for the simulation with Δx = 10 μm). These differences disappear for the simulation with Δx = 5 μm.

We conclude that all VoS methods are significantly faster than the ALE method. The iVoS provides a result that converges toward a solution with an error less than 1% when the mesh resolution reaches 5 μm. The VoS-ψ results in a small error which is due to the reduction of the overall surface area and is corrected in the VoS-ψ′ method by fitting the constant λ for each time-step. The VoS-ψ’ method gives accurate results for all mesh resolution.

In benchmark 2, we use the three VoS methods to model the various dissolution regimes (i.e compact, wormholes, and uniform) that occur during mineral dissolution in a 2D porous media model. ALE simulations are also run for reference. The model is constructed from a homogeneous domain with discs radius 270 μm by adding a random deviation of magnitude 270 μm in disc radius and center position. The micromodel generation code is available open source (https://github.com/hannahmenke/DrawMicromodels) and the method is described in Patsoukis-Dimou et al. (2021). The geometry is presented in Figure 3 and a high resolution image can be found in the Supplementary Material.

The domain is meshed with a Cartesian mesh with uniform resolution Δx = 3 μm, which is snapped on the solid surfaces for the ALE method. A band of 2 cells width is added on each side of the model to avoid dissolution at the boundaries. The final meshes include 1,008,016 cells for the VoS methods and 457,455 for the ALE method. The porosity ϕ and the permeability K of the full domain can be numerically calculated as ϕ = 0.45 and K = 5.6 × 10^–10 m².

At t = 0, acid is injected from the left boundary at constant flow rate, extrapolated from a zero-gradient pressure velocity field (OpenCFD, 2016), and flows out of the domain from the right boundary at constant pressure. The top and bottom boundaries are no-flow, no-slip conditions. The fluid and solid properties are summarized in Table 3. The Kozeny-Carman constant is fitted to obtain the same permeability as in the direct method at t = 0. For each simulation, the inlet flow rate Q _i and the chemical reaction constant k _c are adapted to obtain the correct Pe and Ki, using the pore-scale length L = 12 K / ϕ as the reference length and the average pore velocity U = U _D /ϕ as the reference velocity. The factor 12 is added so that the pore-scale length corresponds to the channel size for an homogeneous bundle of straight channels (Pavuluri et al., 2018). Four cases are considered that characterize four different dissolution regimes: Pe = 0.01 and Ki = 0.1 (compact dissolution), Pe = 0.3 and Ki = 10 (conical wormhole), Pe = 1 and Ki = 1 (dominant wormhole), and Pe = 10 and Ki = 0.01 (uniform dissolution). The simulations are performed until 20% of the solid has been dissolved.

Figure 4 shows the concentration map at the end of the simulations. We observe that all methods are able to model all regimes qualitatively. In particular, the iVoS method reproduces the ALE results almost exactly. However, there are several inaccuracies in the VoS-ψ method. Due to the use of a centered scheme for the computation of the gradient of ɛ, the VoS-ψ method results in several grains with incomplete dissolution in the compact and wormhole regimes. This problem can be resolved by applying a decentered gradient, but this is done at the expense of accuracy, as a decentered gradient generates more numerical diffusion. In addition, the wormhole obtained with VoS-ψ for Pe = 1 and Ki = 1 is more diffused and ramified that the ones obtained with ALE and iVoS, due to the reduction of reaction rate induced by ψ. For the uniform regime, the patterns are almost identical, but the acid concentration obtained with the VoS-ψ method is higher than with ALE and iVoS, suggesting that the reaction rate is lower for VoS-ψ. Most of these errors are corrected in the VoS-ψ’ method, although there are still grains with incomplete dissolution in the compact and wormhole regime, and the wormhole is still slightly more diffuse and ramified for Pe = 1 and Ki = 1.

Figures 5, 6 show the evolution of the solid volume and solid surface area for all regimes obtained with all numerical methods. We observe that the iVoS method reproduces the ALE results with good accuracy for all regimes. For all cases, the solid volume and surface area evolutions obtained with the ALE and the iVoS methods are similar. The VoS-ψ method underpredicts the amount of dissolution occurring in the systems compared to the ALE and the iVoS methods. This is particularly true for the cases with Ki ≤ 0.1, i.e. in the compact and uniform regime. In these cases, the concentration gradient on the reactive surface is small and the reaction rate is mostly dependent on the surface area, which is significantly reduced by the localization function ψ. For the dominant wormhole regime, although the solid volume evolution obtained with the VoS-ψ method is similar to the ones obtained with the ALE and iVoS methods, the evolution of the solid area is significantly different. This suggests that, although the total amount of dissolution is correctly predicted, it occurs in a slightly different pattern than with the ALE and iVoS methods. The dissolution front is less sharp, and the reaction occurs in the vicinity of the wormhole rather than at the tip of the wormhole as predicted by the ALE and iVoS methods. This leads to a more ramified dissolution pattern, as observed in Figure 4. The VoS-ψ’ method corrects some of these errors and the overall amount of dissolution is similar to the ALE and iVoS results for all cases. However, the surface area is slightly higher for all cases. This is because, since the dissolution is high when ɛ ≈ 0.5 and low when ɛ ≈ 0 and ɛ ≈ 1, the interface becomes artificially sharp and leads to a larger interfacial area.

Figure 7 shows the evolution of the permeability K and the correction factor α as a function of porosity for all regimes obtained with the four methods. In all cases, the order of the permeability evolution is similar for the three methods, except for the dominant wormhole regime, for which the order is 19 for ALE and 17 for iVoS, typical of the dominant wormhole regime, but only 8 for VoS-ψ and 12 for VoS-ψ′, which is more typical of the ramified wormhole regime. Compared to the ALE method, the difference of permeability (in log scale) at the end of the simulation is 13% for the iVoS method, 40% for the VoS-ψ′ method and 60% for the VoS-ψ method. The macro-scale reaction constant is consistently lower with VoS-ψ than with ALE and iVoS, especially for the uniform regime where the difference with the ALE result at the end of the simulation is 2% for the iVoS method, 8% for the VoS-ψ’ method and 78% for the VoS-ψ method.

Table 4 shows the CPU time for all cases for all three methods. We observe that the VoS methods are between 3 and 12 time faster than the ALE method. The VoS methods are particularly efficient compared to the ALE method for the cases with localized dissolution front, i.e., for the compact and the wormholing regime, for which the ALE method performs a large number of remeshing steps. The VoS-ψ method is slightly faster than the iVoS method, but this is mostly due to larger time steps resulting from the lower dissolution rates. The VoS-ψ’ method is slightly slower than the iVoS method, which we attribute to a slower convergence of the transport equation due to a more localized dissolution.

We conclude that the iVoS method is capable of modelling all regimes during dissolution in a 2D micromodel and calculating macro-scale coefficients with similar results as the ones obtained with ALE. In addition, the iVoS method is significantly faster than the ALE method. The VoS-ψ method, using ψ = 4ɛ _f (1 − ɛ _f ) underpredicts the dissolution in all cases due to a reduction of overall surface area, but this error is mostly corrected by using the VoS-ψ’ method.

In benchmark 3 we simulate the dissolution of a 3D micro-CT image of Ketton limestone and compare the numerical results to in-situ dissolution experiments. The main advantage of the micro-continuum approach is its computational efficiency compared to the ALE method. In this part, we take advantage of this to perform simulations with a VoS method, while the ALE method is too computationally expensive to perform on a 3D image of this size. The iVoS method was selected since it gave the most accurate results in benchmark 2. The experiment is the one conducted in Menke et al. (2015), where CO₂-saturated brine is injected in a Ketton carbonate core. A core of 4 mm diameter was flooded with a brine solution representative of a typical saline aquifer consisting of 1% KCl and 5% NaCl by weight, pre-equilibrated with supercritical CO₂ at 10 MPa and 50 C. Micro-CT images were acquired after 17, 33, 50 and 67 min.

The simulated conditions are similar to the ones presented in Pereira Nunes et al. (2016). The sample is a cylinder of radius 1.7 mm and length 3.5 mm. The image has 911 × 902 × 922 voxels with resolution 3.8 μm, however, we ran it on two-level adpative mesh with maximum resolution of 7.6 μm to decrease computation time. The initial grid includes 17 million cells. Figure 8A shows the initial porosity field. The velocity field is then initialised by solving Equ. 17) and 11) in the domain using k ₀ = 2 × 10^–13 m² which was fitted to obtain a permeability of 16D. The streamlines and corresponding magnitude of velocity are shown on Figure 8B. The reaction rate is described in terms of the concentration of calcium cations Ca²⁺, following R = k eff C e q − C , where C _eq and k _eff depend on the concentrations of C O 3 2 − , H ⁺ and C O ₂. Because our model only tracks one component, the reaction constant and equilibrium concentration are fitted to match the reaction-limited constant obtained for a flat pure crystal of calcite (Peng et al., 2015) R = 8.1 × 10^–7 kmol/m²/s and the initial reaction rate in the experiment R = 8.8 × 10^–8 kmol/m²/s. Such rates are obtained using a reactant concentration of 0.0035 kmol/m³ at the inlet and a reaction constant k _c = 2.314 × 10^–4 m/s. Figure 8C shows the concentration field in the pores at T = 0 s obtained by solving Equ. (A.4).

The simulation parameters are summarized in Table 5. We use for reference velocity the average pore velocity defined as U = U D ϕ . (38) The reference length is associated to the initial permeability so that L = 8 K ϕ , (39) where the constant 8 is added in order to obtain the throat radius for a capillary bundle of uniform size. With these definitions, we obtain U = 5.3 × 10^–3 m/s and L = 27 μm. Using the diffusion coefficient and reaction constant in Table 5, we obtain Pe = 190 and Da = 4.4 × 10^–2.

The simulation is run with an adaptive time-step until t = 4,000 s on 128 CPUs using Oracle cloud computing. The total CPU time was 65 h. Figure 9 shows a comparison of the evolution of the dissolution pattern between experiment and simulation. The colors show the part of the rock that is dissolved for each time interval. We observe that the simulation captures the correct patterns, with pores being enlarged in the direction of the flow. Some expected differences are observed, primarily due to uncertainties in the experimental conditions, segmentation error due to reaction occurring during image acquisition, and the fact that the simulation is performed on a sub-image of the full sample used in the experiment. Furthermore, some of the differences in the local porosity at different times in the experiment are due to imperfect alignment of the images. Table 6 shows the evolution of the porosity during dissolution for both experiment and simulation. We observe that the porosity is accurately predicted by the simulation until T > 3,000 s, at which time the simulation starts diverging from the experiment slightly. This can be explained by the fact that the simulation is done on a sub-image of the full sample, starting at 2 mm away from the inlet (Menke et al., 2015), and thus the concentration of acid will have varied as dissolution occured in the unimaged portion of the core.

In addition, permeability can be calculated during the simulation. Table 6 shows the porosity and permeability at different times for the simulation and the experiment and we observe a good correspondence. The permeability is then plotted as a function of the porosity in Figure 10. We see that the permeability in the simulation follows the same trend that of the experiment. Finally, the average macro-scale reaction rate for each time interval can be calculated as R Ω = ρ s Δ ϕ M w Δ t , (40) for both simulation and experiment. The correction factor α can then be calculated by dividing by the average concentration in the pore space, calculated at the end of the time interval from the simulation results. The values are summarized in Table 7. We observe a good correspondence between simulation and experiment. The overall trend of decreasing reaction rates and correction factors observed in the experiment is reproduced in the simulation, albeit with a slightly slower rate. This decrease is not a typical characteristic of the uniform regime, and is an indication that the dissolution might be occuring in the channeling regime, identified in Menke et al. (2017).

Two additional simulations are run in different regimes by dividing the reaction rate by 100 and 10,000, which gives Péclet number of 1.9 and 0.019. The simulations are run on 128 CPUs using Oracle cloud computing until the porosity reaches an approximate value of 0.28. The CPU time was 75 h for Pe = 1.9 and 80 h for Pe = 0.019. Figure 11 shows the reactant concentration along reconstructed streamlines at the end of the simulation. At Pe = 0.019, the dissolution is in the compact regime. The reactant is consumed and dissolves the rock close to the inlet face. At Pe = 1.9, the dissolution is in the wormholing regime. The reactant penetrates further in the domain and flow instabilities kick in, leading to a preferential dissolution pathway. The streamlines and reactant concentration for Pe = 190 are also shown in Figure 11 and the dissolution appears to be in the uniform regime.

Figure 12A shows the evolution of permeability as a function of porosity during the simulations for the three regimes. The permeability of the compact and uniform can be fitted with power laws, and we obtain K = 30ϕ ^0.36 for the compact regime and K = 5.5 × 10⁴ ϕ ^4.7 for the uniform regime. For the wormhole regime, the permeability curve has a strong inflection point around ϕ = 0.25, which corresponds to porosity at reactant breakthrough time. This is typical of the wormholing regime. Before breakthrough occurs, the wormhole does not provide a connected pathway throughout the sample, so the permeability increases with low order. When acid reaches the outlet, the sample becomes fully connected by the wormhole and the permeability increases with high order. The permeability curve can then be fitted with two power law curves, K = 74ϕ ^2.1 before breakthrough and K = 3.4 × 10¹⁰ ϕ ¹⁵ after breakthrough. The orders of the uniform and wormhole regimes correspond those observed in the literature (Menke et al., 2017).

Similarly, Figure 12B shows the evolution of the correction factor α as a function of porosity during the simulations for the three regimes. For each regime, the evolution of α can be fitted with a power law (Seigneur et al., 2019) and we obtain α = 1.7 × 10⁻⁶ ϕ ^−5.4 for the compact regime, α = 3.1 × 10⁻⁴ ϕ ^−2.8 for the wormholing regime and α = 7.8 × 10⁻³ ϕ ^−2.4 for the uniform regime. We observe that the correction factor for the uniform regime is not constant but decreases with a power-law. This indicates that the dissolution might not be in the uniform regime, but could be in the channeling regime identified in (Menke et al., 2017).

We conclude that the accuracy and computational efficiency of the iVoS method enables simulation of dissolution in a 3D micro-CT image of a real carbonate sample, reproduce experimental results with good precision and can be used to investigate upscaling parameters.