key: cord-0058878-n2gdc1jm
authors: Lisitsa, Vadim; Khachkova, Tatyana
title: Numerical Simulation of the Reactive Transport at the Pore Scale
date: 2020-08-24
journal: Computational Science and Its Applications - ICCSA 2020
DOI: 10.1007/978-3-030-58799-4_9
sha: ce0028c77b61376ca6beb898c0878e87ce5bb978
doc_id: 58878
cord_uid: n2gdc1jm

We present a numerical algorithm for simulation of the chemical fluid-solid interaction at the pore scale, focusing on the core matrix dissolution and secondary mineral precipitation. The algorithm is based on the explicit use of chemical kinetics of the heterogeneous reactions to determine the rate of the fluid-solid interface changes. Fluid flow and chemical transport are simulated by finite differences. In contrast, level-set methods with immersed boundary conditions are applied to account for arbitrary interface position (not aligned to the grid lines).

Reactive transport in porous media has a number of applications, especially related to the Earth sciences. Among others, it is worth mentioning enhanced oil recovery (EOR), especially for unconventional reservoirs [6] , acid hydrofracturing [38] . This reservoir treatment increases the permeability of the reservoir due to the partial dissolution of the rock matrix and the dissolution channels forming. Also, CO 2 sequestration is, probably, the most well-known example of reactive transport. Typically CO 2 is injected in carbonate formations, where it dissolves in a formation fluid, increasing its acidity; thus, initiating chemical reaction with the rock matrix [16, 35] .

The efficient organization of EOR procedures and CO 2 sequestration requires predictive hydrodynamic modeling accounting for the reactive transport at the reservoir scale. There are a number of models of these processes at the reservoir scale [1, 8, 13, 25, 29] . However, these models use empirical relations between the permeability and tortuosity of the pore space, which determine the convective and diffusion components of the process, respectively, with porosity. Moreover, these relations may vary during the dissolution process. We suggest using numerical simulation directly at the pore scale to recover these relations for a wide range of dissolution scenarios. Mathematical model, describing the reactive transport at the pore scale includes the equations describing the fluid flow, the transfer of chemical components (convection-diffusion) and the reaction at the interface between the solid and liquid phases, taking into account the chemicalkinetic relations [5, 8, 10, 15, 21] . The main difference between the approaches is the approximation of the boundary conditions on a moving boundary. Three different methods can be mentioned. The first is based on the use of rectangular grids and the corresponding finite-difference approximations of the Stokes or Navier-Stokes equations and convection-diffusion. The reaction is taken into account by changing the "concentration" or "relative mass" of the solid phase at the boundary points [14, 37] . The advantage of this approach is the ease of implementation. The disadvantages of the method include the semi-empirical specification of the rate of dissolution of the solid phase (change in "concentration") based on some averaging of the chemical properties of the medium in the vicinity of the boundary. The second group of models is characterized by an explicit treatment of the interface. In this case, the position of the boundary is set explicitly with the application of the finite volume method with truncated cells [23, 24, 30, 33] . However, approaches of this type are difficult to implement when the topology of the computational domain changes. The third group of methods includes those with the implicit treatment of the interface position. Such methods include the level-set [17, 18, 26] and the phase-field method [36] . These approaches allow treating the chemical kinetics of the process explicitly. Moreover, level-set and phase-field methods easily handle topology changes. Such models for solving the Stokes equations and convection-diffusion allow using the finite difference method on a regular rectangular grid in combination with immersed boundary conditions [20, 22, 27, 28, 34] . This paper presents an algorithm for modeling chemical interaction with the rock, based on the level-set method to take into account changes in the position of the boundary using the immersed boundary approach. Section 2 provides a statement of the problem and a mathematical model. The numerical algorithms used to solve the problem are discussed in Sect. 3. The main attention is paid to the approximation of boundary conditions for a complex and changing topology of pore space. Numerical experiments are presented in Sect. 4.

To simulate reactive transport at the pore scale, we use the splitting concerning physical scales. We assume that the interface changes with the slowest rate among all considered processes. Thus, this rate defines the time scale of the problem in general. Also, the fluid flow rate is considered low, which is always a case for underground flows. At the same time, the flow comes to a steady state almost instantly in case of small geometry changes. As a result, three different time scales are considered in the model separately.

We consider a bounded domain D ⊆ R 2 , composed of the nonintersecting time-dependent subdomains D p (t) and D m (t) corresponding the pore space and the matrix, respectively. The boundary of the domain is ∂D = Γ outlet ∪ Γ inlet ∪ Γ nf . We denote the interface between the pore space and the matrix asD p (t) ∩ D m (t) = Γ (t) which is a union of the sufficiently smooth curves. A sketch of hte model is presented in Fig. 1 .

The fluid flow is defined in the pore space D p (t) satisfying the steady state Stokes equation:

with boundary conditions

where u = (u 1 , u 2 ) T ∈ R 2 is the velocity vector, p is the pressure, μ is the dynamic viscosity, p bc (x) is the pressure at the inlet boundary, x = (x 1 , x 2 ) T is the vector of spatial coordinates. The active component propagates in the pore space D p (t) according to the convection-diffusion equation:

with boundary conditions:

In these notations D is the diffusion coefficient, C is the reactant concentration, n is the inner (with respect to D p (t)) normal vector, k r is the reaction rate coefficient, Γ inlet is the inlet boundary. We consider only the first order reactions, one active reactant and no reactant sources inside the computational domain; thus the right-hand side in (3) equals to zero.

The changes in the interface Γ (t) caused by the chemical reactions governed by the equation

where v n is the normal component of the boundary velocity, K c is the stoichiometric coefficient, ρ is the mass density of the matrix mineral, C s reactant concentration at equilibrium.

We use the finite differences to solve Eqs. (1) and (3). However, the interface Γ (t) vary and may not align to the grid lines. To deal with such an irregular geometry of the interface, we use the level-set method where the interface Γ (t) is defined implicitly as the constant level line of a function ϕ(x, t); i.e.

It is natural to require Γ (t) as the zero-level. In this case, the subdomains D p and D m are defined as

Moreover, the level-set function ϕ(x) is constructed as the signed distance to the interface; i.e. ∇ x ϕ(x, t) = 1. This leads to a natural definition of the normal vector

Using the level-set function one may rewrite the equation for the interface changes as follows [9, 26] :

where v n is the normal velocity, defined by the Eq. (5).

To solve Eqs. (1), (3), and (6) we introduce the staggered grid, where the gridfunctions are defied according to the rule

where h 1 and h 2 are the spatial steps, τ is the grid step with respect to time (either physical or artificial). Note, that the interface Γ (t) defined as a constant level of function ϕ is approximated by a piecewise linear line. Note also that the model is a micro-CT scan of the real rock sample [2] [3] [4] ; thus, it is extremely detailed for either description of the pore space structure or the physical processes. As a result the low-order finite difference schemes can be used in application to pore scale simulation.

We use the second order finite-difference scheme to approximate the Stokes equation, where pressure is defined in the centers of the grid cells; i.e. in points with integer indices; whereas the components of the velocity vector are defined at the faces of the cells; i.e., in points with one half-integer index:

In these notations D c 1 and D 2 2 are the central differences approximating first derivatives with respect to spatial variables; whereas D 2 1 and D 2 2 are the second order approximations of the second derivatives:

Here indices I and J can be either integer or half-integer, but i and j are integers only. The operators approximating derivatives with respect to the other spatial direction can be obtained by the permutations of the role of spatial indices. Note that level-set function ϕ i,j is defined at integer points, and its values in half-integer points are interpolated. We are not specifying the iterative method to solve system (1), but focus on the computation of the action of the system matrix on a vector. Equations (1) are valid only if all the points from a stencil, centered in particular point belong to pore space D p (t). If a stencil contains points from D m (t), where the solution is not defined we suggest using the immersed boundary method [12, 17, 19, 20, 22, 28, 34] . The idea of this approach is to extrapolate the solution from D p (t) to the points from D m (t) which are needed for simulation. The extrapolation is based on the boundary conditions stated at the interface;

i.e., u = 0 on Γ (t); thus, velocity vector can be considered as odd function with respect to the interface, whereas p is even. Assume one needs to extrapolate the solution to the point ((x 1 ) I , (x 2 ) J ). The distance from this point to the interface is defined by the level-set function ϕ I,J and the normal direction to the interface is n = ∇ x ϕ(x)| I,J . Thus, the projection of the point to the interface is

and the orthogonal reflection has the coordinates ((x 1 ) n , (x 2 ) n ) = ((x 1 ) I , (x 2 ) J ) + 2ϕ I,J n, and moreover u I,J = −u((x 1 ) n , (x 2 ) n ) + O(ϕ 2 I,J ). A standard approach to compute u((x 1 ) n , (x 2 ) n ) is the use of bilinear interpolation using the four nearest grid points from the regular grid. However, if a pore space is considered, it is not guaranteed that the four nearest points belong to the pore space D p . Thus, we use all available points to construct the interpolation. We pick the four regular points (those of them which belong to the pore space) and also the point at the interface ((x 1 ) c , (x 2 ) c ). It is clear that to compute linear interpolation; three points are enough; thus, if more than three points are available, the problem is overdetermined, and we compute the minimal-norm solution. If less than three points are available, we can only achieve the first order of accuracy. We can use the interpolation weights proportional to the distances from the considered points to point ((x 1 ) n , (x 2 ) n ).

To approximate the convection-diffusion equation we use the first-order scheme

where

Operators, approximating the derivatives with respect to the other spatial direction, can be obtained by the permutation of the spatial indices. Operators D 2 1 and D 2 2 are introduced above in Eq. (8) . Note that due to very fine discretization, the first-order scheme provides suitable accuracy.

Same as before, the hardest part of the approximation is the dealing with the boundary conditions at Γ (t). The concentration satisfies Robin boundary conditions on Γ (t), which can be considered a linear combination of the Dirichlet and Neumann ones to apply the immersed boundary conditions

We assume that the point ((x 1 ) I , (x 2 ) J ) ∈ D m belongs to a stencil centered in D p . Thus an extrapolation is required. Same as before, ((x 1 ) n , (x 2 ) n ) is the orthogonal reflection of the immersed point ((x 1 ) I , (x 2 ) J ) over the interface; thus, the boundary condition can be approximated as follows

where C((x 1 ) n , (x 2 ) n ) is interpolated using available points in the domain D p (t).

To simulate the interface movement one needs to solve the Eq. (6), which can be done by the finite difference method:

where (v n ) n i,j is the rate of the interface changes. Equation (6) is defined everywhere in D, thus its right-hand side should be defined accordingly. However, function v n (x) is defined only at Γ and should be continued inside subdomains D p and D m . Following [7, 26] the rate v n can be continued inside subdomains as a constant along the normal direction to the interface; i.e., the steady-state solution of the following equation should be computed:

whereṽ n are the initial conditions which coincide with the velocity at the interface and trivial elsewhere. We use the cental differences to approximate the gradient of ϕ whereas the upwind scheme is used to approximate the derivatives of q [7] . Note, that accurate solution of Eq. (13) is only needed in a vicinity of the interface; i.e., inside a strip D s = {x : |ϕ(x)| ≤ 2 h 2 1 + h 2 2 }. Thus, only a few iterations are sufficient to get accurate solution there. After that additional redistancing is applied to ensure that ϕ is the signed distance [9, 26, 31, 32] .

We use the developed algorithm to simulate chemical fluid-solid interaction. We apply statistical simulation to generate the models. In particular, we used the truncated Gaussian field [11] , where the model is represented as a realization of the Gaussian random field G(v, l) with prescribed correlation length l, and point x 0 belongs to D p , if values of the field in this point exceeds a threshold G(x 0 , l) > R. Otherwise, the point belongs to the matrix. Porosity or the ratio of the pore space to the whole volume is defined uniquely by the threshold value R. We consider models with a fixed correlation length of 5 · 10 −5 m, which is the typical size of the heterogeneities of carbonate rocks. The porosity is fixed as 65%, which is too high for realistic rocks but ensures the percolation of the pore space in 2D. The size of the model is [0, 50l 1 ] × [0, 50l 2 ], where l 1 and l 2 are the correlation lengths along x 1 and x 2 respectively. We consider the models where l 1 = l 2 . The grid step is h 1 = h 2 = 10 −5 m.

The core matrix is calcite with mass density ρ = 2710 kg/m 3 , and stoichiometric coefficient of the reaction equal to one; i.e., K = 1. The fluid is the reservoir water under assumption that changes in the reactant concentration do not affect the physical properties of the fluid, thus we fix the dynamic viscosity μ = 0.00028 Pa·s. To vary dissolution scenarios we change the pressure drop at the outer boundary; i.e., p bc in Eq. (2); diffusion coefficient D and the reaction rate k r in Eqs. (3) and (4) . The values of the parameters are provided in the Table 1 . The active component is the cations H + with equilibrium concentration corresponding to pH = 7, whereas acidity at the inlet was pH = 2.3 following laboratory experiments presented in [16] .

It is convenient to use dimensionless variables to describe the results of the simulation. These variables are Reynolds number Re = LU ν ; Peclet number P e = UL D , indicating if convection of diffusion dominates in the reactant transport; Damkohler number Da = krL D , indicating if the reaction is diffusion-limited or kinetic limited. In these notations, U is characteristic flow velocity, L is a characteristic length scale, we assume it to be equal to the correlation length of the Gaussian field. During chemical fluid-solid interaction, the matrix of the core dissolved. Thus both L and U may change. In our further description, we will refer to the values corresponding to the initial time.

To illustrate the applicability of the algorithm for simulation of the chemical fluid-solid interaction, we present the results of six numerical experiments corresponding different combinations of the Reynolds, Peclet, and Damkohler numbers, as presented in Table 1 . Figures 2, 3 and 4 represent the distribution of the pH at different time instants for six experiments. Localized dark spots correspond to the core matrix. Models 1,3,5 correspond to low Damkohler numbers (diffusion-controlled reactions) if compared with models 2,4,6. The diffusioncontrolled reactions are slower, and the reactant is delivered in all channels of the pore space, but its access to the interface is limited. Note that the increase of the Reynolds number causes deeper penetration of the reactant to the sample (experiment 5). If the kinetic-controlled reactions are considered (experiments 2,4,6), the reactant can not penetrate the sample and interacts with the matrix close to the inlet regardless of the values of the Re and P e numbers. 5 .3 · 10 −3 8 · 10 −2 5 · 10 −2 9.3 · 10 −8 1.5 · 10 −4 2 1.2 · 10 −6 5.3 · 10 −3 8 · 10 0 5 · 10 −2 9.3 · 10 −8 1.5 · 10 −2 3 1.2 · 10 −5 5.3 · 10 −3 8 · 10 −2 5 · 10 −1 9.3 · 10 −7 1.5 · 10 −3 4 1.2 · 10 −5 5.3 · 10 −3 8 · 10 0 5 · 10 −1 9.3 · 10 −7 1.5 · 10 −1 5 1.2 · 10 −5 5.3 · 10 −2 8 · 10 −2 5 · 10 −1 9.3 · 10 −8 1.5 · 10 −4 6 1.2 · 10 −5 5.3 · 10 −2 8 · 10 0 5 · 10 −1 9.3 · 10 −8 1.5 · 10 −2

We presented a numerical algorithm for the simulation of chemical fluid-solid interaction at the pore scale. The algorithm is based on splitting with respect to the physical processes. It is assumed that the fluid flow gets steady instantly after small changes of the pore space geometry. Moreover, the flow velocity is low; thus, flow is simulated as a solution to the Stokes equation. The convectiondiffusion equation governs chemical transport with Robin boundary conditions. The pore space -matrix interface is defined implicitly by the level-set method where the boundary conditions are approximated using the immersed boundary method. Presented numerical experiments illustrate the applicability of the suggested approach for simulation chemical fluid-solid interaction for a wide range of models, including different types of rocks and artificial porous materials.

Modelling and simulation of reactive transport phenomena

Digital rock physics benchmarks -part I: imaging and segmentation

Digital rock physics benchmarks -part II: computing effective properties

Effect of CT image size and resolution on the accuracy of rock property estimates

A numerical and analytical study on calcite dissolution and gypsum precipitation

Geochemical monitoring of fluid-rock interaction and CO2 storage at the weyburn CO2-injection enhanced oil recovery site

A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)

Carbonate acidizing: modeling, analysis, and characterization of wormhole formation and propagation

A review of level-set methods and some recent applications

CO2-induced dissolution of low permeability carbonates. Part II: numerical modeling of experiments

Stochastic generation of explicit pore structures by thresholding Gaussian random fields

A cartesian grid embedded boundary method for poisson's equation on irregular domains

Effect of medium heterogeneities on reactive dissolution of carbonates

Pore-scale study of dissolution-induced changes in permeability and porosity of porous media

A robust and efficient numerical method for multiphase equilibrium calculations: application to CO2-brine-rock systems at high temperatures, pressures and salinities

Carbon geosequestration in limestone: pore-scale dissolution and geomechanical weakening

Level set simulation of coupled advection-diffusion and pore structure evolution due to mineral precipitation in porous media

A three-dimensional level set simulation of coupled reactive transport and precipitation/dissolution

A ghost-cell immersed boundary method for simulations of heat transfer in compressible flows under different boundary conditions

Sharp interface cartesian grid method i: an easily implemented technique for 3D moving boundary computations

Mesoscopic dynamics of solid-liquid interfaces: a general mathematical model

Immersed boundary methods

An investigation of the effect of pore scale flow on average geochemical reaction rates using direct numerical simulation

Pore-scale controls on calcite dissolution rates from flow-through laboratory and numerical experiments

Modeling acid leakoff during multistage alternate injection of pad and acid in acid fracturing

Level set methods: an overview and some recent results

Flow patterns around heart valves: a numerical method

Immersed boundary methods for simulating fluidstructure interaction

Reactive transport codes for subsurface environmental simulation

A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems

An efficient, interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow

An improved level set method for incompressible two-phase flows

High-resolution simulation of pore-scale reactive transport processes associated with carbon sequestration

A ghost-cell immersed boundary method for flow in complex geometry

Rock physics analysis and time-lapse rock imaging of geochemical effects due to the injection of CO2 into reservoir rocks

Phase-field modeling of solute precipitation and dissolution

Pore-scale simulation of mixinginduced calcium carbonate precipitation and dissolution in a microfluidic pore network

Rock specific hydraulic fracturing and matrix acidizing to enhance a geothermal system -concepts and field results