A Locally Conservative Mixed Finite Element Framework for Coupled Hydro-Mechanical-Chemical Processes in Heterogeneous Porous Media

This paper presents a mixed finite element framework for coupled hydro-mechanical-chemical processes in heterogeneous porous media. The framework combines two types of locally conservative discretization schemes: (1) an enriched Galerkin method for reactive flow, and (2) a three-field mixed finite element method for coupled fluid flow and solid deformation. This combination ensures local mass conservation, which is critical to flow and transport in heterogeneous porous media, with a relatively affordable computational cost. A particular class of the framework is constructed for calcite precipitation/dissolution reactions, incorporating their nonlinear effects on the fluid viscosity and solid deformation. Linearization schemes and algorithms for solving the nonlinear algebraic system are also presented. Through numerical examples of various complexity, we demonstrate that the proposed framework is a robust and efficient computational method for simulation of reactive flow and transport in deformable porous media, even when the material properties are strongly heterogeneous and anisotropic.


Introduction
Hydro-mechanical-chemical (HMC) processes in porous media, in which fluid flow, solid deformation, and chemical reactions are tightly coupled, appear in a variety of problems 1. To formulate a robust numerical approximation scheme for coupled HMC processes in heterogeneous porous media, employing a combination of locally conservative finite element methods. 2. To reduce the computational cost for solving an advection-diffusion-reaction equation by using the EG method, which requires approximately two and three times fewer degrees of freedom than the DG method for 2D and 3D geometries, respectively [58]. 3. To demonstrate the performance and capabilities of the proposed framework for modeling tightly coupled HMC problems with homogeneous to heterogeneous, isotropic to anisotropic permeability fields with local conservation.
The rest of the paper is organized as follows. Section 2 describes the governing equations of coupled HMC processes. Section 3 explains the discretization methods, linearization techniques, and solution algorithms of the proposed framework. Section 4 presents several numerical examples of various complexity and discusses key points found in this paper. Section 5 concludes the work.

Governing equations
This section briefly describes all the equations used in this study, namely poroelasticity and advection-diffusion-reaction equations.

Poroelasticity
To begin, we adopt Biot's poroelasticity theory for coupled hydro-mechanical processes in porous media [59,60]. Although poroelasticity may oversimplify deformations in soft porous materials such as soils [61][62][63][64], it would be reasonably good for stiff materials such as rocks, which is the focus of this work. The poroelasticity theory provides two coupled governing equations, namely linear momentum and mass balance equations. Under quasistatic conditions, the linear momentum balance equation can be written as where f is the body force term defined as ρφg + ρ s (1 − φ)g, where ρ is the fluid density, ρ s is the solid density, φ is the porosity, g is the gravitational acceleration vector. The gravitational force will be neglected in this study, but the body force term will be kept in the succeeding formulations for a more general case. Further, σ is the total stress tensor, which may be related to the effective stress tensor σ and the pore pressure p as σ(u, p) = σ (u) − αpI.
Here, I is the second-order identity tensor, and α is the Biot coefficient defined as [65]: with K and K s being the bulk moduli of the solid matrix and the solid grain, respectively. According to linear elasticity, the effective stress tensor has a constitutive relationship with the displacement vector, which can be written as σ (u) = λ l tr(ε(u))I + 2µ l ε(u).
For this solid deformation problem, the domain boundary ∂Ω is assumed to be suitably decomposed into displacement and traction boundaries, ∂Ω u and ∂Ω t , respectively. Then the linear momentum balance equation is supplemented by the boundary and initial conditions as: ∇ · σ (u) + α∇ · (pI) + f = 0 in Ω × T, u = u D on ∂Ω u × T, σ(u) · n = t D on ∂Ω t × T, u = u 0 in Ω at t = 0, where u D and t D are prescribed displacement and traction values at the boundaries, respectively, and n is the unit normal vector to the boundary. Next, the mass balance equation is given as [4,19,43,66]: is the Biot modulus.
Here, c f is the fluid compressibility, φ 0 is the initial porosity, ε v := tr(ε) = ∇ · u is the volumetric strain, and g is a sink/source term. Because we will introduce chemical effects later on, we have added ∂φc ∂t to the standard poroelasticity equation [3,4,43,67]. This term will be discussed again after introducing chemical effects. Also, q is the superficial velocity vector, which is given by Darcy's law as Note that here the fluid viscosity µ is considered a function of concentration c i . Again, the gravitational force, ρg, will be neglected in this work, without loss of generality. In addition, k(φ) is the matrix permeability tensor defined as The k xx , k yy , and k zz represent the matrix permeability in x-, y-, and z-direction, respectively. The k mult (φ) is a multiplier used to update k when φ is altered, which will be described later.
For the fluid flow problem, the domain boundary ∂Ω is also suitably decomposed into the pressure and flux boundaries, ∂Ω p and ∂Ω q , respectively. In what follows, we apply the fixed stress split scheme [19,68], assuming (σ v − σ v,0 ) + α (p − p 0 ) = Kε v . Then we write the fluid flow problem with boundary and initial conditions as where σ v := 1 3 tr(σ) is the volumetric stress, and p D and q D are the given boundary pressure and flux, respectively.

Reactive flow
An advection-diffusion-reaction system for N c number of the miscible species is given by the following equations. For all i = 1, . . . , N c , where q i (c i ) is a reaction term coupled with sink/source for each component, and the mass flux η(q, c i ) is defined as Here D e,i (φ) is the effective diffusion coefficient tensor defined as where τ = φ − 1 2 [69,70] and D i is the given diffusion coefficient tensor. The boundary for the advection-diffusion-reaction system is decomposed into inflow and outflow boundaries, denoted by ∂Ω in and ∂Ω out , respectively, which are defined as ∂Ω in := {x ∈ ∂Ω : q · n < 0} and ∂Ω out := {x ∈ ∂Ω : q · n ≥ 0}. (16) In what follows, we specialize the model to calcite precipitation and dissolution reactions, which requires us to solve a calcite-carbonic acid system. In general, the system requires eight transport equations to solve the concentration values of the following main species/ions: [3,67,71,72]. For simplicity, in this paper we consider a reduced system based on the empirical relationship presented in [3,4,67], in which N c decreases to 1. Thus, letting c := c 1 , we write the advection-diffusion-reaction system with its boundary and initial conditions as follows: 6 where c in is the inflow concentration, c 0 is the initial concentration, and q represents a source term reflecting the calcite dissolution/precipitation reactions. For this term, here we adopt the term in [3,4,67], given by where A s is the specific surface of the porous medium, and R c is the reaction rate calculated as with and Here, τ is the medium temperature, and a 0 , a 1 , a 2 , a 3 , a 4 , and a 5 are defined in Table 1. Before closing this section, we describe physical properties that are coupled with primary variables, u, q, p, and c. The porosity change due to solid deformation may be expressed as [18,59,73]: where v 0 is the initial volumetric strain. The porosity alteration due to calcite dissolution/precipitation is calculated as where ω is the number of moles of total precipitated species per kilogram of rock (assumed to be 10.0 in this study following [3,4]), and ρ s = 2500 kg/m 3 is used for throughout this paper.
Note that this term, (24), enters (12). Also, the terms φ m and φ c are used to distinguish between the changes in φ due to solid deformation as in (23), and chemical reactions as in (24), respectively. The changes in porosity due to (23) and (24), also affect the specific surface for porous medium (A s ) as where A 0 is the initial value of A s , and it is set as 5000 throughout this study [74]. Furthermore, the porosity change influences the matrix permeability as [75][76][77]: where k 0 is the initial matrix permeability and b is an empirical parameter determined experimentally. In this work, we set b = 22.2 following [75]. The change in c also affects µ, and we adopt the specific form from [78][79][80], given by where c l and c h are lower and higher bounds of the concentration, and µ l and µ h are fluid viscosity corresponding to c l and c h , respectively. Table 2 summarizes the effects of physical processes on material properties considered in this study. Note that the numbers, e.g., (23), point out the equations used to represent these effects, while a hyphen means the absence of a relationship.

Numerical methods
In this section, we describe the numerical methods for the governing system described in the previous sections. Here, we utilize a combination of a mixed finite element method for spatial discretization, and employ both a backward differentiation formula and an explicit Runge-Kutta method for temporal discretization.

Domain discretization and geometrical quantities
We begin by introducing the notations used throughout this paper. Let T h be a shaperegular triangulation obtained by a partition of Ω into d-simplices (triangles in d = 2, tetrahedra in d = 3). For each cell T ∈ T h , we denote by h T the diameter of T , and we set h = max T ∈T h h T and h l = min T ∈T h h T . We further denote by E h the set of all faces (i.e., d − 1 dimensional entities connected to at least a T ∈ T h ) and by E I h and E ∂ h the collection of all interior and boundary facets, respectively. The boundary set E ∂ h is decomposed into two disjoint subsets associated with the Dirichlet boundary faces, and the Neumann boundary faces for each of (7) and (12). In particular, E D,u h and E N,u h correspond to the faces on ∂Ω u and ∂Ω tr , respectively, for (7). On the other hand, for (12), E D,m h and E N,m h conform to ∂Ω p and ∂Ω q , respectively. Lastly, for (17), E ∂ h is decomposed into E In h and E Out h . We also define e = ∂T + ∩ ∂T − , e ∈ E I h , where T + and T − are the two neighboring elements to e. We denote by h e the characteristic length of e calculated as h e := meas (T + ) + meas (T − ) 2 meas(e) , depending on the argument, meas(·) represents the measure of a cell or of a facet. Let n + and n − be the outward unit normal vectors to ∂T + and ∂T − , respectively. For any given scalar function ζ : T h → R and vector function τ : T h → R d , we denote by ζ ± and τ ± the restrictions of ζ and τ to T ± , respectively. Subsequently, we define the weighted average operator as and where δ e is calculated by [81,82]: Here, k + e := n + · k + n + , and where k e is a harmonic average of k + e and k − e which reads and k is defined as in (11). The jump across an interior edge will be defined as

Temporal discretization
The time domain T = (0, T] is partitioned into N subintervals such that 0 =: t 0 < t 1 < · · · < t N := T. The length of each subinterval ∆t n−1 is defined as ∆t n−1 = t n − t n−1 where n represents the current time step. We assume that the user provides the initial ∆t 0 , while an adaptive procedure is carried out to choose ∆t n−1 , n > 1, as follows: where CFL is a constant that the user can provide according to the Courant-Friedrichs-Lewy condition [83], · ∞ is the maximum norm of a vector function, and ∆t max is a maximum allowed time step. Note that we use ∆t max as a tool to control ∆t n as the model approaches a steady-state condition since q n−1 ∞ may approach zero, which would lead to a very large ratio Let ϕ(·, t) be a scalar function and ϕ n be its approximation at time t n , i.e. ϕ n ≈ ϕ (t n ). We employ the following backward differentiation formula [84][85][86] BDF m (ϕ n ) := for the discretization of the time derivative of ϕ(·, t) at time t n . We also utilize the explicit Runge-Kutta methods [24,87]: for the first order Runge-Kutta method corresponding to the explicit Euler method, and for the forth order Runge-Kutta method, F (X n , Y n ) is any functions with independent variable X and dependent variable Y [8,87], which we will specify in the linearization and solving processes in Section 3.4.
Finally, we define an extrapolation operator as follows [8,24]: and in the following we will adopt the notationφ n+1 to denote an extrapolation value of {ϕ n , ϕ n−1 }.

Spatial discretization
In this framework, the displacement field is approximated by the classical continuous Galerkin method (CG) method, and the fluid velocity and pressure fields are discretized by the Brezzi-Douglas-Marini (BDM) element [88] and the piecewise constants discontinuous Galerkin (DG) method, respectively, to ensure local mass conservation. Lastly, the concentration field is discretized by the enriched Galerkin (EG) method [46,49].
To begin, we define the finite element space for the CG function space for a vector-valued function: where C 0 (Ω; R d ) denotes the space of vector-valued piecewise continuous polynomials, P k (T ; R d ) is the space of polynomials of degree at most k over each element T , and ψ u denotes a generic function of U CG k h (T h ). In addition, the CG space for scalar-valued functions is defined as: where C 0 (Ω) := C 0 (Ω; R) and P k (T ) := P k (T ; R). Next, we define the following DG function space: where L 2 (Ω) is the space of square-integrable scalar functions. This non-conforming finite element space allows us to consider discontinuous functions and coefficients rigorously. We then define the EG finite element space with polynomial order k as: i.e., a CG finite element space enriched by the space P DG 0 h (T h ) of piecewise constant functions. In the following we denote ψ c a generic function of P EG k h (T h ). Lastly, we define the BDM function space as follows [88]: where ψ v denotes a generic function of V BDM k h (T h ) and BDM(T ) is defined according to [88].

Fully discrete form
We now present the fully discrete form of the coupled HMC problem using the abovedescribed combination of finite element spaces. In particular, we seek the approximated displacement solution u h ∈ U CG 2 h (T h ) as done in [20,51,89] We multiply the linear momentum balance equation (7) by a test function ψ u ∈ U CG 2 h (T h ). The fully discretized linear momentum balance equation thus has the following form: at each time step t n , where Here T · dV and e · dS refer to volume and surface integrals, respectively, and ∇ s is the symmetric gradient operator. Furthermore, the notation for N u (ψ u ; u n h , p n h ) in (44) highlights before the semicolon the test function, and after the semicolon the (possibly nonlinear) dependence on discrete solutions to the coupled problem. The same notation will be used hereafter for the remaining equations.
Next, the weak form of the mass balance equation (12) is obtained multiplying by ψ p ∈ P DG 0 h (T h ) and integrating by parts, resulting in: for each time step t n , where For the Darcy velocity equation (10), we obtain where Lastly, for the advection-diffusion-reaction equations of species transport we write: for each time step t n , where We note that the D e as defined in [90][91][92]. The γh q n h I term is often referred as the first order artificial diffusivity coefficient [93,94]. In our paper, we set the tuning parameter γ = 0.25. Alternative stabilization strategies including streamline diffusion and crosswind diffusion, or entropy viscosity methods could be also utilized to reduce oscillations in the numerical solution to the concentration field [90][91][92][94][95][96][97][98].
Also, c up h is an upwind value of c n h defined as [45,99]: where c n+ h and c n− h correspond to c n h of T + and T − , respectively. Lastly, the two parameters θ and β define corresponding interior penalty methods. The discretization becomes the symmetric interior penalty Galerkin method (SIPG) when θ = −1, the incomplete interior penalty Galerkin method (IIPG) when θ = 0, and the nonsymmetric interior penalty Galerkin method (NIPG) when θ = 1 [45]. In this study, we set θ = −1 for the simplicity and β = 1.1 throughout this paper.
Remark 1. For the momentum balance equation (7), the traction boundary condition t D (traction) is applied weakly on each e ∈ E N,u h in (44), while the displacement boundary condition u D is strongly enforced on each e ∈ E D,u h . For the mass balance equation (12), since we use a mixed formulation, the flux boundary condition q D is strongly applied on each e ∈ E N,m h , but the pressure boundary condition p D is weakly applied on each E D,m h in (46). Finally, for the transport equation (17), all boundary conditions are weakly applied in (47).

Remark 2.
In our computational framework, we provide a flexible choice of the time discretization schemes for each equation. We use BDF 1 for the time discretization of the mass balance equation (12) since it is sufficient to provide the optimal error convergence rate, see [17]. For the time discretization of the transport equation (17), we use BDF 4 to capture a sharp front in the advection dominated regime [45].

Splitting algorithm
The coupled system obtained from the discrete governing equations (44), (45), (46), and (47) is nonlinear. Although the coupled nonlinear system may be solved in a monolithic manner, here we focus on developing a splitting algorithm for sequential solution to the coupled system, which can provide more flexibility especially when different software packages need to be combined. The overall computational strategy is summarized in Algorithm 1. for each fixed stress iteration step (·) n,ι until δφ n,ι < TOL do 9: Solve (45) and (46)  Calculate φ n,ι m see (23) 13: Evaluate F (σ v n,ι ) see (52) 14: Evaluate δφ n,ι see (50) 15: Update k n,ι see (23), (26) 16: end for 17: In Algorithm 1, we separate our algorithm into two parts. The first part (lines 8 to 17) focuses on solving the coupled hydro-mechanical problem, (44), (45), and (46), using the fixed stress method which is an unconditionally stable splitting scheme [18,19,68,73]. At each iteration ι we solve (45) and (46) for the velocity q n,ι h and the pressure p n,ι h using a monolithic method (line 9) based on given displacement u n,ι−1 h from previous nonlinear iteration and concentration extrapolatedĉ n h from previous time step. Then, we couple with (44) using the fixed-stress split scheme based on the pressure p n,ι h computed at the current nonlinear iteration (line 11). The convergence criterion is based on δφ n,ι (Algorithm 1 line 8), which is defined as: 15 Here, φ n,ι m is the porosity resulting from the solid deformation (23) and φ n,ι f is the porosity resulting from the fluid flow problem defined as [18,68,73]: where (·) ι represents iteration counter inside the fixed-stress loop. From the fixed stress split concept (51) is the φ predictor, while (23) is the φ corrector [18,19,68,73]. Hence, when φ n,ι m and φ n,ι f converge, i.e., δφ n,ι < TOL, the fixed-stress loop is completed. The tolerance TOL is set as 1 × 10 −6 throughout this study. Note that the flow equations, (45) and (46), are solved by assuming that ∂σv ∂t = 0, i.e., F (σ v n,ι ) is frozen; therefore, this term is evaluated explicitly after the momentum equation (44) is solved, as illustrated in Algorithm 1 line 13 [19,68], and F (σ v n,ι ) is defined as: The second part (from line 18) focuses on solving advection-diffusion-reaction equation (47), using q n h , φ n , D n e , and A n s obtained from the first part. One could view this strategy as a one-way coupling scheme between coupled hydro-mechanical and advection-diffusionreaction equations. Next where (·) n represents an extrapolation value based on the extrapolation described in (38). Subsequently, we evaluate F φ c n,ι and F (q n,ι ), which are defined as and F q n+1 := usingĉ n+1 h calculated by (53). We note that the equation (47) becomes linear by employinĝ c n+1 h to calculate F (q). Also, the porosity alteration as a result of calcite dissolution/precipitation (Algorithm 1 line 24) is computed bŷ Note that the porosity change due to the calcite dissolution/precipitation reactions is additional to the porosity change by solid deformation, (23). Subsequently,k n+1 ,D n+1 t , and A s n+1 are determined usingφ n+1 . Lastly, we also calculateμ n+1 usingĉ n+1 h , see (53) and (27).
For all the computations, matrices and vectors are built using the FEniCS form compiler [100]. The block structure is assembled by using the multiphenics toolbox [101]. Solvers are employed from PETSc package [102]. All simulations are computed on XeonE5 2650v4 with a single thread.
Remark 3. We note that the EG method, which is used to approximate the advectiondiffusion-reaction (17), is based on the Galerkin method, which could be extended to consider adaptive meshes that contain hanging nodes. Besides, an adaptive enrichment, i.e., the piecewise-constant functions only added to the elements where the sharp material discontinuities are observed, can be developed.

Numerical examples
In this section, we demonstrate the performance and capabilities of the proposed numerical method through various numerical examples. We begin with a single-layer model comparing the performance for single-phase flow with chemical dissolution/precipitation and solid deformation. Then we illustrate the performance of the developed model for a layered medium as well as a heterogeneous single-layer medium. Lastly, we test the proposed framework using an example with an anisotropic permeability field. All four examples and their mesh are illustrated in Figure 1. More detailed setup, including the input parameters and the boundary conditions of each example, are described in the beginning of each example.  (7), we assume u D · n = 0 on ∂Ω 1 , ∂Ω 3 , and ∂Ω 4 . Furthermore, t D = [0.0, −2.0 × 10 6 ] Pa is applied on ∂Ω 2 . Therefore, the medium is under compression. For the mass balance equation (12), the boundary condition q D = 0 is set on ∂Ω 2 and ∂Ω 4 and we impose p D = 1 × 10 5 Pa on ∂Ω 3 . Here, for the mass balance equation (12), we test two different scenarios on ∂Ω 1 , where scenario (a) corresponds to q D = 2 × 10 −4 m/s and scenario (b) is characterized by q D = 1 × 10 −4 m/s. Thus, scenarios (a) and (b) will be referred to as the high and low injection rate cases, respectively. Since we want to compare the results of the above scenarios at the same total injected volume (I.V.), which is defined as where A d = 30m 2 is the surface area of ∂Ω 1 , the time t n of the scenario (b) is twice to scenario (a). For the advection-diffusion-reaction equation (17), we impose the inflow condition c in = 0.5 on ∂Ω 1 . The initial pressure p 0 is 1 × 10 6 Pa, the initial concentration c 0 is calculated by (22) using p = p 0 and τ = 20 C, and the initial displacement u 0 is calculated as stated in Algorithm 1. The penalty parameter (β) is set to be 1.1 for the EG method. The CFL is used as 0.1 for calculating ∆t n , see (34).
Here, we compare the transient distribution of the concentration achieved with the developed HMC coupled numerical scheme in a homogeneous porous medium for two different injection rates. The aim is to illustrate the impact of different processes on the advance of the flow path and reactive solute transport. Initially, the composition of the pore fluid within the porous medium is in equilibrium with calcite. Note that c eq calculated by (22) is a function of temperature and pressure. In this example, assuming constant temperature, pressure deviates from the initial fluid pressure in time and space. The changes in the pressure field as a result of fluid injection on the left boundary and the fluid production on the right boundary varies the c eq resulting in precipitation or dissolution in the domain. The injected water is also unsaturated with respect to calcite. Therefore, the injected fluid, as advances into the domain, will dissolve the calcite mineral. Figure 2 shows the concentration fields at different injected fluid volumes (I.V.) and for both scenarios associated to q D . There are three main observations from these figures. The first one is the flow instability, or fingering, emerged as a result of the difference between the injected fluid viscosity and the in-situ fluid viscosity. The second observation is that for the higher injection rate scenario, the fingers are more developed at a later time compared to that of the low injection scenario. The third one is that most of the fingers developed initially either merge or vanishes at the later stage, forming one or two main fingers. Next, we present the interaction among different processes including mechanical deformation, calcite dissolution/precipitation, and viscosity alteration in Figure 3 for two different time steps. Note that the results of the low injection rate case are similar (for the same volume of injected fluid); hence, we present here only the results of the high injection rate case. First, one could observe that the effect of mechanical deformation is dictated by both p and u, see Figure 3b and g. Figure 3a and f illustrate the reduction of φ by the solid deformation as the model is under compression. The increased fluid pressure by fluid injection, however, limits the porosity reduction. This is reflected in Figures 3a and f  The ∂φc ∂t result is shown in Figure 3c and h. Since the injected concentration c in = 0.5 is lower than c 0 (initial c eq ), the porous medium is dissolved in places to which the injected fluid is transported. Note that ∂φc ∂t is positive where the dissolution occurs and negative where the perception occurs. At this time step, the maximum magnitude of ∂φc ∂t is 10 −7 , which is much less compared to that of ∂φm ∂t , which is around 10 −4 . We note this magnitude could be varied with different input parameters and boundary conditions of each equation, (7), (12), or (17). The value of µ is also altered, see Figure 3d and i, as the concentration front progresses. This alteration causes the flow instability discussed previously and establishes a preferential flow path. The impact of ∂φm ∂t , ∂φc ∂t , and µ alteration can be seen in q field shown in Figure 3e and j. Interestingly, as the first finger reaches the outlet boundary the second finger gradually disappears resulting in only one preferential path between the inlet and outlet of the model. Thus, we have confirmed that the proposed framework can well simulate the expected physical and chemical phenomena including solid deformation, viscous fingering, and dissolution/precipitation. The key ingredients of this method are the capability for tracking the interface of the concentration species approximated by the high order methods with numerical stabilization, the computation of reaction terms with the EG method, and the locally conservative flux from BDM.
The concentration field c for two different injection scenarios (as discussed in example 1) for the three-layer porous medium are presented in Figure 4. Unlike the single-layer porous medium, even though the concentration fields at the early time are similar between the high, q D = 2 × 10 −4 m/s, and low, q D = 1 × 10 −4 m/s, injection rates, the progression of concentration field is different at the later time. It appears that the dynamic of the coupled processes controlled by the injection rate can impact the development of the dominant finger in the middle layer. Note that since the top and bottom layers, Ω 501 and Ω 502 , have lower permeability than the middle layer, Ω 500 , the flow mainly goes through the middle layer. Similar to the previous example, one of the two initial fingers becomes the main path connecting the inlet and the outlet boundaries. In Figure 5, the behavior of the concentration and velocity fields, together with temporal porosity alteration ( ∂φc ∂t ), are illustrated for both injection scenarios. As mentioned earlier, due to the difference of viscosity (µ) between that of the injected c and the in-situ c, two fingers developed at the beginning, see Figure 5a and e. For the high injection rate, the top finger, however, disappeared while the bottom finger progresses until it reaches the outlet ∂Ω 3 , see Figures 5b-d. One could see that the reaction front shown by ∂φc ∂t progresses similarly to the concentration front shown by the black contours. Besides, as the concentration field develops, the change in µ enhances the flow channeling illustrated by velocity arrows. For the low injection rate case shown in Figure 5e-h, the development of the concentration field is dissimilar to that of the high injection rate case as the top finger becomes a preferable path instead of the bottom one. Note that the dissolution and precipitation are shown in Figure 5 are a combined effect of injecting water that is unsaturated with respect to calcite and fluid pressure changes. It is clear that the majority of the dissolution occurs due to the transport of the injected water in the porous domain. For the animated version of Figure  5, please refer to Videos 1 and 2. These videos represent the flow and concentration field as well as ∂φc ∂t and illustrate the applicability of the presented coupled model for heterogeneous porous media.
Importantly, this example has illustrated the capability of our proposed method-which is equipped with the EG method-for handling discontinuous material properties across different layers and the sharp interface of the concentration species. Moreover, we have again observed the expected physical and chemical phenomena, including solid deformation, viscous fingering, and dissolution/precipitation.

Example 3
In the given computational domain Ω 500 = [0, 100] × [0, 30], we investigate the setup with the heterogeneous k values as shown in Figure 1c. A random field generator [103] is utilized to generate a heterogeneous permeability field with a given mean permeability of k = 1 × 10 −10 I m 2 , variance of 0.5, and correlation lengths in x-and y-direction of 5 and 1 m, respectively. The heterogeneous permeability field varies in two orders of magnitude. All other physical parameters are the same as in the previous examples.
Here, we focus on the interplay between the heterogeneous permeability and the HMC coupled processes. Similar to the previous examples, two different injection rates are applied. In Figure 6   Note that the magnitude of the mechanical deformation is higher than that of the calcite dissolution/precipitation and similar to what was observed in example 1. Therefore changes in porosity due to chemical reaction have a second-order effect on permeability compared to that of induced by the mechanical deformation. Videos 3 and 4 representing the flow and concentration field as well as ∂φc ∂t illustrate the applicability of the presented coupled model for heterogeneous porous media. n · n| e dS, (58) and the discrete numerical flux approximated by BDM,q n · n| e , is defined bȳ q n · n| e := −q n h · n ∀e ∈ E D,m h .
In Figure 8, the values of r n mass are illustrated for each case and time. One could see that the magnitude of r n mass is always less than 1 × 10 −5 , which is the tolerance set for the fixed-stress loop, see Algorithm 1; therefore, the framework is locally mass conservative. We note that the high injection rate case tends to the higher value of the magnitude of r n mass than that of the low injection rate case.

Example 4
Lastly, we investigate the performance of the proposed framework when the permeability field is anisotropic, and the grid is unstructured, as shown in Figure 1d. In the computational domain Ω 500 = [0, 100] × [0, 30], see Figure 1a, we consider the anisotropic permeability field to emphasize the capability of our proposed algorithm. The permeability tensor of this example is defined as follows: where k xx = 8.8 × 10 −10 m 2 and all other parameters are similar to all other cases. Figure 9 shows the reactive flow dynamics and the residual of mass. We observe that the flow in the horizontal direction dominates the flow in the vertical direction since the permeability in the horizontal direction is ten times higher than that of the vertical direction. Figure 9d-f illustrate that the proposed framework is locally mass conservative as the residual of mass values are always less than 1 × 10 −5 , which is the tolerance set for the fixed-stress loop.

Discussion
The main observations of the foregoing numerical examples can be summarized as follows: 1. The injection rate supplied at the inlet boundary is critical in defining flow behavior. The preferential flow paths developed through time are significantly different with different injection rates. Besides, the injection flow rate also controls the development of the advection and reaction fronts. 2. Using the applied set of the input parameters resulted in a more noticeable mechanical effect on the change in φ (and subsequently in k) compared to that of the calcite dissolution/precipitation effect. We note that this observation could vary with different sets of input parameters and required to be further investigated. The change in µ resulted from the change in c is significant, resulting in the development of preferential flow paths. 3. The results of both homogeneous and heterogeneous as well as isotropic and anisotropic permeability field show that our framework preserves mass locally. This property is essential for the coupled HMC system.
In terms of computational efficiency, it is noted that the iteration number for the fixedstress iteration was around three (four for the example 3) at the initial time stage, but it only required two iterations for the rest of the time for all the presented examples. For all examples, we have 31934, 23818, 7852, 11910 degrees of freedom for the displacement, flux, pressure, and concentration fields, respectively. The computational time was around 4.78×10 −5 second per degrees of freedom per each time step. All simulations were computed on XeonE5 2650v4 with a single thread.

Conclusion
This paper has presented a mixed finite element framework for coupled hydro-mechanicalchemical processes in heterogeneous porous media. The main advantage of the proposed framework is its relatively affordable cost to attain local conservation regardless of material anisotropy, thanks particularly to the use of the EG method. Through several numerical examples, we have demonstrated the performance and capabilities of the proposed framework with a focus on local conservation. The numerical results have highlighted how the overall behavior is influenced by different processes, including solid deformation, calcite dissolution, and fluid viscosity alteration. The developed numerical model can provide insight into how the interactions among HMC processes and heterogeneity manifest themselves at a larger scale. Future work includes an extension of the modeling framework to coupled thermohydro-mechanical-chemical processes in heterogeneous and/or fractured porous media.

Acknowledgements
This research has received financial support from the Danish Hydrocarbon Research and Technology Centre under the Advanced Water Flooding program. The computational results in this work have been produced by the multiphenics library [101], which is an extension of FEniCS [100] for multiphysics problems. We acknowledge the developers of and contributors to these libraries.