Water–rock interactions can alter rock properties through chemical reactions during subsurface transport processes like geological CO₂ sequestration (GCS), matrix acidizing, and waterflooding in carbonate formations. Dynamic changes in rock properties cause a failure of waterflooding and GCS and could also dramatically affect the efficiency of the acidizing. Efficient numerical simulations are thus essential to the optimized design of those subsurface processes. In this paper, we develop a three-dimensional (3D) numerical model for simulating the coupled processes of fluid flow and chemical reactions in fractured carbonate formations. In the proposed model, we employ the Stokes–Brinkman equation for momentum balance, which is a single-domain formulation for modeling fluid flow in fractured porous media. We then couple the Stokes–Brinkman equation with reactive-transport equations. The model can be formulated to describe linear as well as radial flow. We employ a decoupling procedure that sequentially solves the Stokes–Brinkman equation and the reactive transport equations. Numerical experiments show that the proposed method can model the coupled processes of fluid flow, solute transport, chemical reactions, and alterations of rock properties in both linear and radial flow scenarios. The rock heterogeneity and the mineral volume fractions are two important factors that significantly affect the structure of conductive channels.

Dissolution rates of porous crystalline materials reflect the superposition of transport and surface control, mainly via the parameters saturation of the ambient fluid and distribution of surface energy. As a result, reacting surfaces evolve over time showing a heterogeneous distribution of surface rates. The spatiotemporal heterogeneity of surface reaction rates is analyzed using the rate map and rate spectra concept. Here, we quantify the dissolution rate variability covering the nm- to mm-scale of dissolving single-crystal and polycrystalline calcite samples, using a combined approach of X-ray micro-computed tomography (μ-CT) and vertical scanning interferometry (VSI). The dissolution experiments cover reaction periods from 15 minutes up to 54 days. The observed rate ranges are remarkably consistent over the entire reaction period but include a variability of about two orders of magnitude (10⁻⁹ − 3 × 10⁻⁷ mol m⁻² s⁻¹). The rate map data underscore the concurrent and superimposing impact of surface- vs. fluid flow controlled rate portions. The impact of fluid flow on reactivity at the mm-scale in the transport-controlled system is confirmed by 2-D reactive transport modeling. The sub-mm spatial heterogeneity of low vs. high reactivity surface portions of polycrystalline calcite is clearly below the mean crystal size. This suggests the dominant impact of highly reactive surface portions irrespective of the orientation of larger crystals on the overall surface reactivity. Correspondingly, the overall range of intrinsic reactivity heterogeneity as observed using singly crystal material is not further expanded for polycrystalline material. As a general conclusion, numerical reactive transport concepts would benefit from the implementation of a reactivity term resembling the experimentally observed existence of multiple rate components.

Accurate modeling and simulation of radionuclide migration in clay rocks such as the Opalinus Clay (OPA) play a key role in the safety assessment of deep geological repositories for nuclear wastes. At the continuum scale, the representative elementary volume (REV) is a fundamental constraint to quantify the effective diffusivity, which is a key parameter in reactive transport (RT) models. Therefore, an accurate estimation of the REV is essential for a meaningful continuum‐scale RT simulation in heterogeneous clay rocks. This study presents a comprehensive analysis of the heterogeneities of porosity and effective diffusivity in clay rocks by using the classical sampling theory and pore‐scale simulations. First, in this study, the two‐dimensional representative elementary area is correlated with the REV for porosity via a characteristic length. Next, it is shown that the REV for diffusivity is larger than the REV for porosity. Moreover, these two REVs can be correlated using Archie's law. In such a way, the REV for diffusivity can be determined by the developed correlations through analyzing two‐dimensional microstructures, thus significantly reducing the computational cost. Finally, the applicability of our approach for clay rocks is validated by experimental data on the diffusion of tritiated water in the heterogeneous sandy facies of OPA. From both the experimental data and the modeling prediction, the REV for diffusivity in the sandy facies of OPA is in the order of cubic centimeters. This study provides critical insights into the diffusion in heterogeneous clay rocks toward an enhanced predictability of radionuclide migration.