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

T. Kadeethum  Technical University of Denmark, Denmark  Cornell University,
New York, USA<EMAIL_ADDRESS>S. Lee  Florida State University, Florida, USA
<EMAIL_ADDRESS>F. Ballarin  mathLab, Mathematics Area, SISSA, Italy
<EMAIL_ADDRESS>J. Choo  The University of Hong Kong, Hong Kong
<EMAIL_ADDRESS>H.M. Nick<EMAIL_ADDRESS>
###### Abstract

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.

###### keywords:

hydro-mechanical-chemical coupling , poroelasticity , reactive flow , mixed
finite element method , enriched Galerkin method , local conservation

††journal: Elsevierfn2fn2footnotetext: Author contribution statement is listed
in Section 7.fn2fn2footnotetext: This research has been awarded the funding
from the 2019 Computers & Geosciences Research grant.

##  1 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 ranging from groundwater and contaminant hydrology to
subsurface energy production [1, 2, 3, 4, 5, 6, 7]. The multiphysical
interactions in these problems give rise to strong heterogeneity in the
material properties. For instance, change in pore pressure perturbs effective
stress in the solid matrix, which can, in turn, alter the conductivity and
storability of the porous medium [8, 9, 10, 11, 12, 13]. Similarly, chemical
processes can result in the precipitation or dissolution of solid minerals,
which decreases or increases the pore volume, respectively, and thus, the
conductivity [14, 3, 4, 15, 16, 6]. Therefore, accurate numerical modeling of
coupled HMC problems requires a computational method that can robustly handle
strong heterogeneity in porous media.

Numerical simulation of multiphysical problems in porous media has been a
subject of extensive research (e.g. [17, 18, 19, 20, 21, 22, 23, 24, 25, 26,
27]), and lots of software packages have been developed for this 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 [28, 29,
15, 30], (2) SIERRA Mechanics, which has simulation capabilities for coupling
thermal, fluid, aerodynamics, solid mechanics and structural dynamics [31],
(3) PyLith, a finite-element code for modeling dynamic and quasi-static
simulations of coupled multiphysics processes [32], (4) OpenGeoSys project,
which is developed mainly based on the finite element method using object-
oriented programming THMC processes in porous media [33], (5) IC-FERST, a
reservoir simulator based on control-volume finite element methods and dynamic
unstructured mesh optimization [34, 35, 36, 37], (6) DYNAFLOW™, a nonlinear
transient finite element analysis platform [38], (7) DARSim, multiscale
multiphysics finite volume based simulator [39, 40, 41], (8) the CSMP, an
object-oriented application program interface, for the simulation of complex
geological processes, e.g. THMC, and their interactions [42, 43], and (9)
PorePy, an open-source modeling platform for multiphysics processes in
fractured porous media [44].

Nevertheless, it remains challenging to simulate coupled HMC processes in
porous media in a robust and efficient manner, especially when the material
properties are highly heterogeneous and/or anisotropic. Because HMC problems
involve transport phenomena in heterogeneous porous media, the numerical
method for these problems must ensure local (element-wise) conservation [45,
46]. The most practical method featuring local mass conservation may be the
finite volume method with a standard two-point flux approximation scheme.
However, this standard finite volume method requires the grid to be aligned
with the principal directions of the permeability/diffusivity tensors [47, 6],
which inhibits the 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 their
implementation is often complicated and onerous [48]. Discontinuous Galerkin
(DG) methods 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 the degrees of freedom.

In this paper, we present a new framework for computational modeling of
coupled HMC processes in porous media, which efficiently provides local mass
conservation even when the material properties are strongly heterogeneous and
anisotropic. The proposed framework combines two types of discretization
methods: (1) an enriched Galerkin (EG) method for reactive flow and transport,
and (2) a three-field mixed finite element method for coupled hydro-mechanical
processes. The EG method, which has recently been developed and advanced in
the literature [49, 46, 50, 51, 48, 52], augments a piecewise constant
function to the continuous Galerkin (CG) function space. This method uses the
same interior penalty type form as the DG method, but it requires a
substantially fewer number of degrees of freedom than the DG method. Thus the
EG method can provide locally conservative solutions to the reactive flow
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, 54, 55], which provides locally conservative, high-order
solutions to the fluid velocity field. Specifically, we employ the Lagrange
finite elements for approximating the displacement field, the Brezzi-Douglas-
Marini (BDM) element for the fluid velocity field, and the piecewise constant
element for the fluid pressure field. It is noted that this combination of
elements is our personal choice, and one may use another combination for the
same three primary fields as in [56, 57, 55].

The purpose of this work is to develop an accurate numerical method for
tackling coupled HMC processes in heterogeneous porous media, with a
practically affordable computational cost. Our specific objectives can be
summarized as follows:

  1. 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. 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. 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.

##  2 Governing equations

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

Let $\Omega\subset\mathbb{R}^{d}$ ($d\in\\{1,2,3\\}$) denote the computational
domain and $\partial\Omega$ denote the boundary. The time domain is denoted by
$\mathbb{T}=\left(0,\mathrm{T}\right]$ with $\mathrm{T}>0$. Primary variables
used in this paper are
$\bm{q}(\cdot,t):\Omega\times\mathbb{T}\to\mathbb{R}^{d}$, which is a vector-
valued Darcy velocity ($\mathrm{m}\mathrm{/}\mathrm{s}$),
$p(\cdot,t):\Omega\times\mathbb{T}\to\mathbb{R}$, which is a scalar-valued
fluid pressure ($\mathrm{P}\mathrm{a}$),
$\bm{u}(\cdot,t):\Omega\times\mathbb{T}\to\mathbb{R}^{d}$, which is a vector-
valued displacement ($\mathrm{m}$),
$c_{i}:\Omega\times\mathbb{T}\rightarrow\mathbb{R}$, which is the $i$-th
component of chemical concentration
($\mathrm{m}\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{/}\mathrm{m}^{3}$).

###  2.1 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 quasi-static
conditions, the linear momentum balance equation can be written as

$\nabla\cdot\bm{\sigma}(\bm{u},p)+\bm{f}=\bm{0},$ (1)

where $\bm{f}$ is the body force term defined as
$\rho\phi\mathbf{g}+\rho_{s}(1-\phi)\mathbf{g}$, where $\rho$ is the fluid
density, $\rho_{s}$ is the solid density, $\phi$ is the porosity, $\mathbf{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, $\bm{\sigma}$ is the
total stress tensor, which may be related to the effective stress tensor
$\bm{\sigma}^{\prime}$ and the pore pressure $p$ as

$\bm{\sigma}(\bm{u},p)=\bm{\sigma}^{\prime}(\bm{u})-\alpha p\mathbf{I}.$ (2)

Here, $\mathbf{I}$ is the second-order identity tensor, and $\alpha$ is the
Biot coefficient defined as [65]:

$\alpha=1-\frac{K}{K_{{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 with the displacement vector, which can
be written as

$\bm{\sigma}^{\prime}(\bm{u})=\lambda_{l}\operatorname{tr}(\bm{\varepsilon}(\bm{u}))\mathbf{I}+2\mu_{l}\bm{\varepsilon}{(\bm{u})}.$
(4)

Here, $\bm{\varepsilon}$ is the infinitesimal strain tensor, defined as

$\bm{\varepsilon}(\bm{u}):=\frac{1}{2}\left(\nabla\bm{u}+(\nabla\bm{u})^{\intercal}\right),$
(5)

and $\lambda_{l}$ and $\mu_{l}$ are the Lamé constants, which are related to
the bulk modulus and the Poisson ratio $\nu$ of the solid matrix as

$\lambda_{l}=\frac{3K\nu}{1+\nu},\text{ and
}\mu_{l}=\frac{3K(1-2\nu)}{2(1+\nu)}.$ (6)

For this solid deformation problem, the domain boundary $\partial\Omega$ is
assumed to be suitably decomposed into displacement and traction boundaries,
$\partial\Omega_{u}$ and $\partial\Omega_{t}$, respectively. Then the linear
momentum balance equation is supplemented by the boundary and initial
conditions as:

$\begin{split}\nabla\cdot\bm{\sigma}^{\prime}(\bm{u})+\alpha\nabla\cdot\left(p\mathbf{I}\right)+\bm{f}=\bm{0}&\text{
\> in \> }\Omega\times\mathbb{T},\\\ \bm{u}=\bm{u}_{D}&\text{ \> on \>
}\partial\Omega_{u}\times\mathbb{T},\\\
\bm{\sigma}{(\bm{u})}\cdot\mathbf{n}=\bm{t}_{D}&\text{ \> on \>
}\partial\Omega_{t}\times\mathbb{T},\\\ \bm{u}=\bm{u}_{0}&\text{ \> in \>
}\Omega\text{ at }t=0,\end{split}$ (7)

where $\bm{u}_{D}$ and ${\bm{t}_{D}}$ are prescribed displacement and traction
values at the boundaries, respectively, and $\mathbf{n}$ is the unit normal
vector to the boundary.

Next, the mass balance equation is given as [66, 19, 4, 43]:

$\frac{1}{M}\dfrac{\partial p}{\partial
t}+\alpha\frac{\partial{\varepsilon_{v}}}{\partial
t}+\dfrac{\partial\phi_{c}}{\partial t}+\nabla\cdot\bm{q}=g\text{ in
}\Omega\times\mathbb{T},$ (8)

where

$\frac{1}{M}=\left(\phi_{0}c_{f}+\dfrac{\alpha-\phi_{0}}{K_{s}}\right)$ (9)

is the Biot modulus. Here, $c_{f}$ is the fluid compressibility, $\phi_{0}$ is
the initial porosity, ${\varepsilon_{v}}$ :=
$\operatorname{tr}(\bm{\varepsilon})=\nabla\cdot\bm{u}$ is the volumetric
strain, and $g$ is a sink/source term. Because we will introduce chemical
effects later on, we have added $\frac{\partial\phi_{c}}{\partial t}$ to the
standard poroelasticity equation [67, 3, 4, 43]. This term will be discussed
again after introducing chemical effects. Also, $\bm{q}$ is the superficial
velocity vector, which is given by Darcy’s law as

$\bm{q}=-\frac{\bm{k}(\phi)}{\mu(c_{i})}(\nabla p-\rho\mathbf{g}).$ (10)

Note that here the fluid viscosity $\mu$ is considered a function of
concentration $c_{i}$. Again, the gravitational force, $\rho\mathbf{g}$, will
be neglected in this work, without loss of generality. In addition,
$\bm{k}(\phi)$ is the matrix permeability tensor defined as

$\bm{k}:=\begin{cases}k_{mult}(\phi)\left[\begin{array}[]{lll}{{k}^{xx}}&{{k}^{xy}}&{{k}^{xz}}\\\
{{k}^{yx}}&{{k}^{yy}}&{{k}^{yz}}\\\
{{k}^{zx}}&{{k}^{zy}}&{k}^{zz}\end{array}\right]&\text{if}\ d=3,\\\ \\\
k_{mult}(\phi)\left[\begin{array}[]{ll}{{k}^{xx}}&{{k}^{xy}}\\\
{{k}^{yx}}&{{k}^{yy}}\\\ \end{array}\right]&\text{if}\ d=2,\\\ \\\
k_{mult}(\phi)\ k&\text{if}\ d=1,\end{cases}$ (11)

The $k^{xx}$, $k^{yy}$, and $k^{zz}$ represent the matrix permeability in
$x$-, $y$-, and $z$-direction, respectively. The $k_{mult}(\phi)$ is a
multiplier used to update $\bm{k}$ when $\phi$ is altered, which will be
described later.

For the fluid flow problem, the domain boundary $\partial\Omega$ is also
suitably decomposed into the pressure and flux boundaries,
$\partial\Omega_{p}$ and $\partial\Omega_{q}$, respectively. In what follows,
we apply the fixed stress split scheme [19, 68], assuming
$\left(\sigma_{v}-\sigma_{v,0}\right)+\alpha\left(p-p_{0}\right)=K\varepsilon_{v}$.
Then we write the fluid flow problem with boundary and initial conditions as

$\begin{split}\left(\frac{1}{M}+\frac{\alpha^{2}}{K}\right)\frac{\partial
p}{\partial t}+\frac{\alpha}{K}\frac{\partial\sigma_{v}}{\partial
t}+\dfrac{\partial\phi_{c}}{\partial t}+\nabla\cdot\bm{q}=g&\text{ \> in \>
}\Omega\times\mathbb{T},\\\ p=p_{D}&\text{ \> on \>
}\partial\Omega_{p}\times\mathbb{T},\\\ \bm{q}\cdot\mathbf{n}=q_{D}&\text{ \>
on \>}\partial\Omega_{q}\times\mathbb{T},\\\ p=p_{0}&\text{ \> in \>
}\Omega\text{ at }t=0,\end{split}$ (12)

where $\sigma_{v}:=\frac{1}{3}\operatorname{tr}(\bm{\sigma})$ is the
volumetric stress, and $p_{D}$ and $q_{D}$ are the given boundary pressure and
flux, respectively.

###  2.2 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,\ldots,N_{c}$,

$\frac{\partial}{\partial t}(\phi
c_{i})+\nabla\cdot\eta(\bm{q},c_{i})=q_{i}(c_{i}),\ \mbox{ in
}\Omega\times{\mathbb{T}},$ (13)

where $q_{i}(c_{i})$ is a reaction term coupled with sink/source for each
component, and the mass flux $\eta(\bm{q},c_{i})$ is defined as

$\eta(\bm{q},c_{i}):=\bm{q}c_{i}-{\bm{D}_{e,i}}(\phi)\nabla c_{i}.$ (14)

Here ${\bm{D}_{e,i}}(\phi)$ is the effective diffusion coefficient tensor
defined as

${\bm{D}_{e,i}}:=\frac{\phi}{\tau}\bm{D}_{i},$ (15)

where $\tau=\phi^{-\frac{1}{2}}$ [69, 70] and $\bm{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
$\partial\Omega_{\rm in}$ and $\partial\Omega_{\rm out}$, respectively, which
are defined as

$\partial\Omega_{\rm
in}:=\\{\bm{x}\in\partial\Omega:\bm{q}\cdot\mathbf{n}<0\\}\quad\mbox{ and
}\quad\partial\Omega_{\rm
out}:=\\{\bm{x}\in\partial\Omega:\bm{q}\cdot\mathbf{n}\geq 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:
$\left\\{\mathrm{H}^{+}\right.$, $\mathrm{Ca}^{2+}$, $\mathrm{CaHCO}_{3}^{+}$,
$\mathrm{OH}^{-}$, $\mathrm{CO}_{3}^{2-}$, $\mathrm{HCO}_{3}^{-}$,
$\mathrm{H}_{2}\mathrm{CO}_{3}^{*}$, $\mathrm{CaCO}_{3}^{*}(\mathrm{Aq})\\}$
[67, 3, 71, 72]. For simplicity, in this paper we consider a reduced system
based on the empirical relationship presented in [67, 3, 4], 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:

$\begin{split}\frac{\partial}{\partial t}(\phi
c)+\nabla\cdot\left(\bm{q}c\right)-\nabla\cdot\left(\bm{D}_{e}(\phi)\nabla
c\right)=q&\mbox{ \> in \> }\Omega\times(0,\mathbb{T}],\\\
\eta(\bm{q},c)\cdot\mathbf{n}={c_{in}}\bm{q}\cdot\mathbf{n}&\mbox{ \> on \>
}\partial\Omega_{{\rm in}}\times(0,\mathbb{T}],\\\ {\bm{D}_{e}}(\phi)\nabla
c\cdot\mathbf{n}=0&\mbox{ \> on \> }\partial\Omega_{{\rm
out}}\times(0,\mathbb{T}],\\\ c={c_{0}}&\text{ \> in \> }\Omega\text{ at
}t=0,\end{split}$ (17)

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
[67, 3, 4], 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 calculated as

$R_{c}=\left\\{\begin{array}[]{cl}10^{r},&\text{ for }\widetilde{c}>0,\\\
-10^{r},&\text{ for }\widetilde{c}<0,\\\ 0.0,&\text{ for
}\widetilde{c}=0,\end{array}\right.$ (19)

with

$\widetilde{c}=\frac{\left(c_{eq}-c\right)}{c_{eq}},$ (20)

$\displaystyle
r=a_{0}+a_{1}\tau+a_{2}\log|\widetilde{c}|+a_{3}\tau^{2}+a_{4}\tau\log|\widetilde{c}|+a_{5}(\log|\widetilde{c}|)^{2},$
(21)

and

$\displaystyle c_{eq}$ $\displaystyle=1.417\times 10^{-3}+3.823\times
10^{-6}p-4.313\times 10^{-7}\tau$ (22) $\displaystyle-2.148\times
10^{-8}p^{2}+4.304\times 10^{-8}p\tau-7.117\times 10^{-8}\tau^{2}.$

Here, $\tau$ is the medium temperature, and $a_{0}$, $a_{1}$, $a_{2}$,
$a_{3}$, $a_{4}$, and $a_{5}$ are defined in Table 1.

Table 1: Coefficients of the (21) for different range of (20) | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$
---|---|---|---|---|---|---
$\widetilde{c}>0.01$ | -5.73 | $1.25\times 10^{-2}$ | 1.38 | $2.61\times 10^{-5}$ | $-4.01\times 10^{-3}$ | $3.26\times 10^{-1}$
$-0.01<\widetilde{c}$ | -6.45 | $2.09\times 10^{-2}$ | $-4.65\times 10^{-2}$ | $3.06\times 10^{-5}$ | $9.25\times 10^{-3}$ | $-4.59\times 10^{-1}$
$-0.01<\widetilde{c}\leq 0.01$ | -5.80 | $1.35\times 10^{-2}$ | $9.97\times 10^{-1}$ | $3.80\times 10^{-5}$ | $1.51\times 10^{-5}$ | $-4.87\times 10^{-4}$

Before closing this section, we describe physical properties that are coupled
with primary variables, $\bm{u}$, $\bm{q}$, $p$, and $c$. The porosity change
due to solid deformation may be expressed as [59, 73, 18]:

$\begin{split}\phi_{m}=\phi_{0}+\left(\alpha-\phi_{0}\right)\left(\epsilon_{v}-\epsilon_{v_{0}}\right)+\frac{\left(\alpha-\phi_{0}\right)(1-\alpha)}{K}\left(p-p_{0}\right),\end{split}$
(23)

where $\epsilon_{v_{0}}$ is the initial volumetric strain. The porosity
alteration due to calcite dissolution/precipitation is calculated as

$\Gamma\left(\bm{u},c\right)=\frac{\partial\phi_{c}}{\partial
t}=\frac{R_{c}A_{s}}{\rho_{s}\omega},$ (24)

where $\omega$ 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
$\rho_{s}=2500$ $\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$ is used for
throughout this paper. Note that this term, (24), enters (12). Also, the terms
$\phi_{m}$ and $\phi_{c}$ are used to distinguish between the changes in
$\phi$ 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

$A_{s}=A_{0}\frac{\phi}{\phi_{0}}\frac{\log(\phi)}{\log\left(\phi_{0}\right)},$
(25)

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]:

$\bm{k}=\bm{k}_{0}k_{mult}(\phi)=\bm{k}_{0}\exp\left(b\left(\frac{\phi}{\phi_{0}}-1\right)\right),$
(26)

where $\bm{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 $\mu$, and we adopt the specific form
from [78, 79, 80], given by

$\mu=\log\left(\mu_{l}\right)+\left(\frac{c-c_{l}}{c_{h}-c_{l}}\right)\left(\log\left(\mu_{h}\right)-\log\left(\mu_{l}\right)\right),$
(27)

where $c_{l}$ and $c_{h}$ are lower and higher bounds of the concentration,
and $\mu_{l}$ and $\mu_{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.

Table 2: Summary of the effects of individual physical processes on physical properties Physical properties | Mechanical deformation | Fluid pressure | Calcite concentration
---|---|---|---
$\phi$ | (23) | (23) | (24)
$\bm{k}$ | (23) + (26) | (23) + (26) | (24) + (26)
$\mu$ | - | - | (27)
$\bm{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 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.

###  3.1 Domain discretization and geometrical quantities

We begin by introducing the notations used throughout this paper. Let
$\mathcal{T}_{h}$ be a shape-regular triangulation obtained by a partition of
$\Omega$ into $d$-simplices (triangles in $d=2$, tetrahedra in $d=3$). For
each cell $T\in\mathcal{T}_{h}$, we denote by $h_{T}$ the diameter of $T$, and
we set $h=\max_{T\in\mathcal{T}_{h}}h_{T}$ and
$h_{l}=\min_{T\in\mathcal{T}_{h}}h_{T}$. We further denote by
$\mathcal{E}_{h}$ the set of all faces (i.e., $d-1$ dimensional entities
connected to at least a $T\in\mathcal{T}_{h}$) and by $\mathcal{E}_{h}^{I}$
and $\mathcal{E}_{h}^{\partial}$ the collection of all interior and boundary
facets, respectively. The boundary set $\mathcal{E}_{h}^{\partial}$ 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,
$\mathcal{E}_{h}^{D,u}$ and $\mathcal{E}_{h}^{N,u}$ correspond to the faces on
$\partial\Omega_{u}$ and $\partial\Omega_{tr}$, respectively, for (7). On the
other hand, for (12), $\mathcal{E}_{h}^{D,m}$ and $\mathcal{E}_{h}^{N,m}$
conform to $\partial\Omega_{p}$ and $\partial\Omega_{q}$, respectively.
Lastly, for (17), $\mathcal{E}_{h}^{\partial}$ is decomposed into
$\mathcal{E}_{h}^{\mathrm{In}}$ and $\mathcal{E}_{h}^{\mathrm{Out}}$.

We also define

$e=\partial T^{+}\cap\partial T^{-},\ \ e\in\mathcal{E}_{h}^{I},$

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}:=\frac{\operatorname{meas}\left(T^{+}\right)+\operatorname{meas}\left(T^{-}\right)}{2\operatorname{meas}(e)},$
(28)

depending on the argument, meas($\cdot$) represents the measure of a cell or
of a facet.

Let $\mathbf{n}^{+}$ and $\mathbf{n}^{-}$ be the outward unit normal vectors
to $\partial T^{+}$ and $\partial T^{-}$, respectively. For any given scalar
function $\zeta:\mathcal{T}_{h}\to\mathbb{R}$ and vector function
$\bm{\tau}:\mathcal{T}_{h}\to\mathbb{R}^{d}$, we denote by $\zeta^{\pm}$ and
$\bm{\tau}^{\pm}$ the restrictions of $\zeta$ and $\bm{\tau}$ to $T^{\pm}$,
respectively. Subsequently, we define the weighted average operator as

$\\{\zeta\\}_{\delta
e}=\delta_{e}\zeta^{+}+\left(1-\delta_{e}\right)\zeta^{-},\ \text{ on
}e\in\mathcal{E}_{h}^{I},$ (29)

and

$\\{\bm{\tau}\\}_{\delta
e}=\delta_{e}\bm{\tau}^{+}+\left(1-\delta_{e}\right)\bm{\tau}^{-},\ \text{ on
}e\in\mathcal{E}_{h}^{I},$ (30)

where $\delta_{e}$ is calculated by [81, 82]:

$\delta_{e}:=\frac{{k}^{-}_{e}}{{k}^{+}_{e}+{k}^{-}_{e}}.$ (31)

Here,

${k}^{+}_{e}:=\left(\mathbf{n}^{+}\right)^{\intercal}\cdot\bm{k}^{+}\mathbf{n}^{+},\
\text{ and
}{k}^{-}_{e}:=\left(\mathbf{n}^{-}\right)^{\intercal}\cdot\bm{k}^{-}\mathbf{n}^{-},$
(32)

where ${k_{e}}$ is a harmonic average of $k^{+}_{e}$ and ${k}^{-}_{e}$ which
reads

${k_{e}}:=\frac{2{k}^{+}_{e}{k}^{-}_{e}}{{k}^{+}_{e}+{k}^{-}_{e}},$ (33)

and $\bm{k}$ is defined as in (11). The jump across an interior edge will be
defined as

$\displaystyle\left[\\!\left[\zeta\right]\\!\right]=\zeta^{+}\mathbf{n}^{+}+\zeta^{-}\mathbf{n}^{-}\quad\mbox{
and
}\quad\left[\\!\left[{\bm{\tau}}\right]\\!\right]=\bm{\tau}^{+}\cdot\mathbf{n}^{+}+\bm{\tau}^{-}\cdot\mathbf{n}^{-}\quad\mbox{on
}e\in\mathcal{E}_{h}^{I}.$

Finally, for $e\in\mathcal{E}^{\partial}_{h}$, we set
$\left\\{\zeta\right\\}_{\delta_{e}}:=\zeta$ and
$\left\\{\bm{\tau}\right\\}_{\delta_{e}}:=\bm{\tau}$ for what concerns the
definition of the weighted average operator, and
$\left[\\!\left[\zeta\right]\\!\right]:=\zeta\mathbf{n}$ and
$\left[\\!\left[\bm{\tau}\right]\\!\right]:=\bm{\tau}\cdot\mathbf{n}$ as
definition of the jump operator.

###  3.2 Temporal discretization

The time domain $\mathbb{T}=\left(0,\mathrm{T}\right]$ is partitioned into $N$
subintervals such that $0=:t^{0}<t^{1}<\cdots<t^{N}:=\mathrm{T}$. The length
of each subinterval $\Delta t^{n-1}$ is defined as $\Delta
t^{n-1}=t^{n}-t^{n-1}$ where $n$ represents the current time step. We assume
that the user provides the initial $\Delta t^{0}$, while an adaptive procedure
is carried out to choose $\Delta t^{n-1}$, $n>1$, as follows:

$\Delta
t^{n-1}:=\begin{cases}\mathrm{CFL}\frac{h_{l}}{\left\|\bm{q}^{n-1}\right\|_{\infty}}&\text{if}\
\Delta t^{n}\leq\Delta t_{\max}\\\ \Delta t_{\max}&\text{if}\ \Delta
t^{n}>\Delta t_{\max},\end{cases}$ (34)

where $\mathrm{CFL}$ is a constant that the user can provide according to the
Courant-Friedrichs-Lewy condition [83], $\left\|\cdot\right\|_{\infty}$ is the
maximum norm of a vector function, and $\Delta t_{\max}$ is a maximum allowed
time step. Note that we use $\Delta t_{\max}$ as a tool to control $\Delta
t^{n}$ as the model approaches a steady-state condition since
$\left\|\bm{q}^{n-1}\right\|_{\infty}$ may approach zero, which would lead to
a very large ratio $\frac{h_{l}}{\left\|\bm{q}^{n-1}\right\|_{\infty}}$.

Let $\varphi(\cdot,t)$ be a scalar function and $\varphi^{n}$ be its
approximation at time $t^{n}$, i.e.
$\varphi^{n}\approx\varphi\left(t^{n}\right)$. We employ the following
backward differentiation formula [84, 85, 86]

$\mathrm{BDF}_{m}\left(\varphi^{n}\right):=\left\\{\begin{array}[]{ll}\frac{1}{\Delta
t^{n}}\left(\varphi^{n}-\varphi^{n-1}\right)&m=1\\\ \frac{1}{2\Delta
t^{n}}\left(3\varphi^{n}-4\varphi^{n-1}+\varphi^{n-2}\right)&m=2\\\
\frac{1}{6\Delta
t^{n}}\left(11\varphi^{n}-18\varphi^{n-1}+9\varphi^{n-2}-2\varphi^{n-3}\right)&m=3\\\
\frac{1}{12\Delta
t^{n}}\left(25\varphi^{n}-48\varphi^{n-1}+36\varphi^{n-2}-16\varphi^{n-3}+3\varphi^{n-4}\right)&m=4\end{array}\right.$
(35)

for the discretization of the time derivative of $\varphi(\cdot,t)$ at time
$t^{n}$. We also utilize the explicit Runge-Kutta methods [87, 24]:

$\mathrm{RK}_{1}(\varphi^{n})=\varphi^{n+1}=\varphi^{n}+\kappa_{1},$ (36)

$\begin{array}[]{l}\kappa_{1}=\Delta
t^{n}F\left(\mathbb{X}^{n},\mathbb{Y}^{n}\right),\end{array}$

for the first order Runge-Kutta method corresponding to the explicit Euler
method, and

$\mathrm{RK}_{4}(\varphi^{n})=\varphi^{n+1}=\varphi^{n}+\frac{1}{6}\kappa_{1}+\frac{1}{3}\kappa_{2}+\frac{1}{3}\kappa_{3}+\frac{1}{6}\kappa_{4},$
(37)

$\begin{array}[]{l}\kappa_{1}=\Delta
t^{n}F\left(\mathbb{X}^{n},\mathbb{Y}^{n}\right),\\\ \kappa_{2}=\Delta
t^{n}F\left(\mathbb{X}^{n}+\frac{1}{2}\Delta
t^{n},\mathbb{Y}^{n}+\frac{1}{2}\kappa_{1}\right),\\\ \kappa_{3}=\Delta
t^{n}F\left(\mathbb{X}^{n}+\frac{1}{2}\Delta
t^{n},\mathbb{Y}^{n}+\frac{1}{2}\kappa_{2}\right),\\\ \kappa_{4}=\Delta
t^{n}F\left(\mathbb{X}^{n}+\Delta
t^{n},\mathbb{Y}^{n}+\kappa_{3}\right),\end{array}$

for the forth order Runge-Kutta method,
$F\left(\mathbb{X}^{n},\mathbb{Y}^{n}\right)$ is any functions with
independent variable $\mathbb{X}$ and dependent variable $\mathbb{Y}$ [87, 8],
which we will specify in the linearization and solving processes in Section
3.4.

Finally, we define an extrapolation operator as follows [24, 8]:

$\mathrm{EX}\left(\varphi\right)=\hat{\varphi}^{n+1}=\begin{cases}\left(1+\dfrac{\Delta
t^{n}}{\Delta t^{n-1}}\right){\varphi}^{n}-\dfrac{\Delta t^{n}}{\Delta
t^{n-1}}{\varphi}^{n-1}&\text{if}\ n\geq 1,\\\ {\varphi}^{0}&\text{if}\
n=0,\end{cases}$ (38)

and in the following we will adopt the notation $\hat{\varphi}^{n+1}$ to
denote an extrapolation value of $\\{\varphi^{n},\varphi^{n-1}\\}$.

###  3.3 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:

$\mathcal{U}_{h}^{\mathrm{CG}_{k}}\left(\mathcal{T}_{h}\right):=\left\\{\bm{\psi_{u}}\in\mathbb{C}^{0}(\Omega{;\mathbb{R}^{d}}):\left.\bm{\psi_{u}}\right|_{T}\in\mathbb{P}_{k}(T{;\mathbb{R}^{d}}),\forall
T\in\mathcal{T}_{h}\right\\},$ (39)

where $\mathbb{C}^{0}(\Omega{;\mathbb{R}^{d}})$ denotes the space of vector-
valued piecewise continuous polynomials, $\mathbb{P}_{k}(T{;\mathbb{R}^{d}})$
is the space of polynomials of degree at most $k$ over each element $T$, and
$\bm{\psi_{u}}$ denotes a generic function of
$\mathcal{U}_{h}^{\mathrm{CG}_{k}}\left(\mathcal{T}_{h}\right)$. In addition,
the CG space for scalar-valued functions is defined as:

$\mathcal{P}_{h}^{\mathrm{CG}_{k}}\left(\mathcal{T}_{h}\right):=\left\\{\psi_{p}\in\mathbb{C}^{0}(\Omega):\left.\psi_{p}\right|_{T}\in\mathbb{P}_{k}(T),\forall
T\in\mathcal{T}_{h}\right\\},$ (40)

where $\mathbb{C}^{0}(\Omega):=\mathbb{C}^{0}(\Omega;\mathbb{R})$ and
$\mathbb{P}_{k}(T):=\mathbb{P}_{k}(T{;\mathbb{R}})$. Next, we define the
following DG function space:

$\mathcal{P}_{h}^{\mathrm{DG}_{k}}\left(\mathcal{T}_{h}\right):=\left\\{\psi_{p}\in
L^{2}(\Omega):\left.\psi_{p}\right|_{T}\in\mathbb{P}_{k}(T),\forall
T\in\mathcal{T}_{h}\right\\},$ (41)

where $L^{2}(\Omega)$ 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:

$\mathcal{P}_{h}^{\mathrm{EG}_{k}}\left(\mathcal{T}_{h}\right):=\mathcal{P}_{h}^{\mathrm{CG}_{k}}\left(\mathcal{T}_{h}\right)\oplus\mathcal{P}_{h}^{\mathrm{DG}_{0}}\left(\mathcal{T}_{h}\right),$
(42)

i.e., a CG finite element space enriched by the space
$\mathcal{P}_{h}^{\mathrm{DG}_{0}}\left(\mathcal{T}_{h}\right)$ of piecewise
constant functions. In the following we denote ${\psi_{c}}$ a generic function
of $\mathcal{P}_{h}^{\mathrm{EG}_{k}}\left(\mathcal{T}_{h}\right)$.

Lastly, we define the BDM function space as follows [88]:

$\mathcal{V}_{h}^{\mathrm{BDM}_{k}}\left(\mathcal{T}_{h}\right):=\left\\{\bm{\psi_{v}}\in
H(\operatorname{div},\Omega):\left.\bm{\psi_{v}}\right|_{T}\in\mathrm{BDM}(T),\forall
T\in\mathcal{T}_{h}\right\\}$ (43)

where $\bm{\psi_{v}}$ denotes a generic function of
$\mathcal{V}_{h}^{\mathrm{BDM}_{k}}\left(\mathcal{T}_{h}\right)$ and
$\mathrm{BDM}(T)$ is defined according to [88].

####  3.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 approximated displacement solution
$\bm{u_{h}}\in\mathcal{U}_{h}^{\mathrm{CG}_{2}}\left(\mathcal{T}_{h}\right)$
as done in [51, 89, 20], fluid pressure
$p_{h}\in\mathcal{P}_{h}^{\mathrm{DG}_{0}}\left(\mathcal{T}_{h}\right)$,
velocity approximation
$\bm{q}_{h}\in\mathcal{V}_{h}^{\mathrm{BDM}_{1}}\left(\mathcal{T}_{h}\right)$,
and concentration approximation
$c_{h}\in\mathcal{P}_{h}^{\mathrm{EG}_{1}}\left(\mathcal{T}_{h}\right)$.

We multiply the linear momentum balance equation (7) by a test function
$\bm{\psi_{u}}\in\mathcal{U}_{h}^{\mathrm{CG}_{2}}\left(\mathcal{T}_{h}\right)$.
The fully discretized linear momentum balance equation thus has the following
form:

$\mathcal{N}_{u}\left(\bm{\psi_{u}};\bm{u}_{h}^{n},p_{h}^{n}\right)=0,\quad\forall\bm{\psi_{u}}\in\mathcal{U}_{h}^{\mathrm{CG}_{2}}\left(\mathcal{T}_{h}\right),$
(44)

at each time step $t^{n}$, where

$\begin{split}\mathcal{N}_{u}\left(\bm{\psi_{u}};\bm{u}_{h}^{n},p_{h}^{n}\right)=&\sum_{T\in\mathcal{T}_{h}}\int_{T}\bm{\sigma}^{\prime}\left(\bm{u}_{h}^{n}\right):\nabla^{s}\bm{\psi_{u}}\>dV-\sum_{T\in\mathcal{T}_{h}}\int_{T}\alpha
p_{h}^{n}\mathbf{I}:\nabla^{s}\bm{\psi_{u}}\>dV\\\
&-\sum_{T\in\mathcal{T}_{h}}\int_{T}\bm{f}\bm{\psi_{u}}\>dV-\sum_{e\in\mathcal{E}_{h}^{N,u}}\int_{e}\bm{t_{D}}\bm{\psi_{u}}\>dS,\quad\forall\bm{\psi_{u}}\in\mathcal{U}_{h}^{\mathrm{CG}_{2}}\left(\mathcal{T}_{h}\right)\end{split}$

Here $\int_{T}\cdot\ dV$ and $\int_{e}\cdot\ dS$ refer to volume and surface
integrals, respectively, and $\nabla^{s}$ is the symmetric gradient operator.
Furthermore, the notation for
$\mathcal{N}_{u}\left(\bm{\psi_{u}};\bm{u}_{h}^{n},p_{h}^{n}\right)$ 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 $\psi_{p}\in\mathcal{P}_{h}^{\mathrm{DG}_{0}}\left(\mathcal{T}_{h}\right)$
and integrating by parts, resulting in:

$\mathcal{N}_{p}\left(\psi_{p};p_{h}^{n},\bm{q}_{h}^{n},c_{h}^{n}\right)=0,\quad\forall\psi_{p}\in\mathcal{P}_{h}^{\mathrm{DG}_{0}}\left(\mathcal{T}_{h}\right),$
(45)

for each time step $t^{n}$, where

$\begin{split}\mathcal{N}_{p}\left(\psi_{p};p_{h}^{n},\bm{q}_{h}^{n},c_{h}^{n}\right)&=\sum_{T\in\mathcal{T}_{h}}\int_{T}\left(\frac{1}{M}+\frac{\alpha^{2}}{K}\right)\mathrm{BDF}_{1}\left(p_{h}^{n}\right)\psi_{p}\>dV+\sum_{T\in\mathcal{T}_{h}}\int_{T}\nabla\cdot\left(\bm{q}_{h}^{n}\right)\psi_{p}\>dV\\\
&+\sum_{T\in\mathcal{T}_{h}}\int_{T}\frac{\alpha}{K}\mathrm{RK}_{1}(\sigma_{v})\psi_{p}\>dV+\sum_{T\in\mathcal{T}_{h}}\int_{T}\mathrm{RK}_{1}(\phi_{c})\psi_{p}\>dV-\sum_{T\in\mathcal{T}_{h}}\int_{T}g\psi_{p}\>dV.\end{split}$

For the Darcy velocity equation (10), we obtain

$\mathcal{N}_{v}\left(\bm{\psi}_{v};\bm{u}_{h}^{n},p_{h}^{n},\bm{q}_{h}^{n},c_{h}^{n}\right)=0,\quad\forall\bm{\psi}_{v}\in\mathcal{V}_{h}^{\mathrm{BDM}_{1}}\left(\mathcal{T}_{h}\right).$
(46)

where

$\begin{split}\mathcal{N}_{v}\left(\bm{\psi}_{v};\bm{u}_{h}^{n},p_{h}^{n},\bm{q}_{h}^{n},c_{h}^{n}\right)&:=\sum_{T\in\mathcal{T}_{h}}\int_{T}p_{h}^{n}\nabla\cdot\bm{\psi_{v}}\>dV\\\
&+\sum_{T\in\mathcal{T}_{h}}\int_{T}\bm{k}(\bm{u}_{h}^{n},c_{h}^{n})^{-1}\mu(c_{h}^{n})\bm{q}_{h}^{n}\bm{\psi_{v}}\>dV\\\
&+\sum_{e\in\mathcal{E}_{h}^{D,m}}\int_{e}p_{D}\bm{\psi}_{v}\cdot\mathbf{n}\>dS.\end{split}$

Lastly, for the advection-diffusion-reaction equations of species transport we
write:

$\mathcal{N}_{c}\left(\psi_{c};\bm{u}_{h}^{n},\bm{q}_{h}^{n},c_{h}^{n}\right)=0,\quad\forall\psi_{c}\in\mathcal{P}_{h}^{\mathrm{EG}_{1}}\left(\mathcal{T}_{h}\right)$
(47)

for each time step $t^{n}$, where

$\begin{split}\mathcal{N}_{c}\left(\psi_{c};\bm{u}_{h}^{n},\bm{q}_{h}^{n},c_{h}^{n}\right)&=\sum_{T\in\mathcal{T}_{h}}\int_{T}\phi\mathrm{BDF}_{4}\left(c_{h}^{n}\right)\psi_{c}\>dV+\sum_{T\in\mathcal{T}_{h}}\int_{T}{\bm{D}_{e}}^{*}(\phi^{n})\nabla
c_{h}^{n}\cdot\nabla\psi_{c}\>dV\\\
&-\sum_{e\in\mathcal{E}_{h}^{I}}\int_{e}\left\\{{\bm{D}_{e}}^{*}(\phi^{n})\nabla
c_{h}^{n}\right\\}_{\delta_{e}}\cdot\llbracket\psi_{c}\rrbracket\>dS\\\
&+\theta\sum_{e\in\mathcal{E}_{h}^{I}}\int_{e}\left\\{{\bm{D}_{e}}^{*}(\phi^{n})\nabla\psi_{c}\right\\}_{\delta_{e}}\cdot\llbracket
c_{h}^{n}\rrbracket\>dS\\\
&+\sum_{e\in\mathcal{E}_{h}^{I}}\int_{e}\frac{\beta}{h_{e}}{{\bm{D}_{e}}^{*}(\phi^{n})}_{{e}}\llbracket
c_{h}^{n}\rrbracket\cdot\llbracket\psi_{c}\rrbracket\>dS\\\
&-\sum_{T\in\mathcal{T}_{h}}\int_{T}\bm{q}_{h}^{n}c_{h}^{n}\cdot\nabla\psi_{c}\>dV+\sum_{e\in\mathcal{E}_{h}^{I}}\int_{e}\bm{q}_{h}^{n}\cdot\mathbf{n}c^{\mathrm{up}}_{h}\llbracket\psi_{c}\rrbracket\>dS\\\
&+\sum_{e\in\mathcal{E}_{h}^{\mathrm{Out}}}\int_{e}\bm{q}_{h}^{n}\cdot\mathbf{n}c_{h}^{n}\psi_{c}\>dS\\\
&-\sum_{T\in\mathcal{T}_{h}}\int_{T}R_{c}A_{s}\psi_{c}\>dV+\sum_{e\in\mathcal{E}_{h}^{\mathrm{In}}}\int_{e}\bm{q}_{h}^{n}\cdot\mathbf{n}c_{\mathrm{in}}\psi_{c}\>dS.\end{split}$

We note that the ${\bm{D}_{e}}^{*}(\phi^{n})$ is redefining
${\bm{D}_{e}}(\phi^{n})$ by including the numerical stabilization term, where

${\bm{D}_{e}}^{*}(\phi^{n}):=\bm{D}_{e}(\phi^{n})+\gamma
h\left\|\bm{q}_{h}^{n}\right\|\mathbf{I},$ (48)

as defined in [90, 91, 92]. The $\gamma
h\left\|\bm{q}_{h}^{n}\right\|\mathbf{I}$ term is often referred as the first
order artificial diffusivity coefficient [93, 94]. In our paper, we set the
tuning parameter $\gamma=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 [95, 90, 91, 92, 94, 96, 97, 98].

Also, $c^{\mathrm{up}}_{h}$ is an upwind value of $c_{h}^{n}$ defined as [99,
45]:

$c^{\mathrm{up}}_{h}=\left\\{\begin{array}[]{ll}c_{h}^{n+}&\text{ if
}\quad\bm{q}_{h}^{n}\cdot\mathbf{n}\geq 0\\\ c_{h}^{n-}&\text{ if
}\quad\bm{q}_{h}^{n}\cdot\mathbf{n}<0\end{array}\quad\forall e=\partial
T^{+}\cap\partial T^{-}\right.$ (49)

where $c_{h}^{n+}$ and $c_{h}^{n-}$ correspond to $c_{h}^{n}$ of $T^{+}$ and
$T^{-}$, respectively.

Lastly, the two parameters $\theta$ and $\beta$ define corresponding interior
penalty methods. The discretization becomes the symmetric interior penalty
Galerkin method (SIPG) when $\theta=-1$, the incomplete interior penalty
Galerkin method (IIPG) when $\theta=0$, and the non-symmetric interior penalty
Galerkin method (NIPG) when $\theta=1$ [45]. In this study, we set $\theta=-1$
for the simplicity and $\beta=1.1$ throughout this paper.

######  Remark 1.

For the momentum balance equation (7), the traction boundary condition
$\bm{t_{D}}$ (traction) is applied weakly on each $e\in\mathcal{E}_{h}^{N,u}$
in (44), while the displacement boundary condition $\bm{u}_{D}$ is strongly
enforced on each $e\in\mathcal{E}_{h}^{D,u}$. For the mass balance equation
(12), since we use a mixed formulation, the flux boundary condition ${q_{D}}$
is strongly applied on each $e\in\mathcal{E}_{h}^{N,m}$, but the pressure
boundary condition $p_{D}$ is weakly applied on each $\mathcal{E}_{h}^{D,m}$
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 $\mathrm{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 $\mathrm{BDF}_{4}$ to
capture a sharp front in the advection dominated regime [45].

###  3.4 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.

Algorithm 1 Splitting algorithm for hydro-mechanical-chemical coupling model

1:Initialize all input parameters $\triangleright$ $p_{0}$ and $c_{0}$ must be
provided.

2:Solve the equilibrium state for $\bm{u_{0}}$ $\triangleright$ see (44)

3:Update $\phi^{0}$, $\bm{k}^{0}$, $\bm{D}_{e}^{0}$, and $A_{s}^{0}$
$\triangleright$ see (23), (26), (15), (25)

4:for each time step $t^{n}$ do

5:  _Part 1: coupling solid and fluid mechanics_

6: Set $\iota\rightarrow 0$ as the nonlinear iterations counter

7: ${p}_{h}^{n-1}\rightarrow p_{h}^{n,\iota=0}$,
$\bm{q}_{h}^{n-1}\rightarrow\bm{q}_{h}^{n,\iota=0}$,
$\bm{u}_{h}^{n-1}\rightarrow\bm{u}_{h}^{n,\iota=0}$

8: for each fixed stress iteration step $\left(\cdot\right)^{n,\iota}$ until
$\delta\phi^{n,\iota}$ $<$ $\mathrm{TOL}$ do

9: Solve (45) and (46) w.r.t. $p_{h}^{n}$ and $\bm{q}_{h}^{n}$ for fixed
$\bm{u}_{h}^{n}:=\bm{u}_{h}^{n,\iota-1},c_{h}^{n}:=\hat{c}_{h}^{n}$ to get
$p_{h}^{n,\iota},\bm{q}_{h}^{n,\iota}$

10: Calculate $\phi_{f}^{n,\iota}$ $\triangleright$ see (51)

11: Solve (44) w.r.t. $\bm{u}_{h}^{n}$ and for fixed
$p_{h}^{n}:=p_{h}^{n,\iota}$ to get $\bm{u}_{h}^{n,\iota}$

12: Calculate $\phi_{m}^{n,\iota}$ $\triangleright$ see (23)

13: Evaluate $\mathrm{F}\left(\dot{\sigma_{v}}^{n,\iota}\right)$
$\triangleright$ see (52)

14: Evaluate $\delta\phi^{n,\iota}$ $\triangleright$ see (50)

15: Update $\bm{k}^{n,\iota}$ $\triangleright$ see (23), (26)

16: end for

17: ${p}_{h}^{n,\iota}\rightarrow{p}_{h}^{n}$,
$\bm{q}_{h}^{n,\iota}\rightarrow\bm{q}_{h}^{n}$,
$\bm{u}_{h}^{n,\iota}\rightarrow\bm{u}_{h}^{n}$

18:  _Part 2: chemical process_

19: Update $\phi^{n}$, $\bm{D}_{e}^{n}$, and $A_{s}^{n}$ $\triangleright$ see
(23), (15), (25)

20: Solve (47) w.r.t. ${c}_{h}^{n}$ for fixed
$\bm{u}_{h}^{n}:=\bm{u}_{h}^{n},\bm{q}_{h}^{n}:=\bm{q}_{h}^{n}$ to get
${c}_{h}^{n}$

21: Extrapolate $\hat{c}_{h}^{n+1}$ $\triangleright$ see (53)

22: Calculate $\Delta t^{n+1}$ $\triangleright$ see (34)

23: Evaluate $\mathrm{F}\left(\hat{q}^{n+1}\right)$ and
$\mathrm{F}\left(\hat{\dot{\phi_{c}}}^{n+1}\right)$ $\triangleright$ see (55),
(54)

24: Update $\hat{\phi}^{n+1}$, $\hat{\bm{k}}^{n+1}$, $\hat{\bm{D}}_{t}^{n+1}$,
$\hat{\mu}^{n+1}$, and $\hat{A_{s}}^{n+1}$ $\triangleright$ see (53), (24),
(26), (15), (27), (25)

25: $p_{h}^{n}\rightarrow p_{h}^{n-1}$,
$\bm{q}_{h}^{n}\rightarrow\bm{q}_{h}^{n-1}$,
$\bm{u}_{h}^{n}\rightarrow\bm{u}_{h}^{n-1}$, $c_{h}^{n}\rightarrow
c_{h}^{n-1}$ $\triangleright$ update time step $n-1$

26: Output

27:end for

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[19, 18, 73, 68]. At each iteration $\iota$ we solve
(45) and (46) for the velocity $\bm{q}_{h}^{n,\iota}$ and the pressure
${p}_{h}^{n,\iota}$ using a monolithic method (line 9) based on given
displacement $\bm{u}_{h}^{n,\iota-1}$ from previous nonlinear iteration and
concentration extrapolated $\hat{c}_{h}^{n}$ from previous time step. Then, we
couple with (44) using the fixed-stress split scheme based on the pressure
${p}_{h}^{n,\iota}$ computed at the current nonlinear iteration (line 11). The
convergence criterion is based on $\delta\phi^{n,\iota}$ (Algorithm 1 line 8),
which is defined as:

$\delta\phi^{n,\iota}:=\frac{\phi_{m}^{n,\iota}-\phi_{f}^{n,\iota}}{\phi_{m}^{n,\iota}}.$
(50)

Here, $\phi_{m}^{n,\iota}$ is the porosity resulting from the solid
deformation (23) and $\phi_{f}^{n,\iota}$ is the porosity resulting from the
fluid flow problem defined as [68, 18, 73]:

$\phi_{f}^{n,\iota}=\phi^{n-1}+\frac{\left(\alpha-\phi^{n-1}\right)}{K}\left(p^{n,\iota}-p^{n-1}\right),$
(51)

where $\left(\cdot\right)^{\iota}$ represents iteration counter inside the
fixed-stress loop. From the fixed stress split concept (51) is the $\phi$
predictor, while (23) is the $\phi$ corrector [19, 18, 73, 68]. Hence, when
$\phi_{m}^{n,\iota}$ and $\phi_{f}^{n,\iota}$ converge, i.e.,
$\delta\phi^{n,\iota}$ $<$ $\mathrm{TOL}$, the fixed-stress loop is completed.
The tolerance $\mathrm{TOL}$ is set as $1\times 10^{-6}$ throughout this
study. Note that the flow equations, (45) and (46), are solved by assuming
that $\frac{\partial\sigma_{v}}{\partial t}=0$, i.e.,
$\mathrm{F}\left(\dot{\sigma_{v}}^{n,\iota}\right)$ 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
$\mathrm{F}\left(\dot{\sigma_{v}}^{n,\iota}\right)$ is defined as:

$\mathrm{F}\left(\dot{\sigma_{v}}^{n,\iota}\right):=\sum_{T\in\mathcal{T}_{h}}\int_{T}\frac{\alpha}{K}\mathrm{RK}_{1}\left(\sigma_{v}(\bm{u}^{n,\iota},\bm{u}^{n-1},{p}^{n,\iota},{p}^{n-1})\right)\psi_{p}\>dV.$
(52)

The second part (from line 18) focuses on solving advection-diffusion-reaction
equation (47), using $\bm{q}_{h}^{n}$, $\phi^{n}$, $\bm{D}_{e}^{n}$, and
$A_{s}^{n}$ obtained from the first part. One could view this strategy as a
one-way coupling scheme between coupled hydro-mechanical and advection-
diffusion-reaction equations. Next, Algorithm 1 line 21 linearizes
${c}_{h}^{n+1}$ by extrapolating ${c}_{h}^{n}$ and ${c}_{h}^{n-1}$ to
$\hat{c}_{h}^{n+1}$ by using (38):

$\hat{c}_{h}^{n+1}=\mathrm{EX}\left({c}_{h}\right)$ (53)

where $\left(\hat{\cdot}\right)^{n}$ represents an extrapolation value based
on the extrapolation described in (38). Subsequently, we evaluate
$\mathrm{F}\left(\dot{\phi_{c}}^{n,\iota}\right)$ and
$\mathrm{F}\left(q^{n,\iota}\right)$, which are defined as

$\mathrm{F}\left(\hat{\dot{\phi_{c}}}^{n+1}\right):=\sum_{T\in\mathcal{T}_{h}}\int_{T}\mathrm{RK}_{1}\left(\phi_{c}(\hat{c}^{n+1},{c}^{n})\right)\psi_{p}\>dV,$
(54)

and

$\mathrm{F}\left(\hat{q}^{n+1}\right):=\sum_{T\in\mathcal{T}_{h}}\int_{T}R_{c}(\hat{c}^{n+1},p^{n})A_{s}(\hat{c}^{n+1},\bm{u}^{n})\psi_{c}\>dV,$
(55)

using $\hat{c}_{h}^{n+1}$ calculated by (53). We note that the equation (47)
becomes linear by employing $\hat{c}_{h}^{n+1}$ to calculate
$\mathrm{F}\left(q\right)$. Also, the porosity alteration as a result of
calcite dissolution/precipitation (Algorithm 1 line 24) is computed by

$\displaystyle\hat{\phi}^{n+1}=\hat{\phi}_{c}^{n+1}$
$\displaystyle=\mathrm{RK}_{4}\left(\Gamma\left(\bm{u}^{n},\hat{c}_{h}^{n+1}\right)\right).$
(56)

Note that the porosity change due to the calcite dissolution/precipitation
reactions is additional to the porosity change by solid deformation, (23).
Subsequently, $\hat{\bm{k}}^{n+1}$, $\hat{\bm{D}}_{t}^{n+1}$, and
$\hat{A_{s}}^{n+1}$ are determined using $\hat{\phi}^{n+1}$. Lastly, we also
calculate $\hat{\mu}^{n+1}$ using $\hat{c}_{h}^{n+1}$, 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 $\mathrm{XeonE5\\_2650v4}$ with a 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 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.

##  4 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.

Figure 1: Geometry and notation used to define material properties;
($\bm{\mathrm{a}}$) example 1: single-layer porous medium ($\Omega_{500}$),
($\bm{\mathrm{b}}$) example 2: three-layer porous medium ($\Omega_{500}$ and
$\Omega_{501}$), ($\bm{\mathrm{c}}$) example 3: heterogeneous porous medium
(the arithmetic mean of $\bm{k}=8.8\times 10^{-10}\mathbf{I}$ $\mathrm{m}^{2}$
with correlation length in $x$\- and $y$-direction of 5 and 1 m,
respectively.), and ($\bm{\mathrm{d}}$) mesh used for all examples (number of
element is 7852).

###  4.1 Example 1

In the first example, the computational domain is defined as $\Omega_{500}$ =
$[0,100]\times[0,30]$, which presents a single layer as shown in Figure 1a.
Following the typical physical properties of rocks [65], we set $K=8.4$
$\mathrm{G}\mathrm{P}\mathrm{a}$, $\alpha=0.74$, $\nu=0.18$, $\phi=0.2$,
$\bm{k}=8.8\times 10^{-10}\mathbf{I}$ $\mathrm{m}^{2}$. The fluid properties
considered in this case are $c_{f}=1.0\times 10^{-10}$
$\mathrm{P}\mathrm{a}^{-1}$, $\rho=1000$
$\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$, $\bm{D}=1.0\times 10^{-12}$
$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$, and $\mu$ is calculated using (27) by
setting $\mu_{h}=5.0$ and $\mu_{l}=1.0\times 10^{-4}$
$\mathrm{P}\mathrm{a}\mathrm{/}\mathrm{s}$ corresponding to $c_{h}=1.68$ and
$c_{l}=0.0$
$\mathrm{m}\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{/}\mathrm{m}^{3}$,
respectively. Next, the boundary conditions for all these examples are applied
as follows. For the momentum balance equation (7), we assume
$\bm{u}_{D}\cdot\mathbf{n}=0$ on $\partial\Omega_{1}$, $\partial\Omega_{3}$,
and $\partial\Omega_{4}$. Furthermore, $\bm{t}_{D}=\left[0.0,-2.0\times
10^{6}\right]$ $\mathrm{P}\mathrm{a}$ is applied on $\partial\Omega_{2}$.
Therefore, the medium is under compression. For the mass balance equation
(12), the boundary condition $q_{D}=0$ is set on $\partial\Omega_{2}$ and
$\partial\Omega_{4}$ and we impose $p_{D}=1\times 10^{5}$
$\mathrm{P}\mathrm{a}$ on $\partial\Omega_{3}$. Here, for the mass balance
equation (12), we test two different scenarios on $\partial\Omega_{1}$, where
scenario (a) corresponds to $q_{D}=2\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ and scenario (b) is characterized by
$q_{D}=1\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{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

$\mathrm{I.V.}=q_{D}t^{n}A_{\mathrm{d}},$ (57)

where $A_{\mathrm{d}}=30$\mathrm{m}^{2}$$ is the surface area of
$\partial\Omega_{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 $\partial\Omega_{1}$. The initial
pressure $p_{0}$ is $1\times 10^{6}$ $\mathrm{P}\mathrm{a}$, the initial
concentration $c_{0}$ is calculated by (22) using $p=p_{0}$ and $\tau=20$
$\mathrm{C}$, and the initial displacement $\bm{u}_{0}$ is calculated as
stated in Algorithm 1. The penalty parameter ($\beta$) is set to be 1.1 for
the EG method. The $\mathrm{CFL}$ is used as 0.1 for calculating $\Delta
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.

Figure 2: Example 1: concentration fields, $c$; $\mathrm{I.V.}=42$
$\mathrm{m}^{3}$ using ($\bm{\mathrm{a}}$) $q_{D}=2\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ at $\partial\Omega_{1}$ and
($\bm{\mathrm{b}}$) $q_{D}=1\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$
at $\partial\Omega_{1}$; and $\mathrm{I.V.}=180$ $\mathrm{m}^{3}$ using
($\bm{\mathrm{c}}$) $q_{D}=2\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$
at $\partial\Omega_{1}$ and ($\bm{\mathrm{d}}$) $q_{D}=1\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ at $\partial\Omega_{1}$. The boundary
conditions shown in this picture are corresponding to $\partial\Omega_{1}$ and
$\partial\Omega_{3}$ of the mass balance equation (12).

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 $\bm{u}$, see Figure 3b and g. Figure 3a and f
illustrate the reduction of $\phi$ 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 in which
$\frac{\partial\phi_{m}}{\partial t}$ is positive in the left part of the
domain and negative in the right part of the domain.

The $\frac{\partial\phi_{c}}{\partial 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 $\frac{\partial\phi_{c}}{\partial t}$ is
positive where the dissolution occurs and negative where the perception
occurs. At this time step, the maximum magnitude of
$\frac{\partial\phi_{c}}{\partial t}$ is $10^{-7}$, which is much less
compared to that of $\frac{\partial\phi_{m}}{\partial 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 $\mu$ 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
$\frac{\partial\phi_{m}}{\partial t}$, $\frac{\partial\phi_{c}}{\partial t}$,
and $\mu$ alteration can be seen in $\bm{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.

Figure 3: Example 1: using $q_{D}=2\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ at $\partial\Omega_{1}$; $\mathrm{I.V.}=42$
$\mathrm{m}^{3}$: ($\bm{\mathrm{a}}$) the rate of change of porosity according
to mechanics deformation, $\frac{\partial\phi_{m}}{\partial t}$, see (23),
($\bm{\mathrm{b}}$) the fluid pressure, $p$, in surface and the displacement,
$\bm{u}$, in grey arrows, ($\bm{\mathrm{c}}$) the rate of change of porosity
according to calcite dissolution/precipitation,
$\frac{\partial\phi_{c}}{\partial t}$, see (24), ($\bm{\mathrm{d}}$) the fluid
viscosity, $\mu$, and ($\bm{\mathrm{e}}$) the magnitude of the fluid velocity,
$\left\|\bm{q}\right\|$. $\mathrm{I.V.}=480$ $\mathrm{m}^{3}$:
($\bm{\mathrm{f}}$) $\frac{\partial\phi_{m}}{\partial t}$, ($\bm{\mathrm{g}}$)
$p$ in surface and $\bm{u}$ in grey arrows, ($\bm{\mathrm{h}}$)
$\frac{\partial\phi_{c}}{\partial t}$, ($\bm{\mathrm{i}}$) $\mu$, and
($\bm{\mathrm{j}}$) $\left\|\bm{q}\right\|$. Note that the magnitude of
$\bm{u}$ is from $0.0$ to $3.0\times 10^{-2}$, and the trend of the results of
the $q_{D}=1\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$ at
$\partial\Omega_{1}$ case are similar.

###  4.2 Example 2

In the second example, we consider three layers
($\Omega_{500}$=$[0,100]\times[10,20]$, $\Omega_{501}$=$[0,100]\times[20,30]$,
and $\Omega_{502}$=$[0,100]\times[0,10]$) as the computational domain. See
Figure 1b. In $\Omega_{500}$, we set $\bm{k}=8.8\times 10^{-10}\mathbf{I}$
$\mathrm{m}^{2}$, while $\bm{k}=8.8\times 10^{-11}\mathbf{I}$ $\mathrm{m}^{2}$
in $\Omega_{501}$ and $\Omega_{502}$. Thus, in this case, the top
$\Omega_{501}$ and bottom $\Omega_{502}$ layers have one order of magnitude of
$\bm{k}$ less than that of the middle layer $\Omega_{500}$. 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 porous medium, even though the concentration
fields at the early time are similar between the high, $q_{D}=2\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$, and low, $q_{D}=1\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{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, $\Omega_{501}$ and $\Omega_{502}$, have lower
permeability than the middle layer, $\Omega_{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.

Figure 4: Example 2: concentration fields, $c$; $\mathrm{I.V.}=42$
$\mathrm{m}^{3}$ using ($\bm{\mathrm{a}}$) $q_{D}=2\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ at $\partial\Omega_{1}$ and
($\bm{\mathrm{b}}$) $q_{D}=1\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$
at $\partial\Omega_{1}$; and $\mathrm{I.V.}=180$ $\mathrm{m}^{3}$ using
($\bm{\mathrm{c}}$) $q_{D}=2\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$
at $\partial\Omega_{1}$ and ($\bm{\mathrm{d}}$) $q_{D}=1\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ at $\partial\Omega_{1}$. The boundary
conditions shown in this picture are corresponding to $\partial\Omega_{1}$ and
$\partial\Omega_{3}$ of the mass balance equation (12).

In Figure 5, the behavior of the concentration and velocity fields, together
with temporal porosity alteration ($\frac{\partial\phi_{c}}{\partial t}$), are
illustrated for both injection scenarios. As mentioned earlier, due to the
difference of viscosity ($\mu$) 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 $\partial\Omega_{3}$, see
Figures 5b-d. One could see that the reaction front shown by
$\frac{\partial\phi_{c}}{\partial t}$ progresses similarly to the
concentration front shown by the black contours. Besides, as the concentration
field develops, the change in $\mu$ 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
$\frac{\partial\phi_{c}}{\partial 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.

Figure 5: Example 2: the results of the rate of change of porosity according
to calcite dissolution/precipitation ($\frac{\partial\phi_{c}}{\partial t}$)
shown in surface plot, concentration ($c$), shown in black contour (10
contours ranging from 0.12 to 1.6), and the fluid velocity ($\bm{q}$) shown in
arrows with $q_{D}=2\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$ at
$\partial\Omega_{1}$; ($\bm{\mathrm{a}}$) $\mathrm{I.V.}=36$ $\mathrm{m}^{3}$
($t=6000$ $\mathrm{s}$), ($\bm{\mathrm{b}}$) $\mathrm{I.V.}=72$
$\mathrm{m}^{3}$ ($t=12000$ $\mathrm{s}$), ($\bm{\mathrm{c}}$)
$\mathrm{I.V.}=156$ $\mathrm{m}^{3}$ ($t=26000$ $\mathrm{s}$), and
($\bm{\mathrm{d}}$) $\mathrm{I.V.}=360$ $\mathrm{m}^{3}$ ($t=60000$
$\mathrm{s}$), and with $q_{D}=1\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ at $\partial\Omega_{1}$; ($\bm{\mathrm{e}}$)
$\mathrm{I.V.}=36$ $\mathrm{m}^{3}$ ($t=1200$ $\mathrm{s}$),
($\bm{\mathrm{f}}$) $\mathrm{I.V.}=72$ $\mathrm{m}^{3}$ ($t=24000$
$\mathrm{s}$), ($\bm{\mathrm{g}}$) $\mathrm{I.V.}=156$ $\mathrm{m}^{3}$
($t=32000$ $\mathrm{s}$), and ($\bm{\mathrm{h}}$) $\mathrm{I.V.}=360$
$\mathrm{m}^{3}$ ($t=120000$ $\mathrm{s}$).

###  4.3 Example 3

In the given computational domain $\Omega_{500}$ = $[0,100]\times[0,30]$, we
investigate the setup with the heterogeneous $\bm{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
$\bm{k}=1\times 10^{-10}\mathbf{I}$ $\mathrm{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, the concentration fields are
illustrated for two different injection rates at different injected fluid
volumes ($\mathrm{I.V.}=42$ $\mathrm{m}^{3}$ and $\mathrm{I.V.}=180$
$\mathrm{m}^{3}$). Unlike the two previous examples, the preferential paths
are established not only because of the flow instability resulting from the
$\mu$ difference but also due to the high $\bm{k}$ channels inherited from the
nature of heterogeneous porous media. During the early time, the concentration
field of the high, $q_{D}=2\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$,
and low, $q_{D}=1\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$, injection
rate cases are similar, see Figure 6a-b. The results of the concentration with
the effects from the reaction are different at a later time (see Figure 6c-d).
During the early time for both cases, the developed fingers follow the high
permeable paths. At a later time, however, the results of the two scenarios
are very different. For the high injection rate case, the top finger continues
developing while the middle 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.

Figure 6: Example 3: concentration fields, $c$; $\mathrm{I.V.}=42$
$\mathrm{m}^{3}$ using ($\bm{\mathrm{a}}$) $q_{D}=2\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ at $\partial\Omega_{1}$ and
($\bm{\mathrm{b}}$) $q_{D}=1\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$
at $\partial\Omega_{1}$; and $\mathrm{I.V.}=180$ $\mathrm{m}^{3}$ using
($\bm{\mathrm{c}}$) $q_{D}=2\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$
at $\partial\Omega_{1}$ and ($\bm{\mathrm{d}}$) $q_{D}=1\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ at $\partial\Omega_{1}$. The boundary
conditions shown in this picture are corresponding to $\partial\Omega_{1}$ and
$\partial\Omega_{3}$ of the mass balance equation (12).

Figure 7 provides further insight into the reactive flow dynamics. It shows
for both injection scenarios how the reaction fronts and flow fields evolve in
time. As mentioned earlier, all the initial fingers at the beginning vanish
except one that reaches the outlet $\partial\Omega_{3}$. The flow velocity
field variations in time depict the emergence of the dominant finger. 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 $\frac{\partial\phi_{c}}{\partial t}$ illustrate the applicability of
the presented coupled model for heterogeneous porous media.

Figure 7: Example 3: the results of the rate of change of porosity according
to calcite dissolution/precipitation ($\frac{\partial\phi_{c}}{\partial t}$)
shown in surface plot, concentration ($c$), shown in black contour (10
contours ranging from 0.12 to 1.6), and the fluid velocity ($\bm{q}$) shown in
arrows with $q_{D}=2\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$ at
$\partial\Omega_{1}$; ($\bm{\mathrm{a}}$) $\mathrm{I.V.}=36$ $\mathrm{m}^{3}$
($t=6000$ $\mathrm{s}$), ($\bm{\mathrm{b}}$) $\mathrm{I.V.}=72$
$\mathrm{m}^{3}$ ($t=12000$ $\mathrm{s}$), ($\bm{\mathrm{c}}$)
$\mathrm{I.V.}=156$ $\mathrm{m}^{3}$ ($t=26000$ $\mathrm{s}$), and
($\bm{\mathrm{d}}$) $\mathrm{I.V.}=360$ $\mathrm{m}^{3}$ ($t=60000$
$\mathrm{s}$), and with $q_{D}=1\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ at $\partial\Omega_{1}$; ($\bm{\mathrm{e}}$)
$\mathrm{I.V.}=36$ $\mathrm{m}^{3}$ ($t=1200$ $\mathrm{s}$),
($\bm{\mathrm{f}}$) $\mathrm{I.V.}=72$ $\mathrm{m}^{3}$ ($t=24000$
$\mathrm{s}$), ($\bm{\mathrm{g}}$) $\mathrm{I.V.}=156$ $\mathrm{m}^{3}$
($t=32000$ $\mathrm{s}$), and ($\bm{\mathrm{h}}$) $\mathrm{I.V.}=360$
$\mathrm{m}^{3}$ ($t=120000$ $\mathrm{s}$).

Next, we investigate the local mass conservation property of the proposed
framework in the heterogeneous domain. The local mass conservation of each
cell at each time step, $\mathrm{r^{n}_{mass}}$, is calculated by

$\mathrm{r^{n}}_{\mathrm{mass}}:=\int_{T}\left(\frac{1}{M}+\frac{\alpha^{2}}{K}\right)\frac{p^{n}-p^{n-1}}{\Delta
t^{n}}+\frac{\alpha}{K}\frac{\sigma_{v}^{n}-\sigma_{v}^{n-1}}{\Delta
t^{n}}+\frac{\hat{\phi_{c}}^{n}-\phi_{c}^{n-1}}{\Delta
t^{n}}dV+\sum_{e\in\mathcal{E}_{h}}\int_{e}\bar{\bm{q}}^{n}\cdot\left.\mathbf{n}\right|_{\mathrm{e}}dS,$
(58)

and the discrete numerical flux approximated by BDM,
$\bar{\bm{q}}^{n}\cdot\left.\mathbf{n}\right|_{\mathrm{e}}$, is defined by

$\displaystyle\bar{\bm{q}}^{n}$ $\displaystyle:=\bm{q}_{h}^{n}\quad\forall
T\in\mathcal{T}_{h},$ (59)
$\displaystyle\bar{\bm{q}}^{n}\cdot\left.\mathbf{n}\right|_{\mathrm{e}}$
$\displaystyle:=-q_{D}\quad\forall e\in\mathcal{E}_{h}^{N,m},$ (60)
$\displaystyle\bar{\bm{q}}^{n}\cdot\left.\mathbf{n}\right|_{\mathrm{e}}$
$\displaystyle:=-\bm{q}_{h}^{n}\cdot\mathbf{n}\quad\forall
e\in\mathcal{E}_{h}^{D,m}.$ (61)

In Figure 8, the values of $\mathrm{r^{n}_{mass}}$ are illustrated for each
case and time. One could see that the magnitude of $\mathrm{r^{n}_{mass}}$ is
always less than $1\times 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 $\mathrm{r^{n}_{mass}}$ than that of the low
injection rate case.

Figure 8: Example 3: the illustration of the local mass conservative property
with $q_{D}=2\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$ at
$\partial\Omega_{1}$; ($\bm{\mathrm{a}}$) $\mathrm{I.V.}=36$ $\mathrm{m}^{3}$
($t=6000$ $\mathrm{s}$), ($\bm{\mathrm{b}}$) $\mathrm{I.V.}=72$
$\mathrm{m}^{3}$ ($t=12000$ $\mathrm{s}$), ($\bm{\mathrm{c}}$)
$\mathrm{I.V.}=156$ $\mathrm{m}^{3}$ ($t=26000$ $\mathrm{s}$), and
($\bm{\mathrm{d}}$) $\mathrm{I.V.}=360$ $\mathrm{m}^{3}$ ($t=60000$
$\mathrm{s}$), and with $q_{D}=1\times 10^{-4}$
$\mathrm{m}\mathrm{/}\mathrm{s}$ at $\partial\Omega_{1}$; ($\bm{\mathrm{e}}$)
$\mathrm{I.V.}=36$ $\mathrm{m}^{3}$ ($t=1200$ $\mathrm{s}$),
($\bm{\mathrm{f}}$) $\mathrm{I.V.}=72$ $\mathrm{m}^{3}$ ($t=24000$
$\mathrm{s}$), ($\bm{\mathrm{g}}$) $\mathrm{I.V.}=156$ $\mathrm{m}^{3}$
($t=32000$ $\mathrm{s}$), and ($\bm{\mathrm{h}}$) $\mathrm{I.V.}=360$
$\mathrm{m}^{3}$ ($t=120000$ $\mathrm{s}$).

###  4.4 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 $\Omega_{500}$ = $[0,100]\times[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:

$\bm{k}:=\left[\begin{array}[]{cc}k_{xx}&0.0\\\
0.0&0.1k_{xx}\end{array}\right],$ (62)

where $k_{xx}=8.8\times 10^{-10}$ $\mathrm{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\times 10^{-5}$, which is the tolerance set for the
fixed-stress loop.

Figure 9: Example 4: the results of the rate of change of porosity according
to calcite dissolution/precipitation ($\frac{\partial\phi_{c}}{\partial t}$)
shown in surface plot, concentration ($c$), shown in black contour (10
contours ranging from 0.12 to 1.6), and the fluid velocity ($\bm{q}$) shown in
arrows with $q_{D}=2\times 10^{-4}$ $\mathrm{m}\mathrm{/}\mathrm{s}$ at
$\partial\Omega_{1}$; ($\bm{\mathrm{a}}$) $\mathrm{I.V.}=180$ $\mathrm{m}^{3}$
($t=30000$ $\mathrm{s}$), ($\bm{\mathrm{b}}$) $\mathrm{I.V.}=480$
$\mathrm{m}^{3}$ ($t=80000$ $\mathrm{s}$), and ($\bm{\mathrm{c}}$)
$\mathrm{I.V.}=720$ $\mathrm{m}^{3}$ ($t=120000$ $\mathrm{s}$), and the local
mass conservative property; ($\bm{\mathrm{d}}$) $\mathrm{I.V.}=180$
$\mathrm{m}^{3}$ ($t=30000$ $\mathrm{s}$), ($\bm{\mathrm{e}}$)
$\mathrm{I.V.}=480$ $\mathrm{m}^{3}$ ($t=80000$ $\mathrm{s}$), and
($\bm{\mathrm{f}}$) $\mathrm{I.V.}=720$ $\mathrm{m}^{3}$ ($t=120000$
$\mathrm{s}$).

###  4.5 Discussion

The main observations of the foregoing numerical examples can be summarized as
follows:

  1. 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. 2.

Using the applied set of the input parameters resulted in a more noticeable
mechanical effect on the change in $\phi$ (and subsequently in $\bm{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 $\mu$ resulted from the
change in $c$ is significant, resulting in the development of preferential
flow paths.

  3. 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 fixed-stress 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\times 10^{-5}$ second per degrees of freedom per each time step. All
simulations were computed on 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 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 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 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. TK also thanks the 2019 Computers & Geosciences Research grant for
the additional support. SL is supported by the National Science Foundation
under Grant No. NSF DMS-1913016. FB thanks Horizon 2020 Program for Grant
H2020 ERC CoG 2015 AROMA-CFD project 681447 that supported the development of
multiphenics. JC acknowledges support from the Research Grants Council of Hong
Kong (Project 27205918).

##  7 CRediT authorship contribution statement

T. 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, Formal 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 main dependencies are Numpy ($\geq$ 1.16.5), FEniCS ($\geq$
2018.1.0) with PETSc ($\geq$ 3.10.5) and petsc4py ($\geq$ 3.10), and
multiphenics ($\geq$ 0.2.0).

## References

  * [1] H. Nick, A. Raoof, F. Centler, M. Thullner, P. Regnier, Reactive dispersive contaminant transport in coastal aquifers: numerical simulation of a reactive henry problem, Journal of contaminant hydrology 145 (2013) 90–104.
  * [2] M. Hu, T. Hueckel, Environmentally enhanced crack propagation in a chemically degrading isotropic shale, Géotechnique 63 (4) (2013) 313–321.
  * [3] S. Pandey, A. Chaudhuri, S. Kelkar, V. Sandeep, H. Rajaram, Investigation of permeability alteration of 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 evolution in a carbonate reservoir: A thermo-hydro-chemical study, Geothermics 69 (2017) 45–54.
  * [5] H. M. Nick, K.-H. Wolf, D. Brhun, Mixed CO2–water injection into geothermal reservoirs: A numerical study, 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-mechanics and 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 during viscous 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-reservoir simulator, 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 productivity in fractured porous media using full factorial experimental design, Journal of Petroleum Science and Engineering 181 (2019) 106233\.
  * [12] T. Kadeethum, S. Salimzadeh, H. Nick, Well productivity evaluation in deformable single-fracture media, Geothermics 87 (2020).
  * [13] M. Nejati, M. Dambly, M. Saar, A methodology to determine the elastic properties of anisotropic rocks from a single uniaxial compression test, Journal of Rock Mechanics and Geotechnical Engineering 11 (6) (2019) 1166–1183.
  * [14] S. Salimzadeh, E. Hagerup, T. Kadeethum, H. Nick, The effect of stress distribution on the shape and direction 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, Computational Geosciences (2020) 1–18.
  * [17] C. Zhang, S. Zarrouk, R. Archer, A mixed finite element solver for natural convection in porous media using 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 coupled flow and deformation in heterogeneous poroelastic media, Computer Methods in Applied Mechanics and Engineering 341 (2018) 788–806.
  * [19] J. Kim, H. Tchelepi, R. Juanes, Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits, Computer Methods in Applied Mechanics and Engineering 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 geomechanics, Computational Geosciences 15 (4) (2011) 647\.
  * [22] H. Nick, S. Matthai, A hybrid finite-element finite-volume method with embedded discontinuities for solute 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 volume finite element method for multi-phase flow in heterogeneous porous media, Journal of Computational Physics 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 porous media, Mathematical Geosciences (2020) 1–39.
  * [26] T. Kadeethum, T. Jørgensen, H. Nick, Physics-informed neural networks for solving nonlinear diffusivity 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 Problems of Nonlinear Biot’s Equations: Batch Training, in: 54th US Rock Mechanics/Geomechanics Symposium, 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 engineered 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, Geothermal 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 of abstracts (2012).
  * [32] B. Aagaard, C. Williams, M. Knepley, PyLith: A finite-element code for modeling quasi-static and dynamic 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 simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media, Environmental Earth 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, Reservoir modelling using parametric surfaces and dynamically adaptive fully unstructured grids, in: ECMOR XV-15th European Conference on the Mathematics of Oil Recovery, European Association of Geoscientists & Engineers, 2016, p. 494.
  * [36] A. Obeysekara, Q. Lei, P. Salinas, D. Pavlidis, J. Latham, J. Xiang, C. Pain, et al., A fluid-solid coupled approach for numerical modeling of near-wellbore hydraulic fracturing and flow dynamics with adaptive mesh refinement, in: 50th US Rock Mechanics/Geomechanics Symposium, American Rock 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 method with mesh adaptivity, Computers and Geotechnics 103 (2018) 229–241.
  * [38] J. H. Prévost, Dynaflow, Princeton University, Princeton, NJ 8544 (1983).
  * [39] M. Tene, M. Al Kobaisi, H. Hajibeygi, Algebraic multiscale method for flow in heterogeneous porous media 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 Computational Physics 299 (2015) 472–486.
  * [41] M. HosseiniMehr, C. Vuik, H. Hajibeygi, Adaptive dynamic multilevel simulation of fractured geothermal 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 complex reservoirs, 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 enhanced geothermal systems, Geothermics 81 (2019) 88–100.
  * [44] E. Keilegavlen, R. Berge, A. Fumagalli, M. Starnoni, I. Stefansson, J. Varela, I. Berre, Porepy: An open-source software for simulation of multiphysics processes in fractured porous media, arXiv preprint arXiv:1908.09869 (2019).
  * [45] B. Riviere, Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation, SIAM, 2008.
  * [46] S. Lee, Y. Lee, M. Wheeler, A locally conservative enriched Galerkin approximation and efficient solver for 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 Mathematik 112 (1) (2009) 115–152.
  * [48] J. Choo, Large deformation poromechanics with local mass conservation: An enriched Galerkin finite element 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 enrichment of 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 capillary pressure, Journal of Computational Physics 367 (2018) 65–86.
  * [51] J. Choo, S. Lee, Enriched Galerkin finite elements for coupled poromechanics with local mass conservation, Computer Methods in Applied Mechanics and Engineering 341 (2018) 311–332.
  * [52] J. Choo, Stabilized mixed continuous/enriched Galerkin formulations for locally mass conservative poromechanics, 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 for poroelasticity 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 for poroelasticity 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 low compressible porous media, International Journal for Numerical and Analytical Methods in Geomechanics 36 (12) (2012) 1507–1522.
  * [56] M. Ferronato, N. Castelletto, G. Gambolati, A fully coupled 3-d mixed finite element model of biot consolidation, 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 flow and 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 by employing 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 with evolving 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 Review Applied 5 (4) (2016) 1–30.
  * [64] Y. Zhao, J. Choo, Stabilized material point methods for coupled large deformation and fluid flow in porous 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 hydrologic system: A coupled thermohydrochemical model incorporating buoyant convection, Water Resources Research 49 (9) (2013) 5880–5899.
  * [68] A. Mikelic, M. Wheeler, Convergence of iterative coupling for coupled flow and geomechanics, Computational Geosciences 17 (3) (2013) 455–461.
  * [69] B. Tjaden, S. Cooper, D. Brett, D. Kramer, P. Shearing, On the origin and application of the Bruggeman correlation for analysing transport phenomena in electrochemical systems, Current Opinion in Chemical Engineering 12 (2016) 44–51.
  * [70] D. Mu, Z. Liu, C. Huang, N. Djilali, Determination of the effective diffusion coefficient in porous media 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 for simulation 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 and poromechanics 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 modified salinity 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 multiphase fluid flow, heat transfer, and deformation in fractured porous rock, International Journal of Rock Mechanics and Mining Sciences 39 (4) (2002) 429–442.
  * [76] J. Rutqvist, O. Stephansson, The role of hydromechanical coupling in fractured rock engineering, Hydrogeology Journal 11 (1) (2003) 7–40.
  * [77] K. Min, J. Rutqvist, C. Tsang, L. Jing, Stress-dependent permeability of fractured rock masses: a numerical 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 colloidal particles in natural porous media, Colloids and Surfaces A: Physicochemical and Engineering Aspects 191 (1-2) (2001) 179–188.
  * [79] B. Bijeljic, M. Blunt, Pore-scale modeling of transverse dispersion in porous media, Water Resources Research 43 (12) (2007).
  * [80] Y. Yortsos, D. Salin, On the selection principle for viscous fingering in porous media, Journal of Fluid Mechanics 557 (2006) 225–236.
  * [81] A. Ern, A. Stephansen, P. Zunino, A discontinuous Galerkin method with weighted averages for advection-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 equations approximated 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, IBM journal of Research and Development 11 (2) (1967) 215–234.
  * [84] Z. Ibrahim, K. Othman, M. Suleiman, Implicit r-point block backward differentiation formula for solving 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 stiff ordinary 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 a poroelastic medium, Multiscale Modeling & Simulation 16 (4) (2018) 1542–1580.
  * [87] J. Dormand, P. Prince, A family of embedded Runge-Kutta formulae, Journal of computational and applied mathematics 6 (1) (1980) 19–26.
  * [88] F. Brezzi, M. Fortin, Mixed and hybrid finite element methods, Vol. 15, Springer Science & Business Media, 2012.
  * [89] T. Kadeethum, H. Nick, S. Lee, C. Richardson, S. Salimzadeh, F. Ballarin, A Novel Enriched Galerkin Method for Modelling Coupled Flow and Mechanical Deformation in Heterogeneous Porous Media, in: 53rd US Rock Mechanics/Geomechanics Symposium, American Rock Mechanics Association, New York, 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 production, 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–diffusion equation, 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 transport and 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 element methods and the Galerkin method with bubble functions, Computer Methods in Applied Mechanics and 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 schemes and convex limiting for hyperbolic systems, Computer Methods in Applied Mechanics and Engineering 347 (2019) 143–175.
  * [97] G. Scovazzi, M. Wheeler, A. Mikelic, S. Lee, Analytical and variational numerical methods for unstable miscible 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 with entropy 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. Eijkhout, 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 http://www.mcs.anl.gov/petsc

  * [103] H. Nick, R. Schotting, M. Gutierrez-Neri, K. Johannsen, Modeling transverse dispersion and variable density flow in porous media, Transport in porous media 78 (1) (2009) 11–35.