A Locally Conservative Mixed Finite Element Framework for Coupled Hydro-Mechanical-Chemical Processes in Heterogeneous Porous Media
AA Locally Conservative Mixed Finite Element Framework forCoupled Hydro-Mechanical-Chemical Processes in HeterogeneousPorous Media
T. Kadeethum a,b, ∗ , S. Lee c , F. Ballarin d , J. Choo e , H.M. Nick a a Technical University of Denmark, Denmark b Cornell University, New York, USA c Florida State University, Florida, USA d mathLab, Mathematics Area, SISSA, Italy e The University of Hong Kong, Hong Kong
Abstract
This paper presents a mixed finite element framework for coupled hydro-mechanical-chemicalprocesses in heterogeneous porous media. The framework combines two types of locally con-servative 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. Thiscombination ensures local mass conservation, which is critical to flow and transport in het-erogeneous porous media, with a relatively affordable computational cost. A particular classof the framework is constructed for calcite precipitation/dissolution reactions, incorporatingtheir nonlinear effects on the fluid viscosity and solid deformation. Linearization schemesand algorithms for solving the nonlinear algebraic system are also presented. Through nu-merical examples of various complexity, we demonstrate that the proposed framework is arobust and efficient computational method for simulation of reactive flow and transport indeformable porous media, even when the material properties are strongly heterogeneous andanisotropic.
Keywords: hydro-mechanical-chemical coupling, poroelasticity, reactive flow, mixed finiteelement method, enriched Galerkin method, local conservation
1. Introduction
Hydro-mechanical-chemical (HMC) processes in porous media, in which fluid flow, soliddeformation, and chemical reactions are tightly coupled, appear in a variety of problems ∗ corresponding author Email addresses: [email protected] (T. Kadeethum), [email protected] (S. Lee), [email protected] (F. Ballarin), [email protected] (J. Choo), [email protected] (H.M. Nick) Author contribution statement is listed in Section 7. This research has been awarded the funding from the 2019 Computers & Geosciences Research grant.
Preprint submitted to Elsevier October 13, 2020 a r X i v : . [ c s . C E ] O c t anging from groundwater and contaminant hydrology to subsurface energy production [1–7]. The multiphysical interactions in these problems give rise to strong heterogeneity inthe material properties. For instance, change in pore pressure perturbs effective stress inthe solid matrix, which can, in turn, alter the conductivity and storability of the porousmedium [8–13]. Similarly, chemical processes can result in the precipitation or dissolutionof solid minerals, which decreases or increases the pore volume, respectively, and thus, theconductivity [3, 4, 6, 14–16]. Therefore, accurate numerical modeling of coupled HMCproblems requires a computational method that can robustly handle strong heterogeneityin porous media.Numerical simulation of multiphysical problems in porous media has been a subject ofextensive research (e.g. [17–27]), and lots of software packages have been developed forthis purpose. Notable examples include: (1) TOUGH software suite, which includes multi-dimensional numerical models for simulating the coupled thermo-hydro-mechanical-chemical(THMC) processes in porous and fractured media [15, 28–30], (2) SIERRA Mechanics, whichhas simulation capabilities for coupling thermal, fluid, aerodynamics, solid mechanics andstructural dynamics [31], (3) PyLith, a finite-element code for modeling dynamic and quasi-static simulations of coupled multiphysics processes [32], (4) OpenGeoSys project, whichis developed mainly based on the finite element method using object-oriented program-ming THMC processes in porous media [33], (5) IC-FERST, a reservoir simulator basedon control-volume finite element methods and dynamic unstructured mesh optimization[34–37], (6) DYNAFLOW ™ , a nonlinear transient finite element analysis platform [38], (7)DARSim, multiscale multiphysics finite volume based simulator [39–41], (8) the CSMP, anobject-oriented application program interface, for the simulation of complex geological pro-cesses, e.g. THMC, and their interactions [42, 43], and (9) PorePy, an open-source modelingplatform for multiphysics processes in fractured porous media [44].Nevertheless, it remains challenging to simulate coupled HMC processes in porous me-dia in a robust and efficient manner, especially when the material properties are highlyheterogeneous and/or anisotropic. Because HMC problems involve transport phenomenain heterogeneous porous media, the numerical method for these problems must ensure lo-cal (element-wise) conservation [45, 46]. The most practical method featuring local massconservation may be the finite volume method with a standard two-point flux approxima-tion scheme. However, this standard finite volume method requires the grid to be alignedwith the principal directions of the permeability/diffusivity tensors [6, 47], which inhibitsthe use of an unstructured grid when the permeability/diffusivity tensors are anisotropic.Multi-point flux-approximation methods have been developed to tackle this issue, but theirimplementation is often complicated and onerous [48]. Discontinuous Galerkin (DG) meth-ods offer an elegant way to handle arbitrarily anisotropic tensor conductivity/diffusivity.However, their computational cost is often impractical as a result of the proliferation of thedegrees of freedom.In this paper, we present a new framework for computational modeling of coupled HMCprocesses in porous media, which efficiently provides local mass conservation even when thematerial properties are strongly heterogeneous and anisotropic. The proposed frameworkcombines two types of discretization methods: (1) an enriched Galerkin (EG) method for2eactive flow and transport, and (2) a three-field mixed finite element method for coupledhydro-mechanical processes. The EG method, which has recently been developed and ad-vanced in the literature [46, 48–52], augments a piecewise constant function to the continuousGalerkin (CG) function space. This method uses the same interior penalty type form as theDG method, but it requires a substantially fewer number of degrees of freedom than theDG method. Thus the EG method can provide locally conservative solutions to the reactiveflow system regardless of the grid–conductivity alignment. For the hydro-mechanical sub-system of the HMC problem, we use a three-field mixed finite element formulation [53–55],which provides locally conservative, high-order solutions to the fluid velocity field. Specifi-cally, we employ the Lagrange finite elements for approximating the displacement field, theBrezzi-Douglas-Marini (BDM) element for the fluid velocity field, and the piecewise con-stant element for the fluid pressure field. It is noted that this combination of elements is ourpersonal choice, and one may use another combination for the same three primary fields asin [55–57].The purpose of this work is to develop an accurate numerical method for tackling coupledHMC processes in heterogeneous porous media, with a practically affordable computationalcost. Our specific objectives can be summarized as follows:1. To formulate a robust numerical approximation scheme for coupled HMC processesin heterogeneous porous media, employing a combination of locally conservative finiteelement methods.2. To reduce the computational cost for solving an advection-diffusion-reaction equationby using the EG method, which requires approximately two and three times fewerdegrees 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 mod-eling tightly coupled HMC problems with homogeneous to heterogeneous, isotropic toanisotropic permeability fields with local conservation.The rest of the paper is organized as follows. Section 2 describes the governing equationsof coupled HMC processes. Section 3 explains the discretization methods, linearizationtechniques, and solution algorithms of the proposed framework. Section 4 presents severalnumerical examples of various complexity and discusses key points found in this paper.Section 5 concludes the work.
2. Governing equations
This section briefly describes all the equations used in this study, namely poroelasticityand advection-diffusion-reaction equations.Let Ω ⊂ R d ( d ∈ { , , } ) denote the computational domain and ∂ Ω denote the bound-ary. The time domain is denoted by T = (0 , T] with T >
0. Primary variables usedin this paper are q ( · , t ) : Ω × T → R d , which is a vector-valued Darcy velocity (m / s), p ( · , t ) : Ω × T → R , which is a scalar-valued fluid pressure (Pa), u ( · , t ) : Ω × T → R d , whichis a vector-valued displacement (m), c i : Ω × T → R , which is the i -th component of chemicalconcentration (mmol / m ). 3 .1. Poroelasticity To begin, we adopt Biot’s poroelasticity theory for coupled hydro-mechanical processesin porous media [59, 60]. Although poroelasticity may oversimplify deformations in softporous materials such as soils [61–64], it would be reasonably good for stiff materials suchas rocks, which is the focus of this work. The poroelasticity theory provides two coupledgoverning equations, namely linear momentum and mass balance equations. Under quasi-static conditions, the linear momentum balance equation can be written as ∇ · σ ( u , p ) + f = , (1)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. Thegravitational force will be neglected in this study, but the body force term will be kept inthe succeeding formulations for a more general case. Further, σ is the total stress tensor,which may be related to the effective stress tensor σ (cid:48) and the pore pressure p as σ ( u , p ) = σ (cid:48) ( u ) − αp I . (2)Here, I is the second-order identity tensor, and α is the Biot coefficient defined as [65]: α = 1 − KK s , (3)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 withthe displacement vector, which can be written as σ (cid:48) ( u ) = λ l tr( ε ( u )) I + 2 µ l ε ( u ) . (4)Here, ε is the infinitesimal strain tensor, defined as ε ( u ) := 12 ( ∇ u + ( ∇ u ) (cid:124) ) , (5)and λ l and µ l are the Lam´e constants, which are related to the bulk modulus and the Poissonratio ν of the solid matrix as λ l = 3 Kν ν , and µ l = 3 K (1 − ν )2(1 + ν ) . (6)For this solid deformation problem, the domain boundary ∂ Ω is assumed to be suitablydecomposed into displacement and traction boundaries, ∂ Ω u and ∂ Ω t , respectively. Then thelinear momentum balance equation is supplemented by the boundary and initial conditionsas: ∇ · σ (cid:48) ( u ) + α ∇ · ( p I ) + f = in Ω × T , u = u D on ∂ Ω u × T , σ ( u ) · n = t D on ∂ Ω t × T , u = u in Ω at t = 0 , (7)4here u D and t D are prescribed displacement and traction values at the boundaries, respec-tively, and n is the unit normal vector to the boundary.Next, the mass balance equation is given as [4, 19, 43, 66]:1 M ∂p∂t + α ∂ε v ∂t + ∂φ c ∂t + ∇ · q = g in Ω × T , (8)where 1 M = (cid:18) φ c f + α − φ K s (cid:19) (9)is the Biot modulus. Here, c f is the fluid compressibility, φ is the initial porosity, ε v := tr( ε ) = ∇ · u is the volumetric strain, and g is a sink/source term. Because we willintroduce 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 q = − k ( φ ) µ ( c i ) ( ∇ p − ρ g ) . (10)Note that here the fluid viscosity µ is considered a function of concentration c i . Again, thegravitational force, ρ g , will be neglected in this work, without loss of generality. In addition, k ( φ ) is the matrix permeability tensor defined as k := k mult ( φ ) k xx k xy k xz k yx k yy k yz k zx k zy k zz if d = 3 ,k mult ( φ ) (cid:34) k xx k xy k yx k yy (cid:35) if d = 2 ,k mult ( φ ) k if d = 1 , (11)The k xx , k yy , and k zz represent the matrix permeability in x -, y -, and z -direction, respec-tively. The k mult ( φ ) is a multiplier used to update k when φ is altered, which will bedescribed later.For the fluid flow problem, the domain boundary ∂ Ω is also suitably decomposed intothe pressure and flux boundaries, ∂ Ω p and ∂ Ω q , respectively. In what follows, we apply thefixed stress split scheme [19, 68], assuming ( σ v − σ v, ) + α ( p − p ) = Kε v . Then we write5he fluid flow problem with boundary and initial conditions as (cid:18) M + α K (cid:19) ∂p∂t + αK ∂σ v ∂t + ∂φ c ∂t + ∇ · q = g in Ω × T ,p = p D on ∂ Ω p × T , q · n = q D on ∂ Ω q × T ,p = p in Ω at t = 0 , (12)where σ v := tr( σ ) is the volumetric stress, and p D and q D are the given boundary pressureand flux, respectively. An advection-diffusion-reaction system for N c number of the miscible species is given bythe following equations. For all i = 1 , . . . , N c , ∂∂t ( φc i ) + ∇ · η ( q , c i ) = q i ( c i ) , in Ω × T , (13)where q i ( c i ) is a reaction term coupled with sink/source for each component, and the massflux η ( q , c i ) is defined as η ( q , c i ) := q c i − D e,i ( φ ) ∇ c i . (14)Here D e,i ( φ ) is the effective diffusion coefficient tensor defined as D e,i := φτ D i , (15)where τ = φ − [69, 70] and D i is the given diffusion coefficient tensor. The boundary forthe 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 < } and ∂ Ω out := { x ∈ ∂ Ω : q · n ≥ } . (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 requireseight transport equations to solve the concentration values of the following main species/ions: { H + , Ca , CaHCO +3 , OH − , CO − , HCO − , H CO ∗ , CaCO ∗ (Aq) } [3, 67, 71, 72]. For sim-plicity, in this paper we consider a reduced system based on the empirical relationshippresented in [3, 4, 67], in which N c decreases to 1. Thus, letting c := c , we write theadvection-diffusion-reaction system with its boundary and initial conditions as follows: ∂∂t ( φc ) + ∇ · ( q c ) − ∇ · ( D e ( φ ) ∇ c ) = q in Ω × (0 , T ] ,η ( q , c ) · n = c in q · n on ∂ Ω in × (0 , T ] , D e ( φ ) ∇ c · n = 0 on ∂ Ω out × (0 , T ] ,c = c in Ω at t = 0 , (17)6here c in is the inflow concentration, c is the initial concentration, and q represents a sourceterm reflecting the calcite dissolution/precipitation reactions. For this term, here we adoptthe term in [3, 4, 67], given by q = R c A s , (18)where A s is the specific surface of the porous medium, and R c is the reaction rate calculatedas R c = r , for (cid:101) c > , − r , for (cid:101) c < , . , for (cid:101) c = 0 , (19)with (cid:101) c = ( c eq − c ) c eq , (20) r = a + a τ + a log | (cid:101) c | + a τ + a τ log | (cid:101) c | + a (log | (cid:101) c | ) , (21)and c eq = 1 . × − + 3 . × − p − . × − τ − . × − p + 4 . × − pτ − . × − τ . (22)Here, τ is the medium temperature, and a , a , a , a , a , and a are defined in Table 1. Table 1: Coefficients of the (21) for different range of (20) a a a a a a (cid:101) c > .
01 -5.73 1 . × − . × − − . × − . × − − . < (cid:101) c -6.45 2 . × − − . × − . × − . × − − . × − − . < (cid:101) c ≤ .
01 -5.80 1 . × − . × − . × − . × − − . × − Before closing this section, we describe physical properties that are coupled with primaryvariables, u , q , p , and c . The porosity change due to solid deformation may be expressedas [18, 59, 73]: φ m = φ + ( α − φ ) ( (cid:15) v − (cid:15) v ) + ( α − φ ) (1 − α ) K ( p − p ) , (23)where (cid:15) v is the initial volumetric strain. The porosity alteration due to calcite dissolution/-precipitation is calculated as Γ ( u , c ) = ∂φ c ∂t = R c A s ρ s ω , (24)where ω is the number of moles of total precipitated species per kilogram of rock (assumed tobe 10.0 in this study following [3, 4]), and ρ s = 2500 kg / m is used for throughout this paper.7ote that this term, (24), enters (12). Also, the terms φ m and φ c are used to distinguishbetween the changes in φ due to solid deformation as in (23), and chemical reactions asin (24), respectively. The changes in porosity due to (23) and (24), also affect the specificsurface for porous medium ( A s ) as A s = A φφ log( φ )log ( φ ) , (25)where A is the initial value of A s , and it is set as 5000 throughout this study [74]. Further-more, the porosity change influences the matrix permeability as [75–77]: k = k k mult ( φ ) = k exp (cid:18) b (cid:18) φφ − (cid:19)(cid:19) , (26)where k is the initial matrix permeability and b is an empirical parameter determinedexperimentally. In this work, we set b = 22 . c also affects µ ,and we adopt the specific form from [78–80], given by µ = log ( µ l ) + (cid:18) c − c l c h − c l (cid:19) (log ( µ h ) − log ( µ l )) , (27)where c l and c h are lower and higher bounds of the concentration, and µ l and µ h are fluidviscosity corresponding to c l and c h , respectively. Table 2 summarizes the effects of physicalprocesses 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 absenceof a relationship. Table 2: Summary of the effects of individual physical processes on physical properties
Physical properties Mechanical deformation Fluid pressure Calcite concentration φ (23) (23) (24) k (23) + (26) (23) + (26) (24) + (26) µ - - (27) D e (23) + (15) (23) + (15) (24) + (15) A s (23) + (25) (23) + (25) (24) + (25)
3. Numerical methods
In this section, we describe the numerical methods for the governing system described inthe previous sections. Here, we utilize a combination of a mixed finite element method forspatial discretization, and employ both a backward differentiation formula and an explicitRunge-Kutta method for temporal discretization.8 .1. Domain discretization and geometrical quantities
We begin by introducing the notations used throughout this paper. Let T h be a shape-regular 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 − T ∈ T h ) and by E Ih and E ∂h the collection ofall interior and boundary facets, respectively. The boundary set E ∂h is decomposed into twodisjoint subsets associated with the Dirichlet boundary faces, and the Neumann boundaryfaces for each of (7) and (12). In particular, E D,uh and E N,uh correspond to the faces on ∂ Ω u and ∂ Ω tr , respectively, for (7). On the other hand, for (12), E D,mh and E N,mh 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 Ih , where T + and T − are the two neighboring elements to e . We denote by h e the characteristiclength of e calculated as h e := meas ( T + ) + meas ( T − )2 meas( e ) , (28)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. Forany 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 weightedaverage operator as { ζ } δe = δ e ζ + + (1 − δ e ) ζ − , on e ∈ E Ih , (29)and { τ } δe = δ e τ + + (1 − δ e ) τ − , on e ∈ E Ih , (30)where δ e is calculated by [81, 82]: δ e := k − e k + e + k − e . (31)Here, k + e := (cid:0) n + (cid:1) (cid:124) · k + n + , and k − e := (cid:0) n − (cid:1) (cid:124) · k − n − , (32)where k e is a harmonic average of k + e and k − e which reads k e := 2 k + e k − e k + e + k − e , (33)and k is defined as in (11). The jump across an interior edge will be defined as[[ ζ ]] = ζ + n + + ζ − n − and [[ τ ]] = τ + · n + + τ − · n − on e ∈ E Ih . e ∈ E ∂h , we set { ζ } δ e := ζ and { τ } δ e := τ for what concerns the definitionof the weighted average operator, and [[ ζ ]] := ζ n and [[ τ ]] := τ · n as definition of the jumpoperator. The time domain T = (0 , T] is partitioned into N subintervals such that 0 =: t < t < · · · < t N := T. The length of each subinterval ∆ t n − is defined as ∆ t n − = t n − t n − where n represents the current time step. We assume that the user provides the initial ∆ t , whilean adaptive procedure is carried out to choose ∆ t n − , n >
1, as follows:∆ t n − := (cid:40) CFL h l (cid:107) q n − (cid:107) ∞ if ∆ t n ≤ ∆ t max ∆ t max if ∆ t n > ∆ t max , (34)where CFL is a constant that the user can provide according to the Courant-Friedrichs-Lewycondition [83], (cid:107)·(cid:107) ∞ is the maximum norm of a vector function, and ∆ t max is a maximumallowed time step. Note that we use ∆ t max as a tool to control ∆ t n as the model approachesa steady-state condition since (cid:107) q n − (cid:107) ∞ may approach zero, which would lead to a very largeratio h l (cid:107) q n − (cid:107) ∞ .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–86]BDF m ( ϕ n ) := t n ( ϕ n − ϕ n − ) m = 1 t n (3 ϕ n − ϕ n − + ϕ n − ) m = 2 t n (11 ϕ n − ϕ n − + 9 ϕ n − − ϕ n − ) m = 3 t n (25 ϕ n − ϕ n − + 36 ϕ n − − ϕ n − + 3 ϕ n − ) m = 4 (35)for the discretization of the time derivative of ϕ ( · , t ) at time t n . We also utilize the explicitRunge-Kutta methods [24, 87]: RK ( ϕ n ) = ϕ n +1 = ϕ n + κ , (36) κ = ∆ t n F ( X n , Y n ) , for the first order Runge-Kutta method corresponding to the explicit Euler method, andRK ( ϕ n ) = ϕ n +1 = ϕ n + 16 κ + 13 κ + 13 κ + 16 κ , (37) κ = ∆ t n F ( X n , Y n ) ,κ = ∆ t n F (cid:0) X n + ∆ t n , Y n + κ (cid:1) ,κ = ∆ t n F (cid:0) X n + ∆ t n , Y n + κ (cid:1) ,κ = ∆ t n F ( X n + ∆ t n , Y n + κ ) , F ( X n , Y n ) is any functions with independentvariable X and dependent variable Y [8, 87], which we will specify in the linearization andsolving processes in Section 3.4.Finally, we define an extrapolation operator as follows [8, 24]:EX ( ϕ ) = ˆ ϕ n +1 = (cid:18) t n ∆ t n − (cid:19) ϕ n − ∆ t n ∆ t n − ϕ n − if n ≥ ,ϕ if n = 0 , (38)and in the following we will adopt the notation ˆ ϕ n +1 to denote an extrapolation value of { ϕ n , ϕ n − } . In this framework, the displacement field is approximated by the classical continuousGalerkin method (CG) method, and the fluid velocity and pressure fields are discretizedby the Brezzi-Douglas-Marini (BDM) element [88] and the piecewise constants discontin-uous Galerkin (DG) method, respectively, to ensure local mass conservation. Lastly, theconcentration 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-valuedfunction: U CG k h ( T h ) := (cid:8) ψ u ∈ C (Ω; R d ) : ψ u | T ∈ P k ( T ; R d ) , ∀ T ∈ T h (cid:9) , (39)where C (Ω; 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 genericfunction of U CG k h ( T h ). In addition, the CG space for scalar-valued functions is defined as: P CG k h ( T h ) := (cid:8) ψ p ∈ C (Ω) : ψ p | T ∈ P k ( T ) , ∀ T ∈ T h (cid:9) , (40)where C (Ω) := C (Ω; R ) and P k ( T ) := P k ( T ; R ). Next, we define the following DG functionspace: P DG k h ( T h ) := (cid:8) ψ p ∈ L (Ω) : ψ p | T ∈ P k ( T ) , ∀ T ∈ T h (cid:9) , (41)where L (Ω) is the space of square-integrable scalar functions. This non-conforming finiteelement space allows us to consider discontinuous functions and coefficients rigorously. Wethen define the EG finite element space with polynomial order k as: P EG k h ( T h ) := P CG k h ( T h ) ⊕ P DG h ( T h ) , (42)i.e., a CG finite element space enriched by the space P DG h ( T h ) of piecewise constant func-tions. 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]: V BDM k h ( T h ) := { ψ v ∈ H (div , Ω) : ψ v | T ∈ BDM( T ) , ∀ T ∈ T h } (43)where ψ v denotes a generic function of V BDM k h ( T h ) and BDM( T ) is defined according to [88].11 .3.1. Fully discrete form We now present the fully discrete form of the coupled HMC problem using the above-described combination of finite element spaces. In particular, we seek the approximateddisplacement solution u h ∈ U CG h ( T h ) as done in [20, 51, 89], fluid pressure p h ∈ P DG h ( T h ),velocity approximation q h ∈ V BDM h ( T h ), and concentration approximation c h ∈ P EG h ( T h ).We multiply the linear momentum balance equation (7) by a test function ψ u ∈ U CG h ( T h ).The fully discretized linear momentum balance equation thus has the following form: N u ( ψ u ; u nh , p nh ) = 0 , ∀ ψ u ∈ U CG h ( T h ) , (44)at each time step t n , where N u ( ψ u ; u nh , p nh ) = (cid:88) T ∈T h (cid:90) T σ (cid:48) ( u nh ) : ∇ s ψ u dV − (cid:88) T ∈T h (cid:90) T αp nh I : ∇ s ψ u dV − (cid:88) T ∈T h (cid:90) T f ψ u dV − (cid:88) e ∈E N,uh (cid:90) e t D ψ u dS, ∀ ψ u ∈ U CG h ( T h )Here (cid:82) T · dV and (cid:82) e · dS refer to volume and surface integrals, respectively, and ∇ s isthe symmetric gradient operator. Furthermore, the notation for N u ( ψ u ; u nh , p nh ) in (44)highlights before the semicolon the test function, and after the semicolon the (possiblynonlinear) dependence on discrete solutions to the coupled problem. The same notation willbe used hereafter for the remaining equations.Next, the weak form of the mass balance equation (12) is obtained multiplying by ψ p ∈P DG h ( T h ) and integrating by parts, resulting in: N p ( ψ p ; p nh , q nh , c nh ) = 0 , ∀ ψ p ∈ P DG h ( T h ) , (45)for each time step t n , where N p ( ψ p ; p nh , q nh , c nh ) = (cid:88) T ∈T h (cid:90) T (cid:18) M + α K (cid:19) BDF ( p nh ) ψ p dV + (cid:88) T ∈T h (cid:90) T ∇ · ( q nh ) ψ p dV + (cid:88) T ∈T h (cid:90) T αK RK ( σ v ) ψ p dV + (cid:88) T ∈T h (cid:90) T RK ( φ c ) ψ p dV − (cid:88) T ∈T h (cid:90) T gψ p dV. For the Darcy velocity equation (10), we obtain N v ( ψ v ; u nh , p nh , q nh , c nh ) = 0 , ∀ ψ v ∈ V BDM h ( T h ) . (46)where 12 v ( ψ v ; u nh , p nh , q nh , c nh ) := (cid:88) T ∈T h (cid:90) T p nh ∇ · ψ v dV + (cid:88) T ∈T h (cid:90) T k ( u nh , c nh ) − µ ( c nh ) q nh ψ v dV + (cid:88) e ∈E D,mh (cid:90) e p D ψ v · n dS. Lastly, for the advection-diffusion-reaction equations of species transport we write: N c ( ψ c ; u nh , q nh , c nh ) = 0 , ∀ ψ c ∈ P EG h ( T h ) (47)for each time step t n , where N c ( ψ c ; u nh , q nh , c nh ) = (cid:88) T ∈T h (cid:90) T φ BDF ( c nh ) ψ c dV + (cid:88) T ∈T h (cid:90) T D e ∗ ( φ n ) ∇ c nh · ∇ ψ c dV − (cid:88) e ∈E Ih (cid:90) e { D e ∗ ( φ n ) ∇ c nh } δ e · (cid:74) ψ c (cid:75) dS + θ (cid:88) e ∈E Ih (cid:90) e { D e ∗ ( φ n ) ∇ ψ c } δ e · (cid:74) c nh (cid:75) dS + (cid:88) e ∈E Ih (cid:90) e βh e D e ∗ ( φ n ) e (cid:74) c nh (cid:75) · (cid:74) ψ c (cid:75) dS − (cid:88) T ∈T h (cid:90) T q nh c nh · ∇ ψ c dV + (cid:88) e ∈E Ih (cid:90) e q nh · n c up h (cid:74) ψ c (cid:75) dS + (cid:88) e ∈E Out h (cid:90) e q nh · n c nh ψ c dS − (cid:88) T ∈T h (cid:90) T R c A s ψ c dV + (cid:88) e ∈E In h (cid:90) e q nh · n c in ψ c dS. We note that the D e ∗ ( φ n ) is redefining D e ( φ n ) by including the numerical stabilizationterm, where D e ∗ ( φ n ) := D e ( φ n ) + γh (cid:107) q nh (cid:107) I , (48)as defined in [90–92]. The γh (cid:107) q nh (cid:107) I term is often referred as the first order artificial diffu-sivity coefficient [93, 94]. In our paper, we set the tuning parameter γ = 0 .
25. Alternativestabilization strategies including streamline diffusion and crosswind diffusion, or entropyviscosity methods could be also utilized to reduce oscillations in the numerical solution tothe concentration field [90–92, 94–98]. 13lso, c up h is an upwind value of c nh defined as [45, 99]: c up h = (cid:26) c n + h if q nh · n ≥ c n − h if q nh · n < ∀ e = ∂T + ∩ ∂T − (49)where c n + h and c n − h correspond to c nh of T + and T − , respectively.Lastly, the two parameters θ and β define corresponding interior penalty methods. Thediscretization becomes the symmetric interior penalty Galerkin method (SIPG) when θ = −
1, the incomplete interior penalty Galerkin method (IIPG) when θ = 0, and the non-symmetric interior penalty Galerkin method (NIPG) when θ = 1 [45]. In this study, we set θ = − β = 1 . Remark 1.
For the momentum balance equation (7), the traction boundary condition t D (traction) is applied weakly on each e ∈ E N,uh in (44), while the displacement boundarycondition u D is strongly enforced on each e ∈ E D,uh . For the mass balance equation (12),since we use a mixed formulation, the flux boundary condition q D is strongly applied oneach e ∈ E N,mh , but the pressure boundary condition p D is weakly applied on each E D,mh 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 dis-cretization schemes for each equation. We use BDF for the time discretization of the massbalance 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 to capture asharp front in the advection dominated regime [45]. 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 monolithicmanner, here we focus on developing a splitting algorithm for sequential solution to the cou-pled system, which can provide more flexibility especially when different software packagesneed to be combined. The overall computational strategy is summarized in Algorithm 1.14 lgorithm 1
Splitting algorithm for hydro-mechanical-chemical coupling model Initialize all input parameters (cid:46) p and c must be provided. Solve the equilibrium state for u (cid:46) see (44) Update φ , k , D e , and A s (cid:46) see (23), (26), (15), (25) for each time step t n do Part 1: coupling solid and fluid mechanics Set ι → p n − h → p n,ι =0 h , q n − h → q n,ι =0 h , u n − h → u n,ι =0 h for each fixed stress iteration step ( · ) n,ι until δφ n,ι < TOL do Solve (45) and (46) w.r.t. p nh and q nh for fixed u nh := u n,ι − h , c nh := ˆ c nh to get p n,ιh , q n,ιh Calculate φ n,ιf (cid:46) see (51) Solve (44) w.r.t. u nh and for fixed p nh := p n,ιh to get u n,ιh Calculate φ n,ιm (cid:46) see (23) Evaluate F ( ˙ σ vn,ι ) (cid:46) see (52) Evaluate δφ n,ι (cid:46) see (50) Update k n,ι (cid:46) see (23), (26) end for p n,ιh → p nh , q n,ιh → q nh , u n,ιh → u nh Part 2: chemical process
Update φ n , D ne , and A ns (cid:46) see (23), (15), (25) Solve (47) w.r.t. c nh for fixed u nh := u nh , q nh := q nh to get c nh Extrapolate ˆ c n +1 h (cid:46) see (53) Calculate ∆ t n +1 (cid:46) see (34) Evaluate F (ˆ q n +1 ) and F (cid:16) ˆ˙ φ n +1 c (cid:17) (cid:46) see (55), (54) Update ˆ φ n +1 , ˆ k n +1 , ˆ D n +1 t , ˆ µ n +1 , and ˆ A sn +1 (cid:46) see (53), (24), (26), (15), (27), (25) p nh → p n − h , q nh → q n − h , u nh → u n − h , c nh → c n − h (cid:46) update time step n − Output end for
In Algorithm 1, we separate our algorithm into two parts. The first part (lines 8 to17) focuses on solving the coupled hydro-mechanical problem, (44), (45), and (46), usingthe 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 usinga monolithic method (line 9) based on given displacement u n,ι − h from previous nonlineariteration and concentration extrapolated ˆ c nh 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 currentnonlinear iteration (line 11). The convergence criterion is based on δφ n,ι (Algorithm 1 line8), which is defined as: δφ n,ι := φ n,ιm − φ n,ιf φ n,ιm . (50)15ere, φ n,ιm is the porosity resulting from the solid deformation (23) and φ n,ιf is the porosityresulting from the fluid flow problem defined as [18, 68, 73]: φ n,ιf = φ n − + ( α − φ n − ) K (cid:0) p n,ι − p n − (cid:1) , (51)where ( · ) ι represents iteration counter inside the fixed-stress loop. From the fixed stress splitconcept (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 toleranceTOL is set as 1 × − throughout this study. Note that the flow equations, (45) and (46),are solved by assuming that ∂σ v ∂t = 0, i.e., F ( ˙ σ vn,ι ) is frozen; therefore, this term is evaluatedexplicitly after the momentum equation (44) is solved, as illustrated in Algorithm 1 line 13[19, 68], and F ( ˙ σ vn,ι ) is defined as:F ( ˙ σ vn,ι ) := (cid:88) T ∈T h (cid:90) T αK RK (cid:0) σ v ( u n,ι , u n − , p n,ι , p n − ) (cid:1) ψ p dV. (52)The second part (from line 18) focuses on solving advection-diffusion-reaction equation(47), using q nh , φ n , D ne , and A ns obtained from the first part. One could view this strategyas a one-way coupling scheme between coupled hydro-mechanical and advection-diffusion-reaction equations. Next, Algorithm 1 line 21 linearizes c n +1 h by extrapolating c nh and c n − h to ˆ c n +1 h by using (38): ˆ c n +1 h = EX ( c h ) (53)where (ˆ · ) n represents an extrapolation value based on the extrapolation described in (38).Subsequently, we evaluate F (cid:16) ˙ φ cn,ι (cid:17) and F ( q n,ι ), which are defined asF (cid:16) ˆ˙ φ n +1 c (cid:17) := (cid:88) T ∈T h (cid:90) T RK (cid:0) φ c (ˆ c n +1 , c n ) (cid:1) ψ p dV, (54)and F (cid:0) ˆ q n +1 (cid:1) := (cid:88) T ∈T h (cid:90) T R c (ˆ c n +1 , p n ) A s (ˆ c n +1 , u n ) ψ c dV, (55)using ˆ c n +1 h calculated by (53). We note that the equation (47) becomes linear by employingˆ c n +1 h to calculate F ( q ). Also, the porosity alteration as a result of calcite dissolution/pre-cipitation (Algorithm 1 line 24) is computed byˆ φ n +1 = ˆ φ n +1 c = RK (cid:0) Γ (cid:0) u n , ˆ c n +1 h (cid:1)(cid:1) . (56)Note that the porosity change due to the calcite dissolution/precipitation reactions is ad-ditional to the porosity change by solid deformation, (23). Subsequently, ˆ k n +1 , ˆ D n +1 t , andˆ A sn +1 are determined using ˆ φ n +1 . Lastly, we also calculate ˆ µ n +1 using ˆ c n +1 h , see (53) and(27). 16or 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 areemployed from PETSc package [102]. All simulations are computed on XeonE5 2650v4 witha single thread. Remark 3.
We note that the EG method, which is used to approximate the advection-diffusion-reaction (17), is based on the Galerkin method, which could be extended to con-sider 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 dis-continuities are observed, can be developed.
4. Numerical examples
In this section, we demonstrate the performance and capabilities of the proposed nu-merical method through various numerical examples. We begin with a single-layer modelcomparing the performance for single-phase flow with chemical dissolution/precipitation andsolid deformation. Then we illustrate the performance of the developed model for a layeredmedium as well as a heterogeneous single-layer medium. Lastly, we test the proposed frame-work using an example with an anisotropic permeability field. All four examples and theirmesh are illustrated in Figure 1. More detailed setup, including the input parameters andthe boundary conditions of each example, are described in the beginning of each example.
In the first example, the computational domain is defined as Ω = [0 , × [0 , K = 8 . α = 0 . ν = 0 . φ = 0 . k = 8 . × − I m .The fluid properties considered in this case are c f = 1 . × − Pa − , ρ = 1000 kg / m , D = 1 . × − m / s, and µ is calculated using (27) by setting µ h = 5 . µ l = 1 . × − Pa / s corresponding to c h = 1 .
68 and c l = 0 . / m , respectively. Next, the boundaryconditions for all these examples are applied as follows. For the momentum balance equation(7), we assume u D · n = 0 on ∂ Ω , ∂ Ω , and ∂ Ω . Furthermore, t D = [0 . , − . × ] Pa isapplied on ∂ Ω . Therefore, the medium is under compression. For the mass balance equation(12), the boundary condition q D = 0 is set on ∂ Ω and ∂ Ω and we impose p D = 1 × Pa on ∂ Ω . Here, for the mass balance equation (12), we test two different scenarios on ∂ Ω , where scenario (a) corresponds to q D = 2 × − m / s and scenario (b) is characterizedby q D = 1 × − m / s. Thus, scenarios (a) and (b) will be referred to as the high andlow injection rate cases, respectively. Since we want to compare the results of the abovescenarios at the same total injected volume (I.V.), which is defined asI . V . = q D t n A d , (57)where A d = 30m is the surface area of ∂ Ω , 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 . ∂ Ω . The initial pressure p is 1 × Pa, the initial concentration c is17 m m ∂ Ω ∂Ω ∂Ω ∂ Ω Ω Ω Ω Ω ( k ) (b)(c)(d)(a) Figure 1: Geometry and notation used to define material properties; ( a ) example 1: single-layer porousmedium (Ω ), ( b ) example 2: three-layer porous medium (Ω and Ω ), ( c ) example 3: heterogeneousporous medium (the arithmetic mean of k = 8 . × − I m with correlation length in x - and y -directionof 5 and 1 m, respectively.), and ( d ) mesh used for all examples (number of element is 7852). calculated by (22) using p = p and τ = 20 C, and the initial displacement u is calculatedas 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 de-veloped HMC coupled numerical scheme in a homogeneous porous medium for two differentinjection rates. The aim is to illustrate the impact of different processes on the advanceof the flow path and reactive solute transport. Initially, the composition of the pore fluidwithin the porous medium is in equilibrium with calcite. Note that c eq calculated by (22) isa 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 pres-sure field as a result of fluid injection on the left boundary and the fluid production on theright boundary varies the c eq resulting in precipitation or dissolution in the domain. Theinjected water is also unsaturated with respect to calcite. Therefore, the injected fluid, asadvances into the domain, will dissolve the calcite mineral.Figure 2 shows the concentration fields at different injected fluid volumes (I.V.) and forboth scenarios associated to q D . There are three main observations from these figures. Thefirst one is the flow instability, or fingering, emerged as a result of the difference betweenthe injected fluid viscosity and the in-situ fluid viscosity. The second observation is that forthe higher injection rate scenario, the fingers are more developed at a later time comparedto that of the low injection scenario. The third one is that most of the fingers developed18nitially either merge or vanishes at the later stage, forming one or two main fingers. (a) (d) c q D = × - q D = × - p D = × p D = × m ( s ) (b) (c) m ( s ) m ( s ) m ( s ) Figure 2: Example 1: concentration fields, c ; I . V . = 42 m using ( a ) q D = 2 × − m / s at ∂ Ω and ( b ) q D = 1 × − m / s at ∂ Ω ; and I . V . = 180 m using ( c ) q D = 2 × − m / s at ∂ Ω and ( d ) q D = 1 × − m / s at ∂ Ω . The boundary conditions shown in this picture are corresponding to ∂ Ω and ∂ Ω of the massbalance equation (12). Next, we present the interaction among different processes including mechanical deforma-tion, calcite dissolution/precipitation, and viscosity alteration in Figure 3 for two differenttime steps. Note that the results of the low injection rate case are similar (for the samevolume of injected fluid); hence, we present here only the results of the high injection ratecase. 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 soliddeformation as the model is under compression. The increased fluid pressure by fluid injec-tion, however, limits the porosity reduction. This is reflected in Figures 3a and f in which ∂φ m ∂t is positive in the left part of the domain and negative in the right part of the domain.The ∂φ c ∂t result is shown in Figure 3c and h. Since the injected concentration c in = 0 . c (initial c eq ), the porous medium is dissolved in places to which the injected fluidis transported. Note that ∂φ c ∂t is positive where the dissolution occurs and negative wherethe perception occurs. At this time step, the maximum magnitude of ∂φ c ∂t is 10 − , which ismuch less compared to that of ∂φ m ∂t , which is around 10 − . We note this magnitude couldbe 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 frontprogresses. This alteration causes the flow instability discussed previously and establishesa preferential flow path. The impact of ∂φ m ∂t , ∂φ c ∂t , and µ alteration can be seen in q fieldshown in Figure 3e and j. Interestingly, as the first finger reaches the outlet boundary thesecond finger gradually disappears resulting in only one preferential path between the inletand outlet of the model.Thus, we have confirmed that the proposed framework can well simulate the expectedphysical and chemical phenomena including solid deformation, viscous fingering, and dis-solution/precipitation. The key ingredients of this method are the capability for tracking19he interface of the concentration species approximated by the high order methods withnumerical stabilization, the computation of reaction terms with the EG method, and thelocally conservative flux from BDM. (f)(h)(i) || v || × - (j) ∂ ø m / ∂ t ∂ ø c / ∂ t µ × - . × - . × - - . × - × - - × - (g) p . × × (a)(c)(d)(e)(b) m ( s) m ( s) Figure 3: Example 1: using q D = 2 × − m / s at ∂ Ω ; I . V . = 42 m : ( a ) the rate of change of porosityaccording to mechanics deformation, ∂φ m ∂t , see (23), ( b ) the fluid pressure, p , in surface and the displacement, u , in grey arrows, ( c ) the rate of change of porosity according to calcite dissolution/precipitation, ∂φ c ∂t , see(24), ( d ) the fluid viscosity, µ , and ( e ) the magnitude of the fluid velocity, (cid:107) q (cid:107) . I . V . = 480 m : ( f ) ∂φ m ∂t ,( g ) p in surface and u in grey arrows, ( h ) ∂φ c ∂t , ( i ) µ , and ( j ) (cid:107) q (cid:107) . Note that the magnitude of u is from 0 . . × − , and the trend of the results of the q D = 1 × − m / s at ∂ Ω case are similar. In the second example, we consider three layers (Ω =[0 , × [10 , =[0 , × [20 , =[0 , × [0 , ,20e set k = 8 . × − I m , while k = 8 . × − I m in Ω and Ω . Thus, in this case,the top Ω and bottom Ω layers have one order of magnitude of k less than that of themiddle layer Ω . All other rock and fluid parameters are the same as in the first example.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 porousmedium, even though the concentration fields at the early time are similar between thehigh, q D = 2 × − m / s, and low, q D = 1 × − m / s, injection rates, the progression ofconcentration field is different at the later time. It appears that the dynamic of the coupledprocesses controlled by the injection rate can impact the development of the dominantfinger in the middle layer. Note that since the top and bottom layers, Ω and Ω , havelower permeability than the middle layer, Ω , the flow mainly goes through the middlelayer. Similar to the previous example, one of the two initial fingers becomes the main pathconnecting the inlet and the outlet boundaries. (a) (d) c q D = × - q D = × - p D = × p D = × m ( s ) (b) (c) m ( s ) m ( s ) m ( s ) Figure 4: Example 2: concentration fields, c ; I . V . = 42 m using ( a ) q D = 2 × − m / s at ∂ Ω and ( b ) q D = 1 × − m / s at ∂ Ω ; and I . V . = 180 m using ( c ) q D = 2 × − m / s at ∂ Ω and ( d ) q D = 1 × − m / s at ∂ Ω . The boundary conditions shown in this picture are corresponding to ∂ Ω and ∂ Ω of the massbalance equation (12). In Figure 5, the behavior of the concentration and velocity fields, together with temporalporosity 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 , twofingers developed at the beginning, see Figure 5a and e. For the high injection rate, the topfinger, however, disappeared while the bottom finger progresses until it reaches the outlet ∂ Ω , see Figures 5b-d. One could see that the reaction front shown by ∂φ c ∂t progresses similarlyto the concentration front shown by the black contours. Besides, as the concentration fielddevelops, the change in µ enhances the flow channeling illustrated by velocity arrows. Forthe low injection rate case shown in Figure 5e-h, the development of the concentration fieldis dissimilar to that of the high injection rate case as the top finger becomes a preferablepath instead of the bottom one. Note that the dissolution and precipitation are shown inFigure 5 are a combined effect of injecting water that is unsaturated with respect to calcite21nd fluid pressure changes. It is clear that the majority of the dissolution occurs due to thetransport of the injected water in the porous domain. For the animated version of Figure5, please refer to Videos 1 and 2. These videos represent the flow and concentration field aswell as ∂φ c ∂t and illustrate the applicability of the presented coupled model for heterogeneousporous media.Importantly, this example has illustrated the capability of our proposed method—whichis equipped with the EG method—for handling discontinuous material properties acrossdifferent layers and the sharp interface of the concentration species. Moreover, we haveagain observed the expected physical and chemical phenomena, including solid deformation,viscous fingering, and dissolution/precipitation. (b) (f) || v || × - m (12000 s) ∂ø c / ∂ t × - × - - × - m (24000 s) m (52000 s) m (120000 s) (a)(c)(d) m (6000 s) m (12000 s) m (26000 s) m (60000 s) (e)(g)(h) q D = × - a t ∂ Ω q D = × - a t ∂ Ω Figure 5: Example 2: the results of the rate of change of porosity according to calcite dissolution/precipita-tion ( ∂φ c ∂t ) shown in surface plot, concentration ( c ), shown in black contour (10 contours ranging from 0.12 to1.6), and the fluid velocity ( q ) shown in arrows with q D = 2 × − m / s at ∂ Ω ; ( a ) I . V . = 36 m ( t = 6000s), ( b ) I . V . = 72 m ( t = 12000 s), ( c ) I . V . = 156 m ( t = 26000 s), and ( d ) I . V . = 360 m ( t = 60000 s),and with q D = 1 × − m / s at ∂ Ω ; ( e ) I . V . = 36 m ( t = 1200 s), ( f ) I . V . = 72 m ( t = 24000 s), ( g )I . V . = 156 m ( t = 32000 s), and ( h ) I . V . = 360 m ( t = 120000 s). .3. Example 3 In the given computational domain Ω = [0 , × [0 , k values as shown in Figure 1c. A random field generator [103] isutilized to generate a heterogeneous permeability field with a given mean permeability of k = 1 × − I m , variance of 0.5, and correlation lengths in x - and y -direction of 5 and1 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 HMCcoupled processes. Similar to the previous examples, two different injection rates are ap-plied. In Figure 6, the concentration fields are illustrated for two different injection rates atdifferent injected fluid volumes (I . V . = 42 m and I . V . = 180 m ). Unlike the two previousexamples, the preferential paths are established not only because of the flow instability re-sulting from the µ difference but also due to the high k channels inherited from the natureof heterogeneous porous media. During the early time, the concentration field of the high, q D = 2 × − m / s, and low, q D = 1 × − m / s, injection rate cases are similar, see Figure6a-b. The results of the concentration with the effects from the reaction are different at alater time (see Figure 6c-d). During the early time for both cases, the developed fingers fol-low the high permeable paths. At a later time, however, the results of the two scenarios arevery different. For the high injection rate case, the top finger continues developing while themiddle and the bottom fingers disappear. The result of the low injection rate case, however,shows that the top and the bottom fingers perish while the middle finger progresses. (a) (d) c q D = × - q D = × - p D = × p D = × m ( s ) (b) (c) m ( s ) m ( s ) m ( s ) Figure 6: Example 3: concentration fields, c ; I . V . = 42 m using ( a ) q D = 2 × − m / s at ∂ Ω and ( b ) q D = 1 × − m / s at ∂ Ω ; and I . V . = 180 m using ( c ) q D = 2 × − m / s at ∂ Ω and ( d ) q D = 1 × − m / s at ∂ Ω . The boundary conditions shown in this picture are corresponding to ∂ Ω and ∂ Ω of the massbalance equation (12). Figure 7 provides further insight into the reactive flow dynamics. It shows for bothinjection scenarios how the reaction fronts and flow fields evolve in time. As mentionedearlier, all the initial fingers at the beginning vanish except one that reaches the outlet ∂ Ω . The flow velocity field variations in time depict the emergence of the dominant finger.23ote that the magnitude of the mechanical deformation is higher than that of the calcitedissolution/precipitation and similar to what was observed in example 1. Therefore changesin porosity due to chemical reaction have a second-order effect on permeability compared tothat of induced by the mechanical deformation. Videos 3 and 4 representing the flow andconcentration field as well as ∂φ c ∂t illustrate the applicability of the presented coupled modelfor heterogeneous porous media. m ( s) m ( s) m ( s) m ( s) m ( s) m ( s) m ( s) m ( s) || v || × - ∂ ø c / ∂ t × - × - - × - (a)(b)(c)(d) (e)(f)(g)(h) q D = × - a t ∂ Ω q D = × - a t ∂ Ω Figure 7: Example 3: the results of the rate of change of porosity according to calcite dissolution/precipita-tion ( ∂φ c ∂t ) shown in surface plot, concentration ( c ), shown in black contour (10 contours ranging from 0.12 to1.6), and the fluid velocity ( q ) shown in arrows with q D = 2 × − m / s at ∂ Ω ; ( a ) I . V . = 36 m ( t = 6000s), ( b ) I . V . = 72 m ( t = 12000 s), ( c ) I . V . = 156 m ( t = 26000 s), and ( d ) I . V . = 360 m ( t = 60000 s),and with q D = 1 × − m / s at ∂ Ω ; ( e ) I . V . = 36 m ( t = 1200 s), ( f ) I . V . = 72 m ( t = 24000 s), ( g )I . V . = 156 m ( t = 32000 s), and ( h ) I . V . = 360 m ( t = 120000 s). Next, we investigate the local mass conservation property of the proposed framework inthe heterogeneous domain. The local mass conservation of each cell at each time step, r nmass ,is calculated byr nmass := (cid:90) T (cid:18) M + α K (cid:19) p n − p n − ∆ t n + αK σ nv − σ n − v ∆ t n + ˆ φ cn − φ n − c ∆ t n dV + (cid:88) e ∈E h (cid:90) e ¯ q n · n | e dS, (58)24nd the discrete numerical flux approximated by BDM, ¯ q n · n | e , is defined by¯ q n := q nh ∀ T ∈ T h , (59)¯ q n · n | e := − q D ∀ e ∈ E N,mh , (60)¯ q n · n | e := − q nh · n ∀ e ∈ E D,mh . (61)In Figure 8, the values of r nmass are illustrated for each case and time. One could seethat the magnitude of r nmass is always less than 1 × − , which is the tolerance set for thefixed-stress loop, see Algorithm 1; therefore, the framework is locally mass conservative. Wenote that the high injection rate case tends to the higher value of the magnitude of r nmass than that of the low injection rate case. m ( s) m ( s) m ( s) m ( s) m ( s) m ( s) m ( s) m ( s) | r m a ss | × - × - (a)(b)(c)(d) (e)(f)(g)(h) q D = × - a t ∂ Ω q D = × - a t ∂ Ω Lorem ipsum
Figure 8: Example 3: the illustration of the local mass conservative property with q D = 2 × − m / s at ∂ Ω ; ( a ) I . V . = 36 m ( t = 6000 s), ( b ) I . V . = 72 m ( t = 12000 s), ( c ) I . V . = 156 m ( t = 26000 s), and( d ) I . V . = 360 m ( t = 60000 s), and with q D = 1 × − m / s at ∂ Ω ; ( e ) I . V . = 36 m ( t = 1200 s), ( f )I . V . = 72 m ( t = 24000 s), ( g ) I . V . = 156 m ( t = 32000 s), and ( h ) I . V . = 360 m ( t = 120000 s). .4. Example 4 Lastly, we investigate the performance of the proposed framework when the permeabilityfield is anisotropic, and the grid is unstructured, as shown in Figure 1d. In the computationaldomain Ω = [0 , × [0 , k := (cid:20) k xx . . . k xx (cid:21) , (62)where k xx = 8 . × − m and all other parameters are similar to all other cases.Figure 9 shows the reactive flow dynamics and the residual of mass. We observe thatthe flow in the horizontal direction dominates the flow in the vertical direction since thepermeability 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 residualof mass values are always less than 1 × − , which is the tolerance set for the fixed-stressloop. (b) (e) || v || × - ∂ø c / ∂ t × - × - - × - (a)(c) m (30000 s) m (80000 s) m (120000 s) (d)(f) q D = × - a t ∂ Ω m (30000 s) m (80000 s) m (120000 s) | r m a ss | × - × - Figure 9: Example 4: the results of the rate of change of porosity according to calcite dissolution/precip-itation ( ∂φ c ∂t ) shown in surface plot, concentration ( c ), shown in black contour (10 contours ranging from0.12 to 1.6), and the fluid velocity ( q ) shown in arrows with q D = 2 × − m / s at ∂ Ω ; ( a ) I . V . = 180 m ( t = 30000 s), ( b ) I . V . = 480 m ( t = 80000 s), and ( c ) I . V . = 720 m ( t = 120000 s), and the local massconservative property; ( d ) I . V . = 180 m ( t = 30000 s), ( e ) I . V . = 480 m ( t = 80000 s), and ( f ) I . V . = 720m ( t = 120000 s). .5. 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 withdifferent injection rates. Besides, the injection flow rate also controls the developmentof the advection and reaction fronts.2. Using the applied set of the input parameters resulted in a more noticeable mechanicaleffect on the change in φ (and subsequently in k ) compared to that of the calcitedissolution/precipitation effect. We note that this observation could vary with differentsets 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 preferentialflow paths.3. The results of both homogeneous and heterogeneous as well as isotropic and anisotropicpermeability field show that our framework preserves mass locally. This property isessential for the coupled HMC system.In terms of computational efficiency, it is noted that the iteration number for the fixed-stress iteration was around three (four for the example 3) at the initial time stage, but itonly required two iterations for the rest of the time for all the presented examples. Forall examples, we have 31934, 23818, 7852, 11910 degrees of freedom for the displacement,flux, pressure, and concentration fields, respectively. The computational time was around4 . × − second per degrees of freedom per each time step. All simulations were computedon XeonE5 2650v4 with a single thread.
5. Conclusion
This paper has presented a mixed finite element framework for coupled hydro-mechanical-chemical processes in heterogeneous porous media. The main advantage of the proposedframework is its relatively affordable cost to attain local conservation regardless of materialanisotropy, thanks particularly to the use of the EG method. Through several numericalexamples, we have demonstrated the performance and capabilities of the proposed frameworkwith a focus on local conservation. The numerical results have highlighted how the overallbehavior is influenced by different processes, including solid deformation, calcite dissolution,and fluid viscosity alteration. The developed numerical model can provide insight into howthe interactions among HMC processes and heterogeneity manifest themselves at a largerscale. Future work includes an extension of the modeling framework to coupled thermo-hydro-mechanical-chemical processes in heterogeneous and/or fractured porous media.
6. Acknowledgements
This research has received financial support from the Danish Hydrocarbon Research andTechnology Centre under the Advanced Water Flooding program. The computational results27n this work have been produced by the multiphenics library [101], which is an extension ofFEniCS [100] for multiphysics problems. We acknowledge the developers of and contributorsto these libraries. TK also thanks the 2019 Computers & Geosciences Research grant forthe additional support. SL is supported by the National Science Foundation under GrantNo. NSF DMS-1913016. FB thanks Horizon 2020 Program for Grant H2020 ERC CoG2015 AROMA-CFD project 681447 that supported the development of multiphenics. JCacknowledges support from the Research Grants Council of Hong Kong (Project 27205918).
7. CRediT authorship contribution statementT. Kadeethum : Conceptualization, Formal analysis, Software, Validation, Writing -original draft, Writing - review & editing.
S. Lee : Conceptualization, Formal analysis,Supervision, Validation, Writing - review & editing.
F. Ballarin : Conceptualization, For-mal analysis, Software, Supervision, Validation, Writing - review & editing.
J. Choo :Conceptualization, Formal analysis, Supervision, Writing - review & editing.
H.M. Nick :Conceptualization, Funding acquisition, Supervision, Writing - review & editing.
8. Computer code availability
The scripts used to produce these results are available at this Git repository. The maindependencies are Numpy ( ≥ ≥ ≥ ≥ ≥ References [1] H. Nick, A. Raoof, F. Centler, M. Thullner, P. Regnier, Reactive dispersive contaminant transport incoastal aquifers: numerical simulation of a reactive henry problem, Journal of contaminant hydrology145 (2013) 90–104.[2] M. Hu, T. Hueckel, Environmentally enhanced crack propagation in a chemically degrading isotropicshale, G´eotechnique 63 (4) (2013) 313–321.[3] S. Pandey, A. Chaudhuri, S. Kelkar, V. Sandeep, H. Rajaram, Investigation of permeability alterationof fractured limestone reservoir due to geothermal heat extraction using three-dimensional thermo-hydro-chemical (THC) model, Geothermics 51 (2014) 46–62.[4] S. Pandey, A. Chaudhuri, The effect of heterogeneity on heat extraction and transmissivity evolutionin a carbonate reservoir: A thermo-hydro-chemical study, Geothermics 69 (2017) 45–54.[5] H. M. Nick, K.-H. Wolf, D. Brhun, Mixed CO –water injection into geothermal reservoirs: A numericalstudy, in: Proceedings of World Geothermal Congress, 2015, pp. 19–25.[6] J. Choo, W. Sun, Cracking and damage from crystallization in pores: Coupled chemo-hydro-mechanicsand phase-field modeling, Computer Methods in Applied Mechanics and Engineering 335 (2018) 347–349.[7] M. Tran, B. Jha, Coupling between transport and geomechanics affects spreading and mixing duringviscous fingering in deformable aquifers, Advances in Water Resources 136 (2020) 103485.[8] Z. Chen, Reservoir simulation: mathematical techniques in oil recovery, Vol. 77, Siam, 2007.[9] J. Du, R. Wong, Application of strain-induced permeability model in a coupled geomechanics-reservoirsimulator, Journal of Canadian Petroleum Technology 46 (12) (2007) 55–61.[10] J. Abou-Kassem, M. Islam, S. Farouq-Ali, Petroleum Reservoir Simulations, Elsevier, 2013.
11] T. Kadeethum, S. Salimzadeh, H. Nick, An investigation of hydromechanical effect on well productivityin fractured porous media using full factorial experimental design, Journal of Petroleum Science andEngineering 181 (2019) 106233.[12] T. Kadeethum, S. Salimzadeh, H. Nick, Well productivity evaluation in deformable single-fracturemedia, Geothermics 87 (2020).[13] M. Nejati, M. Dambly, M. Saar, A methodology to determine the elastic properties of anisotropic rocksfrom a single uniaxial compression test, Journal of Rock Mechanics and Geotechnical Engineering11 (6) (2019) 1166–1183.[14] S. Salimzadeh, E. Hagerup, T. Kadeethum, H. Nick, The effect of stress distribution on the shape anddirection of hydraulic fractures in layered media, Engineering Fracture Mechanics 215 (2019) 151–163.[15] J. Rutqvist, An overview of TOUGH-based geomechanics models, Computers & Geosciences 108(2017) 56–63.[16] M. Ahkami, A. Parmigiani, P. Di Palma, M. Saar, X. Kong, A lattice-boltzmann study of permeability-porosity relationships and mineral precipitation patterns in fractured porous media, ComputationalGeosciences (2020) 1–18.[17] C. Zhang, S. Zarrouk, R. Archer, A mixed finite element solver for natural convection in porous mediausing automated solution techniques, Computers & Geosciences 96 (2016) 181–192.[18] S. Dana, M. Wheeler, Convergence analysis of two-grid fixed stress split iterative scheme for coupledflow and deformation in heterogeneous poroelastic media, Computer Methods in Applied Mechanicsand Engineering 341 (2018) 788–806.[19] J. Kim, H. Tchelepi, R. Juanes, Stability and convergence of sequential methods for coupled flowand geomechanics: Fixed-stress and fixed-strain splits, Computer Methods in Applied Mechanics andEngineering 200 (13-16) (2011) 1591–1606.[20] H. Vik, S. Salimzadeh, H. Nick, Heat recovery from multiple-fracture enhanced geothermal systems:The effect of thermoelastic fracture interactions, Renewable Energy (2018).[21] J. White, R. I. Borja, Block-preconditioned Newton–Krylov solvers for fully coupled flow and geome-chanics, Computational Geosciences 15 (4) (2011) 647.[22] H. Nick, S. Matthai, A hybrid finite-element finite-volume method with embedded discontinuities forsolute transport in heterogeneous media, Vadose Zone Journal 10 (1) (2011) 299–312.[23] P. Salinas, D. Pavlidis, Z. Xie, H. Osman, C. Pain, M. Jackson, A discontinuous control volumefinite element method for multi-phase flow in heterogeneous porous media, Journal of ComputationalPhysics 352 (2018) 602–614.[24] Z. Chen, G. Huan, Y. Ma, Computational methods for multiphase flows in porous media, Vol. 2, Siam,2006.[25] T. Kadeethum, S. Lee, H. Nick, Finite element solvers for biot’s poroelasticity equations in porousmedia, Mathematical Geosciences (2020) 1–39.[26] T. Kadeethum, T. Jørgensen, H. Nick, Physics-informed neural networks for solving nonlinear diffu-sivity and Biot’s equations, PLoS ONE 15(5):e0232683 (2020).[27] T. Kadeethum, T. Jørgensen, H. Nick, Physics-informed Neural Networks for Solving Inverse Prob-lems of Nonlinear Biot’s Equations: Batch Training, in: 54th US Rock Mechanics/GeomechanicsSymposium, American Rock Mechanics Association, Golden, CO, USA, 2020.[28] K. Pruess, TOUGH user’s guide (1987).[29] J. Taron, D. Elsworth, Thermal-hydrologic-mechanical-chemical processes in the evolution of engi-neered geothermal reservoirs, International Journal of Rock Mechanics and Mining Sciences 46 (5)(2009) 855–864.[30] G. Danko, D. Bahrami, A new THMC model development for discrete-fracture EGS studies, Geother-mal Resources Council Transactions 36 (2012) 383–392.[31] J. Bean, M. Sanchez, J. Arguello, Sierra mechanics, an emerging massively parallel hpc capability,for use in coupled thmc analyses of hlw repositories in clay/shale, 5th International meeting Book ofabstracts (2012).[32] B. Aagaard, C. Williams, M. Knepley, PyLith: A finite-element code for modeling quasi-static and ynamic crustal deformation, Eos Trans. AGU 89 (53) (2008).[33] O. Kolditz, S. Bauer, L. Bilke, N. Bottcher, J. Delfs, T. Fischer, U. Gorke, T. Kalbacher,G. Kosakowski, C. McDermott, et al., OpenGeoSys: an open-source initiative for numerical sim-ulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media, EnvironmentalEarth Sciences 67 (2) (2012) 589–599.[34] A. Adam, D. Pavlidis, J. Percival, P. Salinas, R. Loubens, C. Pain, A. Muggeridge, M. Jackson, et al.,Dynamic mesh adaptivity for immiscible viscous fingering, in: SPE Reservoir Simulation Conference,Society of Petroleum Engineers, 2017.[35] Y. Melnikova, C. Jacquemyn, H. Osman, P. Salinas, G. Gorman, G. Hampson, M. Jackson, Reservoirmodelling using parametric surfaces and dynamically adaptive fully unstructured grids, in: ECMORXV-15th European Conference on the Mathematics of Oil Recovery, European Association of Geosci-entists & Engineers, 2016, p. 494.[36] A. Obeysekara, Q. Lei, P. Salinas, D. Pavlidis, J. Latham, J. Xiang, C. Pain, et al., A fluid-solidcoupled approach for numerical modeling of near-wellbore hydraulic fracturing and flow dynamicswith adaptive mesh refinement, in: 50th US Rock Mechanics/Geomechanics Symposium, AmericanRock Mechanics Association, 2016.[37] A. Obeysekara, Q. Lei, P. Salinas, D. Pavlidis, J. Xiang, J. Latham, C. Pain, Modelling stress-dependent single and multi-phase flows in fractured porous media based on an immersed-body methodwith mesh adaptivity, Computers and Geotechnics 103 (2018) 229–241.[38] J. H. Pr´evost, Dynaflow, Princeton University, Princeton, NJ 8544 (1983).[39] M. Tene, M. Al Kobaisi, H. Hajibeygi, Algebraic multiscale method for flow in heterogeneous porousmedia with embedded discrete fractures (F-AMS), Journal of Computational Physics 321 (2016) 819–845.[40] M. Cusini, A. Lukyanov, J. Natvig, H. Hajibeygi, Constrained pressure residual multiscale (CPR-MS)method for fully implicit simulation of multiphase flow in porous media, Journal of ComputationalPhysics 299 (2015) 472–486.[41] M. HosseiniMehr, C. Vuik, H. Hajibeygi, Adaptive dynamic multilevel simulation of fractured geother-mal reservoirs, Journal of Computational Physics: X (2020) 100061.[42] S. Matthai, S. Geiger, S. Roberts, A. Paluszny, M. Belayneh, A. Burri, A. Mezentsev, H. Lu,D. Coumou, T. Driesner, et al., Numerical simulation of multi-phase fluid flow in structurally complexreservoirs, Geological Society, London, Special Publications 292 (1) (2007) 405–429.[43] S. Salimzadeh, H. Nick, A coupled model for reactive flow through deformable fractures in enhancedgeothermal systems, Geothermics 81 (2019) 88–100.[44] E. Keilegavlen, R. Berge, A. Fumagalli, M. Starnoni, I. Stefansson, J. Varela, I. Berre, Porepy: Anopen-source software for simulation of multiphysics processes in fractured porous media, arXiv preprintarXiv:1908.09869 (2019).[45] B. Riviere, Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory andimplementation, SIAM, 2008.[46] S. Lee, Y. Lee, M. Wheeler, A locally conservative enriched Galerkin approximation and efficient solverfor elliptic and parabolic problems, SIAM Journal on Scientific Computing 38 (3) (2016) A1404–A1429.[47] K. Lipnikov, M. Shashkov, I. Yotov, Local flux mimetic finite difference methods, Numerische Math-ematik 112 (1) (2009) 115–152.[48] J. Choo, Large deformation poromechanics with local mass conservation: An enriched Galerkin finiteelement framework, International Journal for Numerical Methods in Engineering 116 (1) (2018) 66–90.[49] S. Sun, J. Liu, A locally conservative finite element method based on piecewise constant enrichmentof the continuous Galerkin method, SIAM Journal on Scientific Computing 31 (4) (2009) 2528–2548.[50] S. Lee, M. Wheeler, Enriched Galerkin methods for two-phase flow in porous media with capillarypressure, Journal of Computational Physics 367 (2018) 65–86.[51] J. Choo, S. Lee, Enriched Galerkin finite elements for coupled poromechanics with local mass conser-vation, Computer Methods in Applied Mechanics and Engineering 341 (2018) 311–332.[52] J. Choo, Stabilized mixed continuous/enriched Galerkin formulations for locally mass conservative oromechanics, Computer Methods in Applied Mechanics and Engineering 357 (2019) 112568.[53] P. Phillips, M. Wheeler, A coupling of mixed and continuous Galerkin finite element methods forporoelasticity I: the continuous in time case, Computational Geosciences 11 (2) (2007) 131.[54] P. Phillips, M. Wheeler, A coupling of mixed and continuous Galerkin finite element methods forporoelasticity II: the discrete-in-time case, Computational Geosciences 11 (2) (2007) 145–158.[55] J. Haga, H. Osnes, H. Langtangen, On the causes of pressure oscillations in low permeable and lowcompressible porous media, International Journal for Numerical and Analytical Methods in Geome-chanics 36 (12) (2012) 1507–1522.[56] M. Ferronato, N. Castelletto, G. Gambolati, A fully coupled 3-d mixed finite element model of biotconsolidation, Journal of Computational Physics 229 (12) (2010) 4813–4830.[57] B. Jha, R. Juanes, A locally conservative finite element framework for the simulation of coupled flowand reservoir geomechanics, Acta Geotechnica 2 (3) (2007) 139–153.[58] T. Kadeethum, H. Nick, S. Lee, F. Ballarin, Flow in porous media with low dimensional fractures byemploying Enriched Galerkin method, Advances in Water Resources (2020).[59] M. Biot, General theory of three-dimensional consolidation, Journal of applied physics 12 (2) (1941)155–164.[60] M. Biot, D. Willis, The elastic coefficients of the theory of consolidation, J. appl. Mech 15 (1957)594–601.[61] J. Choo, J. White, R. Borja, Hydromechanical modeling of unsaturated flow in double porosity media,International Journal of Geomechanics 16 (6) (2016) D4016002.[62] R. Borja, J. Choo, Cam-Clay plasticity, Part VIII: A constitutive framework for porous materials withevolving internal structure, Computer Methods in Applied Mechanics and Engineering 309 (2016) 653–679.[63] C. Macminn, E. Dufresne, J. Wettlaufer, Large deformations of a soft porous material, Physical ReviewApplied 5 (4) (2016) 1–30.[64] Y. Zhao, J. Choo, Stabilized material point methods for coupled large deformation and fluid flow inporous materials., Computer Methods in Applied Mechanics and Engineering 362 (2020) 112742.[65] J. Jaeger, N. G. Cook, R. Zimmerman, Fundamentals of rock mechanics, John Wiley & Sons, 2009.[66] O. Coussy, Poromechanics, John Wiley & Sons, 2004.[67] A. Chaudhuri, H. Rajaram, H. Viswanathan, Early-stage hypogene karstification in a mountain hy-drologic system: A coupled thermohydrochemical model incorporating buoyant convection, WaterResources Research 49 (9) (2013) 5880–5899.[68] A. Mikelic, M. Wheeler, Convergence of iterative coupling for coupled flow and geomechanics, Com-putational Geosciences 17 (3) (2013) 455–461.[69] B. Tjaden, S. Cooper, D. Brett, D. Kramer, P. Shearing, On the origin and application of the Brugge-man correlation for analysing transport phenomena in electrochemical systems, Current Opinion inChemical Engineering 12 (2016) 44–51.[70] D. Mu, Z. Liu, C. Huang, N. Djilali, Determination of the effective diffusion coefficient in porousmedia including knudsen effects, Microfluidics and Nanofluidics 4 (3) (2008) 257–260.[71] A. Raoof, H. Nick, S. M. Hassanizadeh, C. Spiers, Poreflow: A complex pore-network model forsimulation of reactive transport in variably saturated porous media, Computers & Geosciences 61(2013) 160–174.[72] F. Morel, J. Hering, Principles and applications of aquatic chemistry, John Wiley & Sons, 1993.[73] S. Dana, B. Ganis, M. Wheeler, A multiscale fixed stress split iterative scheme for coupled flow andporomechanics in deep subsurface reservoirs, Journal of Computational Physics 352 (2018) 1–22.[74] M. Taheriotaghsara, M. Bonto, A. Eftekhari, H. Nick, Prediction of oil breakthrough time in modifiedsalinity water flooding in carbonate cores, Fuel 274 (2020) 117806.[75] J. Rutqvist, Y. Wu, C. Tsang, G. Bodvarsson, A modeling approach for analysis of coupled multiphasefluid flow, heat transfer, and deformation in fractured porous rock, International Journal of RockMechanics and Mining Sciences 39 (4) (2002) 429–442.[76] J. Rutqvist, O. Stephansson, The role of hydromechanical coupling in fractured rock engineering, ydrogeology Journal 11 (1) (2003) 7–40.[77] K. Min, J. Rutqvist, C. Tsang, L. Jing, Stress-dependent permeability of fractured rock masses: anumerical study, International Journal of Rock Mechanics and Mining Sciences 41 (7) (2004) 1191–1210.[78] D. Grolimund, M. Elimelech, M. Borkovec, Aggregation and deposition kinetics of mobile colloidalparticles in natural porous media, Colloids and Surfaces A: Physicochemical and Engineering Aspects191 (1-2) (2001) 179–188.[79] B. Bijeljic, M. Blunt, Pore-scale modeling of transverse dispersion in porous media, Water ResourcesResearch 43 (12) (2007).[80] Y. Yortsos, D. Salin, On the selection principle for viscous fingering in porous media, Journal of FluidMechanics 557 (2006) 225–236.[81] A. Ern, A. Stephansen, P. Zunino, A discontinuous Galerkin method with weighted averages foradvection-diffusion equations with locally small and anisotropic diffusivity, IMA J. Numer. Anal.29 (2) (2009) 235–256.[82] A. Ern, A. Stephansen, A posteriori energy-norm error estimates for advection-diffusion equationsapproximated by weighted interior penalty methods, Journal of Computational Mathematics (2008)488–510.[83] R. Courant, K. Friedrichs, H. Lewy, On the partial difference equations of mathematical physics, IBMjournal of Research and Development 11 (2) (1967) 215–234.[84] Z. Ibrahim, K. Othman, M. Suleiman, Implicit r-point block backward differentiation formula forsolving first-order stiff ODEs, Applied Mathematics and Computation 186 (1) (2007) 558–565.[85] O. Akinfenwa, S. Jator, N. Yao, Continuous block backward differentiation formula for solving stiffordinary differential equations, Computers & Mathematics with Applications 65 (7) (2013) 996–1005.[86] S. Lee, A. Mikelic, M. Wheeler, T. Wick, Phase-field modeling of two phase fluid filled fractures in aporoelastic medium, Multiscale Modeling & Simulation 16 (4) (2018) 1542–1580.[87] J. Dormand, P. Prince, A family of embedded Runge-Kutta formulae, Journal of computational andapplied mathematics 6 (1) (1980) 19–26.[88] F. Brezzi, M. Fortin, Mixed and hybrid finite element methods, Vol. 15, Springer Science & BusinessMedia, 2012.[89] T. Kadeethum, H. Nick, S. Lee, C. Richardson, S. Salimzadeh, F. Ballarin, A Novel Enriched GalerkinMethod for Modelling Coupled Flow and Mechanical Deformation in Heterogeneous Porous Media,in: 53rd US Rock Mechanics/Geomechanics Symposium, American Rock Mechanics Association, NewYork, NY, USA, 2019.[90] R. Araya, E. Behrens, R. Rodriguez, An adaptive stabilized finite element scheme for the advection–reaction–diffusion equation, Applied Numerical Mathematics 54 (3-4) (2005) 491–503.[91] I. Harari, T. Hughes, Stabilized finite element methods for steady advection—diffusion with produc-tion, Computer Methods in Applied Mechanics and Engineering 115 (1-2) (1994) 165–191.[92] A. Masud, R. Khurram, A multiscale/stabilized finite element method for the advection–diffusionequation, Computer Methods in Applied Mechanics and Engineering 193 (21-22) (2004) 1997–2018.[93] E. Onate, Derivation of stabilized equations for numerical solution of advective-diffusive transportand fluid flow problems, Computer methods in applied mechanics and engineering 151 (1-2) (1998)233–265.[94] F. Brezzi, M. Bristeau, L. Franca, M. Mallet, G. Roge, A relationship between stabilized finite elementmethods and the Galerkin method with bubble functions, Computer Methods in Applied Mechanicsand Engineering 96 (1) (1992) 117–129.[95] A. Bonito, J. Guermond, B. Popov, Stability analysis of explicit entropy viscosity methods for non-linear scalar conservation equations, Mathematics of Computation 83 (287) (2014) 1039–1062.[96] J. Guermond, B. Popov, I. Tomas, Invariant domain preserving discretization-independent schemesand convex limiting for hyperbolic systems, Computer Methods in Applied Mechanics and Engineering347 (2019) 143–175.[97] G. Scovazzi, M. Wheeler, A. Mikelic, S. Lee, Analytical and variational numerical methods for unstable iscible displacement flows in porous media, Journal of Computational Physics 335 (2017) 444–496.[98] S. Lee, M. Wheeler, Adaptive enriched Galerkin methods for miscible displacement problems withentropy residual stabilization, Journal of Computational Physics 331 (2017) 19–37.[99] B. Riviere, M. Wheeler, A discontinuous Galerkin method applied to nonlinear parabolic equations,in: Discontinuous Galerkin methods, Springer, 2000, pp. 231–244.[100] M. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. Rognes,G. Wells, The FEniCS Project Version 1.5, Archive of Numerical Software 3 (100) (2015).[101] F. Ballarin, G. Rozza, multiphenics - easy prototyping of multiphysics problems in FEniCS (2019).URL https://mathlab.sissa.it/multiphenics [102] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Ei-jkhout, W. Gropp, D. Kaushik, M. Knepley, D. May, L. McInnes, R. Mills, T. Munson, K. Rupp,P. Sanan, B. Smith, S. Zampini, H. Zhang, H. Zhang, PETSc Users Manual, Tech. Rep. ANL-95/11- Revision 3.10, Argonne National Laboratory (2018).URL [103] H. Nick, R. Schotting, M. Gutierrez-Neri, K. Johannsen, Modeling transverse dispersion and variabledensity flow in porous media, Transport in porous media 78 (1) (2009) 11–35.[103] H. Nick, R. Schotting, M. Gutierrez-Neri, K. Johannsen, Modeling transverse dispersion and variabledensity flow in porous media, Transport in porous media 78 (1) (2009) 11–35.