∎

11institutetext: Saumik Dana 22institutetext: University of Southern
California, Los Angeles, CA 90007
22email<EMAIL_ADDRESS>
33institutetext: Xiaoxi Zhao 44institutetext: University of Southern
California, Los Angeles, CA 90007
44email<EMAIL_ADDRESS>
55institutetext: Birendra Jha 66institutetext: University of Southern
California, Los Angeles, CA 90007
66email<EMAIL_ADDRESS>

# Two-grid method on unstructured tetrahedra: Applying computational geometry
to staggered solution of coupled flow and mechanics problems

Saumik Dana  Xiaoxi Zhao  Birendra Jha

(Received: date / Accepted: date)

###### Abstract

We develop a computational framework that leverages the features of
sophisticated software tools and numerics to tackle some of the pressing
issues in the realm of earth sciences. The algorithms to handle the physics of
multiphase flow, concomitant geomechanics all the way to the earth’s surface
and the complex geometries of field cases with surfaces of discontinuity are
stacked on top of each other in a modular fashion which allows for easy use to
the end user. The current focus of the framework is to provide the user with
tools for assessing seismic risks associated with energy technologies as well
as for use in generating forward simulations in inversion analysis from data
obtained using GPS and InSAR. In this work, we focus on one critical aspect in
the development of the framework: the use of computational geometry in a two-
grid method for unstructured tetrahedral meshes

###### Keywords:

Energy technologies Coupled multiphase flow and geomechanics Computational
geometry Staggered solution algorithm

##  1 Introduction and motivation

Figure 1: Direct imposition of heuristic overburden pressures on the flow
domain completely disregards the mechanical behavior of the surrounding rock
and obviates a study of fault slip away from the reservoir as well as
deformation of earth’s surface Figure 2: The two-grid method enables the study
of induced seismicity in faults away from reservoir, takes into account the
mechanical behavior of surrounding rock and allows for determination of
surface deformation due to subsurface pressure perturbation Figure 3:
Depiction of couple of intersecting flow-geomechanics element pairs. We refer
to the flow element as $E^{f}$ and geomechanics element as $E^{p}$. Depiction
is to reiterate the point that there is no restriction of whether the flow or
geomechanics element needs to be smaller than the other

The development of a computationally inexpensive framework with the capability
to model deformation of the earth’s surface and fault slip due to deep
subsurface pressure perturbations associated with multiphase flow is critical
from two different standpoints

  * •

Assessing the seismic risks associated with carbon sequestration, enhanced
geothermal systems, waste water disposal, enhanced oil recovery thereby
offering guiding protocols for the design of these operations induced1 ;
induced2 ; induced3 ; one ; two ; three ; geo1 ; geo2

  * •

Estimating subsurface properties using inversion analysis of ground
deformation data obtained from Global Positioning System (GPS) and
interferometric synthetic aperture radar (InSAR). This typically requires
multiple forward simulations with the geomechanical domain extending all the
way to the earth’s surface 1 ; 2 ; 3 ; 4 ; 5 ; 6 ; 7 ; 8 ; 9 ;
jha2015reservoir

The major bugbear to developing such a framework is the computationally
intractable size of the geomechanical grid no matter which numerical method is
used to resolve the coupled system of equations. The typical approaches to
avoid such issues are to impose an overburden pressure directly on the
reservoir thus treating it as a coupled problem domain (see Fig. 1) or to
model flow on huge domain with zero permeability cells mimicking the no flow
boundary condition between the flow and non-flow region. The former approach
precludes a study of surface deformation, does not mimic the true effect of
the overburden on the stress sensitive reservoir, and is incapable of
capturing induced seismicity inside faults away from reservoir whereas the
latter approach is computationally intractable for large field scale problems
due to memory requirements. In order to address this, we develop a two-grid
coupled multiphase flow and geomechanics framework which allows for spatial
decoupling of the flow and geomechanics domains with the geomechanics
subproblem being resolved on a separate grid with a larger spatial extent
going all the way to the free surface (see Fig. 2). This computational
framework is built on top of a staggered solution algorithm that solves the
flow and mechanics subproblems sequentially and iteratively.

Typically, in such problems, the geomechanics mesh is expected to be coarser
than the flow mesh everywhere the two meshes exist, but we generalize the
method to the cases where the geomechanics elements can be smaller than the
flow elements in small localized region where capturing the mechanics is more
pertinent. Furthermore, we would like to generalize the notion of the two-grid
from structured hexahedral meshes (as was done in dana-2018 ) to unstructured
tetrahedral meshes. A depiction of the intersection of two tetrahedral
elements in given in Fig. 3. In this document, we focus on demonstrating the
convergence of the two-grid method for unstructured tetrahedral grids using
the classical Mandel’s problem mandel-1953 ; ref27  analytical solution. We
work two cases, one in which we have a finer mesh for flow, and the other in
which we have a finer mesh for mechanics.

This paper is structured as follows: the governing equations are provided in
section 2, the solution strategy is explained in section 3, the salient
features of our software are elucidated in section 4, the numerical
simulations for Mandel’s problem are provided in section 5, and the
conclusions and outlook are provided in section 6.

##  2 Governing equations

We use a classical continuum representation in which the fluids and the solid
skeleton are viewed as overlapping continua BeaJ1972 ; CouO2005 . The
governing equations for coupled flow and geomechanics are obtained from
conservation of mass and balance of linear momentum. We assume that the
deformations are small, that the geomaterial is isotropic, and that the
conditions are isothermal.

Under the quasistatic assumption for earth displacements, the governing
equation for linear momentum balance of the solid/fluid system can be
expressed as

$\nabla\cdot\boldsymbol{\sigma}+\rho_{b}\boldsymbol{g}=\boldsymbol{0},$ (2.1)

where $\boldsymbol{\sigma}$ is the Cauchy total stress tensor,
$\boldsymbol{g}$ is the gravity vector, and
$\rho_{b}=\phi\sum_{\beta}^{n_{\text{phase}}}\rho_{\beta}S_{\beta}+(1-\phi)\rho_{s},$
is the bulk density, $\rho_{\beta}$ and $S_{\beta}$ are the density and
saturation of fluid phase $\beta$, $\rho_{s}$ is the density of the solid
phase, $\phi$ is the true porosity, and $n_{\text{phase}}$ is the number of
fluid phases. The true porosity is defined as the ratio of the pore volume
$V_{p}$ to the bulk volume $V_{b}$ in the current (deformed) configuration.

Assuming that the fluids are immiscible, the mass-conservation equation for
each phase $\alpha$ is

$\frac{dm_{\alpha}}{dt}+\nabla\cdot{\boldsymbol{w}}_{\alpha}=\rho_{\alpha}f_{\alpha},$
(2.2)

where the accumulation term ${dm_{\alpha}}/{dt}$ describes the time variation
of fluid mass relative to the motion of the solid skeleton,
${\boldsymbol{w}}_{\alpha}$ is the mass-flux of fluid phase $\alpha$ relative
to the solid skeleton, and $f_{\alpha}$ is the volumetric source term for
phase $\alpha$. Balance equations (2.1) and (2.2) are coupled by virtue of
poromechanics. On one hand, changes in the pore fluid pressure lead to changes
in effective stress, and induce deformation of the porous material—such as
ground subsidence caused by groundwater withdrawal. On the other hand,
deformation of the porous medium affects fluid mass content and fluid
pressure. The simplest model of this two-way coupling is Biot’s macroscopic
theory of poroelasticity BioM1941 ; GeeJ1957 ; CouO1995 . The specific
equations of single phase poromechanics are provided in Appendix A

##  3 Solution strategy

Figure 4: Two grid staggered solution algorithm

The two-grid approach is built on top of a staggered solution algorithm (Fig.
4) in which the flow and geomechanics subproblems are solved sequentially and
iteratively using a fixed stress split iterative scheme jha2014coupled ; kim-
stability ; ref12 ; mikelic-2014 ; nicola-fixedstress ; white-block ; scalable
; danacg ; danacmame ; dana2020 ; almani ; almanicg ; both2019 ; white2019 ;
borregales2019 ; storvik2019 ; radu-robust ; raduadaptive ; lu_wheeler .

###  3.1 Space and time discretization

We use the finite volume method for the discretization of the flow problem
AziK1979 , and the nodal-based finite element method for the discretization of
the mechanics problem HugT1987 ; ZieO2005 . This space discretization is
locally mass conservative at the element level, and enjoys excellent stability
properties JhaB2007 ; PhiP2007-01 ; PhiP2007-02 ; KimJ2011-spej . The
pressures and saturations degrees of freedom are located at the element
centers of the flow grid, and the displacement vector degrees of freedom are
located at the element nodes of the mechanics grid. In quasi-static
poromechanics, the time derivative appears only in the accumulation term of
the fluid mass balance equation, and we treat it using a fully implicit
Backward Euler time integration scheme.

###  3.2 Two-grid operators

for Loop over all geomechanics elements do

for Loop over all flow elements do

count1 $\leftarrow$ 0

for Loop over four vertices of flow element do

if Vertex inside geomechanics element then

count1 $\leftarrow$ count1 + 1

if count1!=0 then

Flow to geomechanics projection operator $\leftarrow Meas(E^{f})$

Geomechanics to flow projection operator $\leftarrow Meas(E^{p})$

count2 $\leftarrow$ 0

for Loop over four vertices of geomechanics element do

if Vertex inside flow element then

count2 $\leftarrow$ count2 + 1

if count2!=0 then

Flow to geomechanics projection operator $\leftarrow Meas(E^{f})$

Geomechanics to flow projection operator $\leftarrow Meas(E^{p})$

Algorithm 1 Determining the relationship between flow-geomechanics element
pairs

The basic idea of the two-grid approach is to transfer pore pressures from the
flow domain to the geomechanics domain and conversely transfer volumetric
strain from the geomechanical domain to the flow domain. The transfered pore
pressure at every geomechanics element is a volume average of the pore
pressures at the flow elements interacting with the geomechanics element.
Likewise, the transfered volumetric strain at each flow element is a volume
average of the volumetric strains at the geomechanics elements interacting
with the flow element. The global two-grid operators are constructed in the
pre-processing step by applying a local two-grid procedure to each pair of
intersecting finite elements of the two grids as shown in Algorithm 1. The
relationship between the flow-geomechanics element pair is obtained by
determining the barycentric coordinates of the vertices of the flow element
with respect to the geomechanics element and vice versa as explained in
Algorithm 1.

##  4 Software implementation

Figure 5: Workflow

We developed a coupled multiphase flow and geomechanical simulator by
integrating our version of Stanford’s General Purpose Research Simulator
(GPRS) CaoH2002 ; GPRSMan2010  called flowsim which operates as the flow
simulator into an open source framework called PyLith pylithman2011 ; AaaB2013
which operates as the mechanics simulator. The workflow of our software are
elucidated in Fig. 5. Available seismic and geological information are used to
construct a geomechanical model and populate it with poromechanical
properties. Reservoir flow properties (porosity and permeability) and fluid
properties (fluid density, viscosity, and compressibility) are obtained from
an older uncoupled model and transferred to the geomechanical model using the
Matlab Reservoir Simulation Toolbox lie2019 .

###  4.1 Flow Simulator

flowsim is a general purpose, object-oriented, field-scale simulator for
multiphase (oil, water, and gas) and multicomponent (e.g. methane, ethane,
decane, CO2, N2, and water) flow through subsurface. It treats element
connections through a general connection list, which allows for both
structured and unstructured grids. It is capable of handling complex
production and injection scenarios in the field, such as wells perforated at
multiple depths and flowing under variable rate and pressure controls.

###  4.2 Mechanics Simulator

Figure 6: PyLith requires a finite-element mesh, simulation parameters, and
spatial databases (defining the spatial variation of various parameters). It
writes the solution output to either VTK or HDF5/Xdmf files, which can be
visualized with ParaView or Visit. Post-processing is generally done using the
HDF5 files with Python or Matlab scripts aagaard1

PyLith is a open-source finite element code for the simulation of static and
dynamic large-scale deformation problems aagaard1 ; aagaard2 . It uses an
implicit formulation to solve quasi-static problems and an explicit
formulation to solve dynamic rupture problems. The salient features of Pylith
are elucidated in Fig. 6. Some of the advantages of PyLith are: (a) It is
written using C++ and Python languages and is extendable, (b) it is suitable
for parallel computing, (c) it allows localized deformation along discrete
features, such as faults, (d) it is well integrated with meshing codes, such
as LaGriT for tetrahedral meshes lagrit  and CUBIT for both tetrahedral and
hexahedral meshes cubit

###  4.3 Gridding

We employ different grids for flowsim and PyLith to go with our two-grid
method. The finite element grids are generated using CUBIT cubit  mesh
generation software and exported in the Exodus-II format. Since the numerical
method employed for the flow model is finite volume instead of finite element,
we employ some further processing of the flow grid file to render it for the
finite volume calculations. We process the Exodus-II grid file using a MATLAB
script to generate the equivalent finite volume grid in the domain with
element centroid coordinates, element bulk volumes, and face
transmissibilities in the Corner Point Geometry format eclipse

###  4.4 Coupling

We design a C++ class, iflowsim, to allow communication between the flow and
the mechanics simulators. The two-grid operators built on top of iflowsim to
project information between the two grids. PyLith supports distributed memory
parallelization (Message Passing Interface or MPI) whereas flowsim’s
parallelization is based on the shared memory architecture (Multiprocessing or
OpenMP)

##  5 Mandel’s problem

Figure 7: Schematic for Mandel’s problem

(a) Flow mesh, 5706 elements

(b) Pylith mesh, 3697 elements

Figure 8: Finer mesh for flow compared to Pylith

(a) Pylith mesh, 5706 elements

(b) Flow mesh, 3697 elements

Figure 9: Finer mesh for Pylith compared to flow

(a) Flow mesh: 5706 elements and pylith mesh: 3697 elements

(b) Flow mesh: 3697 elements and pylith mesh: 5706 elements

Figure 10: Comparison of non-monotonic pore pressure response

The Mandel’s problem mandel-1953 ; ref27  has been used as a benchmark problem
for testing the validity of numerical codes of coupled poroelasticity. Its
main feature, the Mandel-Cryer effect, is that the pore pressure at the center
of a loaded specimen rises above its initial value because of the two-way
coupling between fluid flow and solid deformation. It involves a long specimen
of rectangular cross section pressed on one side with an impermeable plate
that applies a constant compressive stress and fixed on two sides using
impermeable roller boundaries (Fig. 7). The fourth side of the cross section
is free from normal and shear stresses (traction-free boundary) and is open to
the atmosphere (constant pressure boundary). The porous medium is saturated
with a slightly compressible fluid, water, with initial pressure set at the
reference value. Since the specimen is long, we assume plane strain
conditions, namely, that the displacement and fluid flux vanish in the z
direction (perpendicular to the 2-D domain).

With these boundary conditions, the three-dimensional equations of
poroelasticity reduce to one-dimensional equations for $\sigma_{xx}(y,t)$ and
$p(y,t)$, which can be solved analytically mandel-1953 ; ref27 . At $t=0^{+}$
, a uniform undrained pressure is generated by the Skempton effect, along with
uniform stress $\sigma_{xx}=-\sigma_{0}$. The specimen expands toward the top
boundary due to the Poisson effect. As time progresses, the pressure near the
top boundary decreases because of fluid drainage, which makes the specimen
more compliant there. If the hydraulic diffusivity is small, the effect of
drainage is not observed immediately near the no-flux bottom boundary. This
results into load transfer of compressive total stress toward the bottom
boundary, in response to which the pressure there continues to rise above its
undrained value. At long times, all excess pressure vanishes and a uniform
horizontal stress, $\sigma_{xx}=-\sigma_{0}$, returns. Hence, the pressure
evolution at points away from the drained boundary is nonmonotonic, a
phenomenon not observed in a purely diffusive process such as that modeled by
the Terzaghi theory, where the pressure is uncoupled from the solid
deformation.

In one simulation, we employ a finer mesh for the flow problem as shown in
Fig. 8 whereas in the other simulation, we employ a finer mesh for the
mechanics problem as shown in Fig. 9. Fig. 10 compares the pressure response
at the right bottom corner $(x=50,y=0)$ at different times.

##  6 Conclusions and outlook

We developed a two-grid method which leverages computational geometry on
unstructured tetrahedral meshes that allows us to solve the flow and
geomechanics subproblems on separate grids. This allows use for flexibility in
the choice of the two separate meshes depending on the needs of the specific
problem. There is no restriction on whether the flow element needs to be finer
or coarser than the geomechanics element in any region of the subdomains. This
method also allows us to have separate domains for flow and geomechanics,
which is critical for large scale field problems typically encountered in
energy technologies. We shall now employ our computational framework to solve
field scale problems with faults in future articles.

###### Acknowledgements.

The first author would like to thank Sreekanth Arikatla (Staff R&D Engineer,
Kitware Inc.) for discussions on computational geometry for unstructured
tetrahedra

## Conflict of interest

The authors declare that they have no conflict of interest.

##  Appendix A Single-phase poromechanics

For isothermal single-phase flow of a slightly compressible fluid in a
poroelastic medium with no stress dependence of permeability, the single-phase
fluid mass conservation equation reduces to

$\frac{dm}{dt}+\nabla\cdot\boldsymbol{w}=\rho_{f}f,$ (A.1)

where $m$ is the fluid mass content (fluid mass per unit bulk volume of porous
medium), $\rho_{f}$ is the fluid density,
$\boldsymbol{w}=\rho_{f}\boldsymbol{v}$ is the fluid mass flux (fluid mass
flow rate per unit area and time), and $\boldsymbol{v}$ is the seepage
velocity relative to the deforming skeleton, given by Darcy’s law:

$\boldsymbol{v}=-\frac{\boldsymbol{k}}{\mu}\left(\nabla
p-\rho_{f}\boldsymbol{g}\right),$ (A.2)

where $\boldsymbol{k}$ is the intrinsic permeability tensor, $\mu$ is the
fluid dynamic viscosity and $p$ is the pore-fluid pressure BeaJ1972 . It is
useful to define the fluid content variation $\zeta$,

$\zeta:=\frac{\delta m}{\rho_{f,0}},$ (A.3)

where $\delta m=m-m_{0}$ is the increment in fluid mass content with respect
to the initial reference state, and $\rho_{f,0}$ is the reference fluid
density. The self-consistent theory of poroelastic behavior proposed by
BioM1941  links the changes in total stress and fluid pressure with changes in
strain and fluid content. Following CouO1995 , the poroelasticity equations
can be written in incremental form as

$\begin{split}\delta\boldsymbol{\sigma}&=\boldsymbol{C}_{dr}:\boldsymbol{\varepsilon}-b\delta
p\boldsymbol{1}\\\ \zeta&=b\varepsilon_{v}+\frac{1}{M}\delta p,\end{split}$
(A.4)

where $\boldsymbol{C}_{dr}$ is the rank-4 drained elasticity tensor,
$\boldsymbol{1}$ is the rank-2 identity tensor, $\boldsymbol{\varepsilon}$ is
the linearized strain tensor, defined as the symmetric gradient of the
displacement vector $\boldsymbol{u}$,

$\boldsymbol{\varepsilon}:=\frac{1}{2}\left(\nabla\boldsymbol{u}+\nabla^{T}\boldsymbol{u}\right),$
(A.5)

and $\varepsilon_{v}=\text{tr}(\boldsymbol{\varepsilon})$ is the volumetric
strain. Note that we use the convention that tensile stress is positive. It is
useful to express the strain tensor as the sum of its volumetric and
deviatoric components:

$\boldsymbol{\varepsilon}=\frac{1}{3}\varepsilon_{v}\boldsymbol{1}+\boldsymbol{e},$
(A.6)

from which it follows that the volumetric stress
$\sigma_{v}=\text{tr}(\boldsymbol{\sigma})/3$ satisfies:

$\delta\sigma_{v}=K_{dr}\varepsilon_{v}-b\delta p.$ (A.7)

Equation (A.4) implies that the effective stress in single-phase
poroelasticity, responsible for skeleton deformation, is defined in
incremental form as

$\delta\boldsymbol{\sigma}^{\prime}:=\delta\boldsymbol{\sigma}+b\delta
p\boldsymbol{1}.$ (A.8)

Biot’s theory of poroelasticity has two coupling coefficients: the Biot
modulus $M$ and the Biot coefficient $b$. They are related to rock and fluid
properties as CouO1995

$\frac{1}{M}=\phi_{0}c_{f}+\frac{b-\phi_{0}}{K_{s}},\quad
b=1-\frac{K_{dr}}{K_{s}},$ (A.9)

where $c_{f}=1/K_{f}$ is the fluid compressibility, $K_{f}$ is the bulk
modulus of the fluid, $K_{s}$ is the bulk modulus of the solid grain, and
$K_{dr}$ is the drained bulk modulus of the porous medium. To set the stage
for the numerical solution strategy of the coupled problem, it is useful to
write the fluid mass balance equation (A.1) (the pressure equation) in a way
that explicitly recognizes the coupling with mechanical deformation. Equations
(A.4) state that the increment in fluid mass content has two components:
increment due to expansion of the pore space and increment due to increase in
the fluid pressure. Assuming small elastic deformations and applying
linearization from the reference state to the current state, we can write Eqs.
(A.4) as

$\displaystyle\boldsymbol{\sigma}-\boldsymbol{\sigma}_{0}$
$\displaystyle=\boldsymbol{C}_{dr}:\boldsymbol{\varepsilon}-b\left(p-p_{0}\right)\boldsymbol{1},$
(A.10) $\displaystyle\frac{1}{\rho_{f,0}}(m-m_{0})$
$\displaystyle=b\varepsilon_{v}+\frac{1}{M}(p-p_{0}).$ (A.11)

Substituting Eq. (A.11) into Eq. (A.1), we obtain the fluid mass balance
equation in terms of the pressure and the volumetric strain:

$\frac{1}{M}\frac{\partial p}{\partial
t}+b\frac{\partial\varepsilon_{v}}{\partial t}+\nabla\cdot\boldsymbol{v}=f.$
(A.12)

Linearizing the relation between volumetric total stress and volumetric strain
with respect to the reference state,

$\sigma_{v}-\sigma_{v,0}=K_{dr}\varepsilon_{v}-b\left(p-p_{0}\right),$ (A.13)

allows us to express the change in porosity as the sum of a volumetric stress
component and a fluid pressure component. From $m=\rho_{f}\phi$ and Eq.
(A.11),

$\frac{\rho_{f}}{\rho_{f,0}}\phi-\phi_{0}=\frac{b}{K_{dr}}\left(\sigma_{v}-\sigma_{v,0}\right)+\left(\frac{b^{2}}{K_{dr}}+\frac{1}{M}\right)\left(p-p_{0}\right).$
(A.14)

Using the effective stress equation, Eq. (A.13), we can rewrite Eq. (A.12) in
terms of pressure and volumetric total stress:

$\left(\frac{b^{2}}{K_{dr}}+\frac{1}{M}\right)\frac{\partial p}{\partial
t}+\frac{b}{K_{dr}}\frac{\partial\sigma_{v}}{\partial
t}+\nabla\cdot\boldsymbol{v}=f.$ (A.15)

## References

  * (1) Aagaard, B., Kientz, S., Knepley, M., Strand, L., Williams, C.: PyLith User Manual, Version 1.8.0.  Computational Infrastructure for Geodynamics, University of California, Davis (2012)
  * (2) Aagaard, B., Kientz, S., Knepley, M., Strand, L., Williams, C.: Pylith user manual: version 2.1. 0.  Davis, CA: Computational Infrastructure of Geodynamics (2013)
  * (3) Aagaard, B.T., Knepley, M.G., Williams, C.A.: A domain decomposition approach to implementing fault slip in finite-element models of quasi-static and dynamic crustal deformation.  J. Geophys. Res. 118, 3059–3079 (2013)
  * (4) Aagaard, B.T., Knepley, M.G., Williams, C.A.: A domain decomposition approach to implementing fault slip in finite-element models of quasi-static and dynamic crustal deformation.  Journal of Geophysical Research: Solid Earth 118(6), 3059–3079 (2013)
  * (5) Abousleiman, Y., Cheng, A.H.D., Cui, L., Detournay, E., Roegiers, J.C.: Mandel’s problem revisited.  Géotechnique 46(2), 187–195 (1996)
  * (6) Ahmed, E., Nordbotten, J.M., Radu, F.A.: Adaptive asynchronous time-stepping, stopping criteria, and a posteriori error estimates for fixed-stress iterative schemes for coupled poromechanics problems.  Journal of Computational and Applied Mathematics 364, 112312 (2020)
  * (7) Almani, T., Kumar, K., Dogru, A., Singh, G., Wheeler, M.F.: Convergence analysis of multirate fixed-stress split iterative schemes for coupling flow with geomechanics.  Computer Methods in Applied Mechanics and Engineering 311, 180–207 (2016)
  * (8) Almani, T., Kumar, K., Wheeler, M.F.: Convergence and error analysis of fully discrete iterative coupling schemes for coupling flow with geomechanics.  Computational Geosciences (2017)
  * (9) Aziz, K., Settari, A.: Petroleum Reservoir Simulation.  Elsevier, London (1979)
  * (10) Bear, J.: Dynamics of Fluids in Porous Media.  Wiley, New York (1972)
  * (11) Biot, M.A.: General theory of three-dimensional consolidation.  J. Appl. Phys. 12, 155–164 (1941)
  * (12) Blacker, T.D., Bohnhoff, W.J., Edwards, T.L.: Cubit mesh generation environment. volume 1: Users manual.  Tech. rep., Sandia National Labs., Albuquerque, NM (United States) (1994)
  * (13) Borregales, M., Kumar, K., Radu, F.A., Rodrigo, C., Gaspar, F.J.: A partially parallel-in-time fixed-stress splitting method for biot’s consolidation model.  Computers and Mathematics with Applications 77(6), 1466 – 1478 (2019).  7th International Conference on Advanced Computational Methods in Engineering (ACOMEN 2017)
  * (14) Both, J.W., Kumar, K., Nordbotten, J.M., Radu, F.A.: Anderson accelerated fixed-stress splitting schemes for consolidation of unsaturated porous media.  Computers and Mathematics with Applications 77(6), 1479 – 1502 (2019).  7th International Conference on Advanced Computational Methods in Engineering (ACOMEN 2017)
  * (15) Cao, H.: Development of techniques for general purpose simulators.  Ph.D. thesis, Stanford University (2002)
  * (16) Castelletto, N., White, J.A., Ferronato, M.: Scalable algorithms for three-field mixed finite element coupled poromechanics.  Journal of Computational Physics 327, 894–918 (2016)
  * (17) Castelletto, N., White, J.A., Tchelepi, H.A.: Accuracy and convergence properties of the fixed-stress iterative solution of two-way coupled poromechanics.  International Journal for Numerical and Analytical Methods in Geomechanics 39(14), 1593–1618 (2015)
  * (18) Chang, H., Chen, Y., Zhang, D., et al.: Data assimilation of coupled fluid flow and geomechanics using the ensemble kalman filter.  SPE Journal 15(02), 382–394 (2010)
  * (19) Council, N.R., et al.: Induced seismicity potential in energy technologies.  National Academies Press (2013)
  * (20) Coussy, O.: Mechanics of Porous Continua.  John Wiley and Sons, Chichester, England (1995)
  * (21) Coussy, O.: Poromechanics of freezing materials.  J. Mech. Phys. Solids 53, 1689–1718 (2005)
  * (22) Dana, S., Ganis, B., Wheeler, M.F.: A multiscale fixed stress split iterative scheme for coupled flow and poromechanics in deep subsurface reservoirs.  Journal of Computational Physics 352, 1–22 (2018)
  * (23) Dana, S., Ita, J., Wheeler, M.F.: The correspondence between voigt and reuss bounds and the decoupling constraint in a two-grid staggered algorithm for consolidation in heterogeneous porous media.  Multiscale Modeling & Simulation 18(1), 221–239 (2020)
  * (24) Dana, S., Wheeler, M.F.: Convergence analysis of fixed stress split iterative scheme for anisotropic poroelasticity with tensor biot parameter.  Computational Geosciences 22(5), 1219–1230 (2018)
  * (25) Dana, S., Wheeler, M.F.: 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, 788–806 (2018)
  * (26) Galloway, D.L., Hoffmann, J.: The application of satellite differential sar interferometry-derived ground displacements in hydrogeology.  Hydrogeology Journal 15(1), 133–154 (2007)
  * (27) Galloway, D.L., Hudnut, K.W., Ingebritsen, S., Phillips, S.P., Peltzer, G., Rogez, F., Rosen, P.: Detection of aquifer system compaction and land subsidence using interferometric synthetic aperture radar, antelope valley, mojave desert, california.  Water Resources Research 34(10), 2573–2585 (1998)
  * (28) Gaucher, E., Schoenball, M., Heidbach, O., Zang, A., Fokker, P.A., van Wees, J.D., Kohl, T.: Induced seismicity in geothermal reservoirs: A review of forecasting approaches.  Renewable and Sustainable Energy Reviews 52, 1473 – 1490 (2015)
  * (29) Geertsma, J.: The effect of fluid pressure decline on volumetric change of porous rocks.  Trans. AIME 210, 331–340 (1957)
  * (30) GeoQuest, S.: Eclipse reference manual.  Schlumberger, Houston, Texas (2014)
  * (31) Hesse, M.A., Stadler, G.: Joint inversion in coupled quasi-static poroelasticity.  Journal of Geophysical Research: Solid Earth 119(2), 1425–1445 (2014)
  * (32) Hughes, T.J.R.: The Finite Element Method: Linear Static and Dynamic Finite Element Analysis.  Prentice-Hall, Englewood Cliffs, NJ (1987)
  * (33) Iglesias, M.A., McLaughlin, D.: Data inversion in coupled subsurface flow and geomechanics models.  Inverse Problems 28(11), 115009 (2012)
  * (34) Jha, B., Bottazzi, F., Wojcik, R., Coccia, M., Bechor, N., McLaughlin, D., Herring, T., Hager, B.H., Mantica, S., Juanes, R.: Reservoir characterization in an underground gas storage field using joint inversion of flow and geodetic data.  International Journal for Numerical and Analytical Methods in Geomechanics 39(14), 1619–1638 (2015)
  * (35) Jha, B., Juanes, R.: A locally conservative finite element framework for the simulation of coupled flow and reservoir geomechanics.  Acta Geotechnica 2, 139–153 (2007)
  * (36) Jha, B., Juanes, R.: Coupled multiphase flow and poromechanics: A computational model of pore pressure effects on fault slip and earthquake triggering.  Water Resources Research 50(5), 3776–3808 (2014)
  * (37) Kim, J., Tchelepi, H., Juanes, R.: Stability, accuracy and efficiency of sequential methods for coupled flow and geomechanics.  SPE Journal 16(2), 249–262 (2011)
  * (38) Kim, J., Tchelepi, H.A., Juanes, R.: Stability, accuracy and efficiency of sequential methods for coupled flow and geomechanics.  Soc. Pet. Eng. J. 16(2), 249–262 (2011)
  * (39) Lackner, K.S.: A guide to co2 sequestration.  Science 300(5626), 1677–1678 (2003)
  * (40) LaGriT: Los alamos grid toolbox,(lagrit) los alamos national laboratory (2013)
  * (41) Lie, K.A.: An introduction to reservoir simulation using MATLAB/GNU Octave: User guide for the MATLAB Reservoir Simulation Toolbox (MRST).  Cambridge University Press (2019)
  * (42) Lu, X., Wheeler, M.F.: Three-way coupling of multiphase flow and poromechanics in porous media.  Journal of Computational Physics 401, 109053 (2020)
  * (43) Majer, E., Nelson, J., Robertson-Tait, A., Savy, J., Wong, I.: Protocol for addressing induced seismicity associated with enhanced geothermal systems.  US Department of Energy p. 52 (2012)
  * (44) Mandel, J.: Consolidation des sols (étude mathématique)*.  Géotechnique 3(7), 287–299 (1953)
  * (45) Mikelić, A., Wang, B., Wheeler, M.F.: Numerical convergence study of iterative coupling for coupled flow and geomechanics.  Computational Geosciences 18(3), 325–341 (2014)
  * (46) Mikelić, A., Wheeler, M.: Convergence of iterative coupling for coupled flow and geomechanics.  Computational Geosciences 17(3), 455–461 (2013)
  * (47) Orr, F.M.: Onshore geologic storage of co2.  Science 325(5948), 1656–1658 (2009)
  * (48) Pan, H., Cao, H.: User Manual for General Purpose Research Simulator.  Stanford University Petroleum Engineering Institute, Stanford, CA (2010)
  * (49) Phillips, P.J., Wheeler, M.F.: A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I: the continuous in time case.  Comput. Geosci. 11, 131–144 (2007)
  * (50) Phillips, P.J., Wheeler, M.F.: A coupling of mixed and continuous Galerkin finite element methods for poroelasticity II: the discrete-in-time case.  Comput. Geosci. 11, 145–158 (2007)
  * (51) Radu, F.A., Nordbotten, J.M., Pop, I.S., Kumar, K.: A robust linearization scheme for finite volume based discretizations for simulation of two-phase flow in porous media.  Journal of Computational and Applied Mathematics 289, 134 – 141 (2015)
  * (52) Rutqvist, J., Rinaldi, A.P., Cappa, F., Moridis, G.J.: Modeling of fault reactivation and induced seismicity during hydraulic fracturing of shale-gas reservoirs.  Journal of Petroleum Science and Engineering 107, 31–44 (2013)
  * (53) Storvik, E., Both, J.W., Kumar, K., Nordbotten, J.M., Radu, F.A.: On the optimization of the fixed‐stress splitting for biot’s equations.  International Journal for Numerical Methods in Engineering p. nme.6130 (2019)
  * (54) Szulczewski, M.L., MacMinn, C.W., Herzog, H.J., Juanes, R.: Lifetime of carbon capture and storage as a climate-change mitigation technology.  Proceedings of the National Academy of Sciences 109(14), 5185–5189 (2012)
  * (55) Vasco, D.: Estimation of flow properties using surface deformation and head data: A trajectory-based approach.  Water Resources Research 40(10) (2004)
  * (56) Vasco, D., Karasaki, K., Kishida, K.: A coupled inversion of pressure and surface displacement.  Water resources research 37(12), 3071–3089 (2001)
  * (57) Vasco, D., Rucci, A., Ferretti, A., Novali, F., Bissell, R., Ringrose, P., Mathieson, A., Wright, I.: Satellite-based measurements of surface deformation reveal fluid flow associated with the geological storage of carbon dioxide.  Geophysical Research Letters 37(3) (2010)
  * (58) White, J.A., Castelletto, N., Klevtsov, S., Bui, Q.M., Osei-Kuffuor, D., Tchelepi, H.A.: A two-stage preconditioner for multiphase poromechanics in reservoir simulation.  Computer Methods in Applied Mechanics and Engineering 357, 112575 (2019)
  * (59) White, J.A., Castelletto, N., Tchelepi, H.A.: Block-partitioned solvers for coupled poromechanics: A unified framework.  Computer Methods in Applied Mechanics and Engineering 303, 55–74 (2016)
  * (60) Wilschut, F., Peters, E., Visser, K., Fokker, P.A., van Hooff, P., et al.: Joint history matching of well data and surface subsidence observations using the ensemble kalman filter: a field study.  In: SPE Reservoir Simulation Symposium. Society of Petroleum Engineers (2011)
  * (61) Zang, A., Oye, V., Jousset, P., Deichmann, N., Gritto, R., McGarr, A., Majer, E., Bruhn, D.: Analysis of induced seismicity in geothermal reservoirs – an overview.  Geothermics 52, 6 – 21 (2014)
  * (62) Zienkiewicz, O.C., Taylor, R.L.: The Finite Element Method for Solid and Structural Mechanics.  Elsevier, Burlington, MA (2005)