# Unexpected fault activation in underground gas storage. Part I: Mathematical
model and mechanisms

Andrea Franceschini<EMAIL_ADDRESS>Claudia Zoccarato  Selena
Baldan  Matteo Frigo  Massimiliano Ferronato  Carlo Janna  Giovanni Isotton
Pietro Teatini

###### Abstract

Underground gas storage (UGS) is a worldwide well-established technology that
is becoming even more important to cope with seasonal peaks of gas consumption
due to the growing uncertainties of the energy market. Safety issues
concerning the reactivation of pre-existing faults might arise if the target
reservoir is located in a faulted basin, where human activities can trigger
(micro-)seismicity events. From a mechanical viewpoint, a fault is activated
when the shear stress exceeds the limiting frictional value. In the
Netherlands, it has been observed that this occurrence can develop somehow
“unexpectedly” after the primary production (PP), i.e., during cushion gas
injection (CGI) and UGS cycles, when the stress regime should be in the
unloading/reloading path and the fault state far from failure. In order to
understand the physical mechanisms responsible for such occurrences and build
reliable simulation tools for predictive purposes, a 3D mathematical model
coupling frictional contact mechanics in faulted porous rocks with fluid flow
is developed, implemented and tested. In particular, the mechanisms and the
critical factors responsible for the fault reactivation during the various UGS
stages are investigated in the real-world setting of the Rotliegend formation
in the Netherlands. The effect of the storage of different fluids for various
purposes, such as the long-term sequestration of CO2, the regular injection
and extraction cycles of CH4, and the highly irregular cycles of H2, is
investigated with respect to fault activation risk. The final aim of this two-
part work is to define a safe operational bandwidth for the pore pressure
range for UGS activities in the faulted reservoirs of the Rotliegend
formation.

Part I of this work concerns the development of the mathematical and numerical
model of frictional contact mechanics and flow in faulted porous rocks. A
mixed discretization of the governing PDEs under frictional contact
constraints along the faults is used, where displacement and pressure in the
porous medium, and traction on the fault surfaces are the main variables. A
slip-weakening constitutive law governing the fault macroscopic behavior is
also presented. The model is tested in the setting of an ideal reservoir
located in the Rotliegend formation. The analyses point out how fault
reactivation during PP can lead to a stress redistribution, giving rise to a
new (deformed) equilibrium configuration. When the fault is reloaded in the
opposite direction during the CGI and/or UGS stages, further activation events
can occur even if the stress range does not exceed either the undisturbed
initial value or the maximum strength ever experienced by the formation.

###### keywords:

Frictional contact , Mixed discretization , Underground gas storage , Slip-
weakening law , Fault reactivation

††journal: J. Comput. Phys.

[UNIPD]organization=Department of Civil, Environmental and Architectural
Engineering, University of Padova, city=Padova, postalcode=35131,
country=Italy [M3E]organization=M3E S.r.l., city=Padova, postalcode=35121,
country=Italy

##  1 Introduction

Seismicity associated to fluid withdrawal from and injection into deep
reservoirs is a geomechanical hazard that is receiving a growing attention in
the scientific literature [1, 2, 3]. Fault reactivation, both aseismic and
seismic, is caused by the change of the natural stress regime on the
discontinuity surface due to the pore pressure $p$ changes in the reservoir
where mining activities are operated. More specifically, the onset and amount
of slip, and the size of the reactivated fault zone depend on how the stress
changes caused by the human operations at depth can interfere with the natural
stress regime [4, 5, 6].

The current state-of-the-art research on this topic focuses on the two main
processes: i) seismicity induced by production of (conventional) hydrocarbon
reservoirs, where pore pressure depletion $\Delta p$ and differential
reservoir compaction are the main factors yielding fault reactivation [4, 7]
(Fig. 1a); ii) fluid injection at depth (CO2 sequestration, production from
unconventional reservoirs, enhanced geothermal systems) where, independently
of the possible thermal processes, the increase of the fluid pressure
(largely) above the natural undisturbed value $p_{i}$ within the faulted zone
crossing or bounding the targeted formation drives the reactivation of rock
discontinuities [8, 9, 10, 11] (Fig. 1b).

Figure 1: Sketches of two (a and b) “expected” induced seismicity scenarios
and one (c) “unexpected”. a) Primary production with large pressure drop, b)
fluid injection (CO2 sequestration, waste water disposal, fracking) with
significant pressure increase), c) pressure in the range already experienced
(UGS with $p<p_{i}$).

Over the last decade, induced seismicity has been observed in some parts of
the world also in reservoirs used for underground gas storage (UGS). Somehow
unexpectedly, fault reactivation occurred not only during primary production
(PP) or gas storage at pressure larger than $p_{i}$ [12, 13, 14], i.e., at a
stress regime that had never been experienced before by the reservoir and the
nearby faults, but also during cushion gas injection (GCI) or producing and
storing phases with a pore pressure smaller than $p_{i}$ and larger than
$p_{\mathrm{min}}$, i.e., the minimum pressure experienced by the field
usually at the end of primary production before its conversion to UGS (Fig.
1c) [15, 16, 17, 18].

The present work is aimed to shed light on these “unexpected” events. Because
of the current CH4 importance for energy production purposes and the
international turbulence on this market, the interest in developing UGS
projects is increasing worldwide. Multiple elements presently characterize
UGS: seasonal and short-term balancing, strategic reserves in case of
interruption of deliveries, optimisation of gas production and gas system
distribution, overcoming of local restrictions of gas grids [19, 20]. More
recently, UGS has also been investigated as a possible method to store green
energy in terms of compressed air and H2 [21, 22, 23]. Sources of green
energy, such as wind, waves, and sun, are characterized by a natural high-
frequency fluctuation (from hours to day/night and to weeks). Excess
electricity can be used to synthesize hydrogen or to compress air, store the
gas in deep aquifers or depleted reservoirs, and use it at a later stage as
fuel to generate electricity. The same technology can be also applied for
long-term geological sequestration activities, such as CO2 capture and
sequestration to reduce carbon dioxide emissions in the atmosphere [24]. In
this case, the targeted reservoir does not undergo a cyclic loading/unloading
strength, but the pressure can increase up to a steady state value usually
smaller or equal to the initial value $p_{i}$. Analysis of the social and
environmental hazards and risks associated with subsurface gas storage is a
recurrent issue whenever a new UGS site is planned. Many different aspects are
involved, such as formation integrity, health and safety as related to public
perception, economic risk, and environmental impact. Among the latter, the
geomechanical effects induced by seasonal gas injection and withdrawal, such
as movements of the land surface, may play an important role [25].

UGS has been rarely associated with induced seismicity. According to data
provided by the European Commission [26], recent works [27], and the HiQuake
database [28], only a few sites have reported human-induced earthquakes out of
160 UGS facilities in Europe and more than 380 in USA [2]. Three cases of
these, i.e., Bergermeer, Norg, and Grijpskerk fields, are located in The
Netherlands. These reservoirs are located in the Carboniferous-Rotliegend
formation, northern Europe, which is one of the most intensively explored
petroleum systems in the world [29], and where a relatively high number of
induced seismicity events has been recorded over the last decades [30, 31].
Several studies have addressed the topic of fault reactivation in Rotliegend
reservoirs, the most famous of which concerns the Groningen field [32]. Most
studies focus on a specific reservoir in the Netherlands and northern Germany,
or, more generally, try to investigate the relationship between the typical
geological features of these reservoirs, their usual production life, and the
possible induced seismicity [7, 33, 34, 35, 6]. The recent literature,
however, is mainly concerned with primary production only and does not
investigate the reasons why fault reactivation can occur during UGS phases.
Moreover, a very simplified geological structure is assumed in such analyses,
with a single fault in a two dimensional (2D) vertical plane, and most likely
this can only partially capture the complex response expected from many
intersecting faults in a fully three dimensional (3D) environment [34]. Only a
few relatively old publications addressed the topic in UGS reservoirs [36,
37]. Nagelhout and Roest [36] developed a 2D geomechanical model by means of
the FLAC simulator for a typical faulted vertical section and concluded that
while the gas field is depleted, fault slip occurs due to compaction of the
reservoir and due to the upward movement of strata underlying the reservoir.
Negligible amounts of additional slip are induced when the reservoir is
subjected to alternating injection/extraction periods. Orlic et al. [37]
simulated the geomechanical behavior of a specific UGS reservoir using the
finite-element package DIANA. Their results highlighted that the critically
stressed section of the central fault affected by the fault slipped … during
gas production. Additional fault slip could be expected during the subsequent
phase of cushion gas injection … During annual cycles of gas injection and
production, the central fault is not critically stressed anymore.

The aim of this work is multifold: i) to develop a robust computational
framework allowing for the simulation of the inception of fault activation in
3D real-world geological settings; ii) to improve the understanding of the
physical mechanisms underlying induced seismicity during UGS activities, with
specific reference to the typical configurations of Dutch UGS reservoirs; iii)
to investigate the factors that can increase the chance of fault reactivation
during UGS activities, identifying the settings, conditions and material
properties that could most likely cause “unexpected” fault reactivation in the
reservoirs located in the Rotliegend formations; iv) to define a set of
practical guidelines allowing for a safe operational bandwidth in such UGS
fields, in consideration also of the different potential storage activities
(CH4, H2, CO2). A few preliminary outcomes were already reported in [38] and
[39]. In order to accomplish such a complex multi-disciplinary task, the
overall work is subdivided into two parts. The present paper (Part I) is
mainly concerned with objectives i) and ii), focusing on the mathematical and
computational aspects of the modeling approach, on its application in a
representative 3D test case of the problem of interest, and on the mechanisms
that can cause “unexpected” fault activation during UGS activities. The
application of the model developed herein to the specific real-world cases of
the Rotliegend formation, with a detailed sensitivity analysis for the
different storage activities and the definition of preliminary guidelines
(aforementioned objectives iii) and iv)), is the target of Part II [40].

The paper is organized as follows. The mathematical model of frictional
contact mechanics and flow in a 3D visco-elasto-plastic porous medium, built
on top of the works [41, 42, 43], is introduced along with its numerical
discretization and solution algorithms. Faults are explicitly simulated within
the porous rock as inner contact boundaries, whose activation is
macroscopically governed by Coulomb’s criterion. Pressure change within the
faults, variation of Coulomb’s parameters due to slip-weakening, and the
rheology of the caprock are properly accounted for. The model is applied to a
synthetic reservoir and fault system that realistically represents the main
geological features of the Rotliegend reservoirs. Two scenarios are simulated
to deepen the understanding of the geomechanical behavior of a faulted UGS
system. Computational results are presented and the mechanisms responsible for
fault reactivation during UGS phases are pointed out. A few conclusive remarks
close the presentation.

##  2 Mathematical and numerical model

In this section, we discuss the development of the mathematical and numerical
model used to investigate the fault activation in the context of UGS
reservoirs. The aim is at solving the frictional contact mechanics problem for
a faulted porous medium, where the constraints are imposed in an exact way by
Lagrange multipliers. The friction behavior of the fracture is governed by
Coulomb’s criterion, with a slip-weakening constitutive law. The variational
formulation, its numerical discretization and the possible related instability
phenomena are discussed. The pore pressure, both in the continuous matrix and
inside the fracture network, is computed by a flow simulator with a one-way
coupled approach [44, 45], which turns out to be fully warranted at the space
and time scale of interest. We use the quasi-static assumption, i.e., no
acceleration contribution is accounted for, under the hypothesis of likely
negligible inertia of the system when small (e.g., centimetric) slip and small
areal extent characterize the fault reactivation [46].

###  2.1 Strong formulation for the contact problem

A fault can be modeled at the macroscale as a lower dimensional internal
boundary $\Gamma_{f}$ embedded in a 3D domain $\Omega\subset\mathbb{R}^{3}$.
The fracture is represented as a pair of surfaces in contact, conventionally
denoted as top and bottom and represented by $\Gamma_{f}^{+}$ and
$\Gamma_{f}^{-}$, respectively. On such surfaces, normal and frictional
contact conditions have to be enforced, like the impenetrability of solid
bodies and the fulfillment of a friction criterion. In this work, we use
Coulomb’s frictional criterion to provide the limiting modulus for the shear
component of traction on the fault surface. To complete the problem setting
definition, we introduce the external domain boundary
$\Gamma\equiv\partial\Omega$, with its outer unit normal vector
$\boldsymbol{n}$, while
$\boldsymbol{n}_{f}=\boldsymbol{n}_{f}^{-}=-\boldsymbol{n}_{f}^{+}$ denotes
the normal direction to the fracture surface $\Gamma_{f}$. Fig. 2 shows a
sketch of the domain $\Omega$, the fault $\Gamma_{f}$ and the related
quantities. Any vector field can be decomposed along the normal and tangential
direction to the fracture, i.e.,
$\boldsymbol{v}=v_{N}\boldsymbol{n}_{f}+\boldsymbol{v}_{T}$, with
$v_{N}=\boldsymbol{n}_{f}^{T}\boldsymbol{v}$ and
$\boldsymbol{v}_{T}=\boldsymbol{v}-v_{N}\boldsymbol{n}_{f}=(\boldsymbol{1}-\boldsymbol{n}_{f}\otimes\boldsymbol{n}_{f})\,\boldsymbol{v}$,
where the subscripts $N$ and $T$ are used to denote the normal and tangential
components, respectively, and $\boldsymbol{1}$ is the identity tensor of order
2.

Figure 2: Sketch of the 3D domain $\Omega$ with its boundary, outer normal and
inner fracture $\Gamma_{f}$ (left), made of the top and bottom contact
surfaces and the normal direction $\boldsymbol{n}_{f}$ (right).

Assuming quasi-static conditions and infinitesimal strains, the strong form of
the linear momentum balance at every instant $t$ in the time interval
$[0,t_{\max}]$ can be stated as follows [47, 48, 49]: find the displacement
vector
${\boldsymbol{u}:\overline{\Omega}\times[0,t_{\max}]\rightarrow\mathbb{R}^{3}}$
such that:

$\displaystyle\nabla\cdot\hat{\boldsymbol{\sigma}}(\boldsymbol{u})+\boldsymbol{b}$
$\displaystyle=0$ $\displaystyle\text{in }\Omega\times[0,t_{\max}],$ (1a)
$\displaystyle\boldsymbol{u}$ $\displaystyle=\bar{\boldsymbol{u}}$
$\displaystyle\text{on }\Gamma_{u}\times[0,t_{\max}],$ (1b)
$\displaystyle\hat{\boldsymbol{\sigma}}(\boldsymbol{u})\cdot\boldsymbol{n}$
$\displaystyle=\bar{\boldsymbol{t}}$ $\displaystyle\text{on
}\Gamma_{\sigma}\times[0,t_{\max}],$ (1c)

where $\hat{\boldsymbol{\sigma}}$ is the total stress tensor, $\boldsymbol{b}$
collects the external body loads and $\Gamma_{u}\cup\Gamma_{\sigma}=\Gamma$,
$\Gamma_{u}\cap\Gamma_{\sigma}=\varnothing$, are the portion of the boundary
where Dirichlet and Neumann conditions are imposed, respectively. On the
fracture $\Gamma_{f}$, normal and friction compatibility conditions need to be
enforced [48, 49]. The normal contact conditions on the fracture read:

$\displaystyle t_{N}=\boldsymbol{t}\cdot\boldsymbol{n}_{f}$ $\displaystyle\leq
0$ $\displaystyle\text{only compressive traction is allowed},$ (2a)
$\displaystyle
g_{N}=\llbracket\boldsymbol{u}\rrbracket\cdot\boldsymbol{n}_{f}$
$\displaystyle\geq 0$ $\displaystyle\text{impenetrability condition},$ (2b)
$\displaystyle t_{N}g_{N}$ $\displaystyle=0$ $\displaystyle\text{either the
fracture is compressed or it is open}.$ (2c)

The conditions for the frictional component are:

$\displaystyle\left\|\boldsymbol{t}_{T}\right\|_{2}\leq\tau_{\max}\left(t_{N},\|\boldsymbol{g}_{T}\|_{2}\right)$
$\displaystyle\text{Coulomb's criterion},$ (3a)
$\displaystyle\dot{\boldsymbol{g}}_{T}\cdot\boldsymbol{t}_{T}-\tau_{\max}\left(t_{N},\|\boldsymbol{g}_{T}\|_{2}\right)\|\dot{\boldsymbol{g}}_{T}\|_{2}=0$
$\displaystyle\text{frictional traction is aligned with sliding rate}.$ (3b)

In Eqs. (2)-(3), we split the traction $\boldsymbol{t}$ on the fracture and
the displacement jump across it into normal and tangential components, i.e.,
$\boldsymbol{t}=t_{N}\boldsymbol{n}_{f}+\boldsymbol{t}_{T}$ and
$\llbracket\boldsymbol{u}\rrbracket=g_{N}\boldsymbol{n}_{f}+\boldsymbol{g}_{T}$,
respectively. The jump is defined as
$\llbracket\boldsymbol{u}\rrbracket=\boldsymbol{u}|_{\text{top}}-\boldsymbol{u}|_{\text{bottom}}$.
To characterize the standard Coulomb frictional criterion, $c$ and $\varphi$
are introduced, i.e., the cohesion and the friction angle, respectively,
obtaining:

$\tau_{\max}(t_{N},\|\boldsymbol{g}_{T}\|_{2})=c-t_{N}\tan\left(\varphi\left(\|\boldsymbol{g}_{T}\|_{2}\right)\right).$
(4)

In Eq. (4), the friction angle generally depends on the modulus of the
tangential component of the displacement jump, i.e., the slippage, so as to
simulate a slip-weakening frictional behavior. Since a quasi-static approach
is used, we can replace the tangential displacement rate
$\dot{\boldsymbol{g}}_{T}$ in Eq. (3b) with the incremental tangential
displacement $\Delta\boldsymbol{g}_{T}$ with respect to the previous time-step
value [50].

The fault surface $\Gamma_{f}$ can be split into three non-intersecting
portions. Each portion is characterized by a different operating mode allowed
for by the possible combinations of the previous conditions (2)-(3):

  * 1.

Stick: the surface is compressed ($t_{N}<0$) and the shear traction modulus
does not exceed the limiting value provided by the criterion in Eq. (4). The
displacement field is continuous across $\Gamma_{f}$;

  * 2.

Slip: the normal traction is still negative, but the surface is free to slip.
In this case, Coulomb’s equality holds [51]:

$\boldsymbol{t}_{T}^{*}=\tau_{\max}(t_{N},\|\boldsymbol{g}_{T}\|_{2})\frac{\Delta\boldsymbol{g}_{T}}{\|\Delta\boldsymbol{g}_{T}\|_{2}}.$
(5)

Only the normal component of the displacement field is continuous across
$\Gamma_{f}$;

  * 3.

Open: the normal traction is non-negative and the two contact surfaces
$\Gamma_{f}^{+}$ and $\Gamma_{f}^{-}$ are free to move on condition to avoid
compenetration. Hence, the displacement field across $\Gamma_{f}$ is
discontinuous and the traction on the fault vanishes, i.e.,
$\boldsymbol{t}=0$.

For additional details on the mathematical formulation, see [47, 48, 49] and
more recently [41, 52, 53].

According to Terzaghi’s principle, the total stress tensor in a saturated
porous medium can be decomposed as the sum of two contributions: the effective
stress tensor acting on the solid skeleton and a volumetric term depending on
the averaged fluid pressure $p$:

$\hat{\boldsymbol{\sigma}}=\begin{cases}\boldsymbol{\sigma}-\boldsymbol{1}p&\text{on
}\Gamma_{\sigma}\cup\Gamma_{f},\\\
\boldsymbol{\sigma}-\alpha\boldsymbol{1}p=\mathbb{C}:\boldsymbol{\varepsilon}-\alpha\boldsymbol{1}p&\text{in
}\Omega,\end{cases}$ (6)

where the fluid pressure is averaged by the saturation indices of the
different phases, $\alpha$ is the Biot coefficient, taking care of the ratio
between the grain and the porous matrix compressibility, and $\boldsymbol{1}$
the identity tensor of order $2$ [54]. The effective stress tensor
$\boldsymbol{\sigma}$ is called effective Terzaghi stress tensor on the domain
boundaries, while it is the effective Biot stress tensor inside the domain
itself. In Eq. (6), the constitutive relationship defining the effective Biot
stress tensor is introduced, where $\mathbb{C}$ is a fourth order elasticity
tensor, generally non-linear, and
$\boldsymbol{\varepsilon}=\nabla^{s}\boldsymbol{u}$ is the strain tensor, with
$\nabla^{s}=(\nabla+\nabla^{T})/2$ the symmetric gradient operator. The
mechanical constitutive law relates a strain variation in the porous medium to
an effective stress variation. Such a law can be described by a simple linear
elastic model (Hooke’s law), with constant or variable parameters, but also
more complex elasto-plastic rules with time-dependent contributions can be
introduced, e.g., a visco-elasto-plastic law. For more details on the
appropriate constitutive laws and their implementations, see [51, 55].

###  2.2 Mass balance equation

The mass conservation of the fluid species $\kappa$ reads [56, 57, 58, 59]:

$\frac{\partial}{\partial
t}(\rho^{\kappa})+\nabla\cdot\boldsymbol{F}^{\kappa}=q_{s}^{\kappa},$ (7)

where $\rho^{\kappa}$ and $\boldsymbol{F}^{\kappa}$ are the density and the
flux, respectively, of the fluid species $\kappa$. The density $\rho^{\kappa}$
represents the mass of $\kappa$ per unit of rock volume and can be written as:

$\rho^{\kappa}=\phi\sum_{\beta}S_{\beta}\rho_{\beta}\chi_{\beta}^{\kappa},$
(8)

with $\phi$ the porosity, $S_{\beta}$ the saturation of phase $\beta$, that
can be either liquid or gas, $\rho_{\beta}$ the density of phase $\beta$, and
$\chi_{\beta}^{\kappa}$ the mass fraction of component $\kappa$ in phase
$\beta$. Usually, assuming isothermal conditions the fluid density is a
function of pressure, but can also depends on other quantities, such as the
mass fraction according to some equation of state. Saturations and mass
fractions are constrained by the well-known conditions:

$\sum_{\beta}S_{\beta}=1\quad\text{and}\quad\sum_{\kappa}\chi_{\beta}^{\kappa}=1.$
(9)

The fluid flux of component $\kappa$ is the sum of the fluxes for each phase:

$\boldsymbol{F}^{\kappa}=\sum_{\beta}\chi_{\beta}^{\kappa}\boldsymbol{F}_{\beta},$
(10)

and each phase flux is described by Darcy’s law as:

$\boldsymbol{F}_{\beta}=\rho_{\beta}\boldsymbol{v}_{\beta}=-\rho_{\beta}\frac{\boldsymbol{k}\
k_{r,\beta}}{\mu_{\beta}}\left(\nabla
p_{\beta}-\rho_{\beta}\boldsymbol{g}\right),$ (11)

with $\mu_{f}$ and $k_{r,\beta}$ the viscosity and the relative permeability
of the phase, $\boldsymbol{k}$ the permeability tensor, and $p_{\beta}$ the
pressure in phase $\beta$.

According to [54, 60], the porosity update accounting for poro-elastic effects
can be expressed as:

$\phi=\phi_{0}+\alpha\varepsilon_{v}+\frac{\left(\alpha-\phi_{0}\right)\left(1-\alpha\right)}{K_{d}}\left(p-p_{0}\right),$
(12)

where $K_{d}$ is the drained bulk modulus,
$\varepsilon_{v}=\text{trace}\left(\boldsymbol{\varepsilon}\right)$ is the
volumetric strain, and $\phi_{0}$ and $p_{0}$ are the reference porosity and
fluid pressure respectively. In Eq. (12), $p$ is the averaged fluid pressure,
computed as:

$p=\sum_{\beta}S_{\beta}p_{\beta}.$ (13)

In oedometric conditions and a constant total stress state, from Terzaghi’s
principle (Eq. (6)) we have $d\sigma_{z}=\alpha\,dp$ and

$\frac{\partial}{\partial
t}\varepsilon_{v}=C_{m}\alpha\frac{\partial}{\partial t}p,$ (14)

with $C_{m}$ is the vertical uniaxial compressibility. In this case, the mass
balance is decoupled from the linear momentum balance and can be solved in
advance, providing a pressure field acting as an external body load for the
structural problem. Even though this assumption is not guaranteed for the
reservoir application considered in the present work, at the (large) space and
time scale of interest coupling is weak and a one-way coupled approach, where
Eqs. (7) are solved first for all the fluid species and the averaged pressure
of Eq. (13) is then introduced into Eq. (1), is fully warranted, see for
instance [44, 61, 62, 63, 64, 65].

The set of Eqs. (7) for obtaining the pressure field in the porous medium is
usually solved applying a finite volume method because it preserves the mass
conservation at the elemental level [66, 67, 68]. Nevertheless, also a finite
element or mixed finite element approach can be successfully used, e.g., [69,
70, 71, 72, 73]. In our analysis, the numerical simulation of the multiphase
flow in the porous matrix has been carried out by using Open Porous Media, an
open-source reservoir simulator based on a classical finite volume
discretization [74, 75]. As to the computation of the pressure field within
the network of faults, two strategies can be employed: either the domain
explicitly contains the faults as thin 3D cells, or the pressure is extended
to the faults from the surrounding 3D cells, according to some physical
treatment of the contact surfaces as inner boundaries. We elect to use the
latter approach, thus allowing to represent the fault at the macroscale as
zero-thickness lower dimensional elements. Generally speaking, the two
limiting cases that can be met in reality are sealing and non-sealing faults.
In the former case, the fault acts as an impermeable barrier and the pressure
change does not propagate from one side to the other of the contact surfaces.
In this situation, we can assume that the pressure variation in the fault is
null. On the contrary, in the latter case, the fault is fully permeable and
does not exhibit any resistance to the fluid flow. In this situation, we
assume the pressure variation in the fault to be equal to the arithmetic
average of the pressure computed on the two side cells.

###  2.3 Variational formulation and discretization

In this section, the variational formulation for the strong form of the linear
momentum balance in Eq. (1), equipped with the constraints of Eqs. (2)-(3), is
presented.

The weak form of the governing equations naturally produces a variational
inequality because of the frictional contact constraints [47]. In order to
avoid this difficulty, it is possible to reduce the original inequality to a
standard variational formulation by an active-set strategy and either a
penalty regularization or the introduction of Lagrange multipliers. We elect
to use the Lagrange multiplier technique, which can be computationally more
expensive, because new primary unknowns are introduced and the resulting
algebraic problem gains a saddle-point nature, but generally much more
accurate, robust and stable. Moreover, though generating saddle-point systems,
this formulation allows to produce a sequence of linear problems less sensible
to ill-conditioning issues [76]. From a physical viewpoint, Lagrange
multipliers represent the traction field on the fault surfaces, thus the
stress evaluation and the related reactivation risk becomes straightforward.

We emphasize that to retrieve the complete variational formulation for the
problem at hand, defining the right functional spaces and conditions, is far
beyond the purposes of this work and we refer the interested readers to more
specific researches, like Kikuchi and Oden [47] and Wohlmuth [50]. Here, we
briefly report the selected function spaces and the residual equations. The
notations $(\cdot,\cdot)_{\Omega}$ and $\langle\cdot,\cdot\rangle_{\Gamma}$
denote the $L^{2}$-inner product of functions in $\Omega$ (3D domain) and
$\Gamma$ (lower dimensional domain), respectively. Let
$\mathcal{V}=[H^{1}(\Omega)]^{3}$ be the Sobolev space of vector functions
whose first derivatives belong to $L^{2}(\Omega)$; let $\mathcal{M}$ be the
dual space of the trace space $\mathcal{W}=[H^{1/2}(\Omega)]^{3}$; and let
$\mathcal{M}(t_{N},\|\boldsymbol{g}_{T}\|)$ be its subspace such that

$\mathcal{M}(t_{N},\|\boldsymbol{g}_{T}\|)=\left\\{\boldsymbol{\mu}\in\mathcal{M}:\langle\boldsymbol{\mu},\boldsymbol{v}\rangle_{\Gamma_{f}}\leq\langle\tau_{\max}(t_{N},\|\boldsymbol{g}_{T}\|_{2}),\|\boldsymbol{v}_{T}\|\rangle_{\Gamma_{f}},\boldsymbol{v}\in\mathcal{W}\text{
with }v_{N}\leq 0\right\\}.$ (15)

Given the finite-dimensional subspaces $\mathcal{V}^{h}\subset\mathcal{V}$ and
$\mathcal{M}^{h}(t^{h}_{N},\|\boldsymbol{g}^{h}_{T}\|)\subset\mathcal{M}(t_{N},\|\boldsymbol{g}_{T}\|)$,
the finite dimensional weak form of the problem in Eq. (1) with Terzaghi
relation of Eq. (6) and conditions Eqs. (2)-(3) can be stated as follows: at
every instant $t\in[0,t_{\max}]$, find
$\\{\boldsymbol{u}^{h},\boldsymbol{t}^{h}\\}\in\mathcal{V}^{h}\times\mathcal{M}^{h}(t^{h}_{N},\|\boldsymbol{g}^{h}_{T}\|)$
such that:

$\displaystyle\mathcal{R}_{u}$
$\displaystyle=\left(\nabla^{s}\boldsymbol{\eta},\hat{\boldsymbol{\sigma}}\right)_{\Omega}-\langle\boldsymbol{\eta},\hat{\boldsymbol{\sigma}}\cdot\boldsymbol{n}\rangle_{\Gamma}-\left(\boldsymbol{\eta},\boldsymbol{b}\right)_{\Omega}$
$\displaystyle=\left(\nabla^{s}\boldsymbol{\eta},\hat{\boldsymbol{\sigma}}\right)_{\Omega}-\langle\boldsymbol{\eta},\hat{\boldsymbol{\sigma}}\cdot\boldsymbol{n}_{f}^{+}\rangle_{\Gamma_{f}^{+}}-\langle\boldsymbol{\eta},\hat{\boldsymbol{\sigma}}\cdot\boldsymbol{n}_{f}^{-}\rangle_{\Gamma_{f}^{-}}-\langle\boldsymbol{\eta},\hat{\boldsymbol{\sigma}}\cdot\boldsymbol{n}\rangle_{\Gamma_{\sigma}}-\left(\boldsymbol{\eta},\boldsymbol{b}\right)_{\Omega}$
$\displaystyle=\left(\nabla^{s}\boldsymbol{\eta},\boldsymbol{\sigma}(\boldsymbol{u}^{h})-\alpha\boldsymbol{1}p\right)_{\Omega}-\langle\llbracket\boldsymbol{\eta}\rrbracket,\boldsymbol{t}^{h}-p\boldsymbol{n}_{f}\rangle_{\Gamma_{f}}-\langle\boldsymbol{\eta},\bar{\boldsymbol{t}}\rangle_{\Gamma_{\sigma}}-\left(\boldsymbol{\eta},\boldsymbol{b}\right)_{\Omega}=0$
$\displaystyle\forall\boldsymbol{\eta}\in\mathcal{V}^{h},$ (16a)
$\displaystyle\mathcal{R}_{t}$ $\displaystyle=\langle
t_{N}^{h}-\mu,g_{N}\rangle_{\Gamma_{f}}+\langle\boldsymbol{t}_{T}^{h}-\boldsymbol{\mu}_{T},\Delta\boldsymbol{g}_{T}\rangle_{\Gamma_{f}}\geq
0$
$\displaystyle\forall\boldsymbol{\mu}\in\mathcal{M}^{h}(t^{h}_{N},\|\boldsymbol{g}^{h}_{T}\|).$
(16b)

We use a Galerkin approach, hence the test functions $\boldsymbol{\eta}$ and
$\boldsymbol{\mu}$ belong to the same function spaces used to define the trial
functions for the displacement and traction field, respectively. To transform
the variational inequality of Eq. (16b) into a variational equality, an
iterative active-set algorithm [77, 78] is applied. According to this
approach, the fault surfaces $\Gamma_{f}$ is subdivided into active and
inactive regions for both components of the traction, i.e., the Stick
($\Gamma_{f}^{\text{stick}}$, active for normal and tangential components),
Slip ($\Gamma_{f}^{\text{slip}}$, active for normal component), and Open
($\Gamma_{f}^{\text{open}}$, inactive) portions of $\Gamma_{f}$. With this
subdivision, the variational inequality of Eq. (16b) becomes:

$\displaystyle\mathcal{R}_{t}=\langle\boldsymbol{\mu},\llbracket\boldsymbol{u}^{h}\rrbracket\rangle_{\Gamma^{\text{stick}}_{f}}+\langle\mu_{N},g_{N}\rangle_{\Gamma^{\text{slip}}_{f}}+\frac{1}{k}\langle\boldsymbol{\mu}_{T},\boldsymbol{t}_{T}^{h}-\boldsymbol{t}_{T}^{*}\rangle_{\Gamma^{\text{slip}}_{f}}+\frac{1}{k}\langle\boldsymbol{\mu},\boldsymbol{t}^{h}\rangle_{\Gamma^{\text{open}}_{f}}$
$\displaystyle=0$
$\displaystyle\forall\boldsymbol{\mu}\in\mathcal{M}^{h}(t^{h}_{N},\|\boldsymbol{g}^{h}_{T}\|),$
(17)

where $k$ is a unitary coefficient introduced to ensure dimensional
consistency. The non-linear system of Eqs. (16a)-(17) is solved by a Newton
linearization and at convergence we check the consistency of the traction
state on the faults with the initial subdivision of $\Gamma_{f}$ into
$\Gamma_{f}^{\text{stick}}$, $\Gamma_{f}^{\text{slip}}$, and
$\Gamma_{f}^{\text{open}}$. If the consistency check is satisfied, the active-
set algorithm is stopped and we can move to the following time instant,
otherwise a new subdivision of $\Gamma_{f}$ is defined, Eq. (17) is re-
computed and the resulting non-linear system solved again.

Introducing $\boldsymbol{u}^{h}=\sum_{i}u_{i}\boldsymbol{\eta}_{i}$ and
$\boldsymbol{t}^{h}=\sum_{j}t_{j}\boldsymbol{\mu}_{j}$, i.e., the discrete
representation of the displacement and traction fields, where
$\\{\boldsymbol{\eta}_{i}\\}$ and $\\{\boldsymbol{\mu}_{j}\\}$ are bases for
$\mathcal{V}^{h}$ and $\mathcal{M}^{h}(t^{h}_{N},\|\boldsymbol{g}^{h}_{T}\|)$,
respectively, the set of variational equalities of Eqs. (16a)-(17) becomes an
algebraic nonlinear system. The bases for the finite-dimensional spaces
$\mathcal{V}^{h}$ and $\mathcal{M}^{h}$ are selected with the aid of the
finite element method. Given the regularity requirements defined above, we use
low-order discretization spaces for both displacement and traction. The
computational domain is subdivided into non-overlapping hexahedral elements,
$\Omega=\bigcup_{i=1}^{n_{e}}\Omega_{i}$ and
$\Omega_{i}\cap\Omega_{j}=\varnothing$ for any $i\neq j$. This choice is done
for the sake of the consistency with the domain discretization used for the
multiphase flow model, which is based on a standard finite volume approach. We
use a conformal representation of the faults, i.e.,
$\Gamma_{f}=\bigcup_{j}\partial\Omega_{j}$, where $\Omega_{j}$ are elements
sharing a face with $\Gamma_{f}$. In such a way, the fault contact surfaces
are composed by pairs of quadrilateral elements. Each pair of quadrilateral
elements is also denoted as a zero-thickness interface finite element [79].
According to the value of the traction, every interface element can change its
status, i.e., it can belong to either the stick, slip or open portion of
$\Gamma_{f}$.

The mixed finite element discretization adopted in this work is described in
Franceschini et al. [52, 53]. It consists of a $\mathbb{Q}_{1}$ first-order
interpolation for the nodal-based displacement field and a $\mathbb{P}_{0}$
piecewise constant interpolation for the element-based traction field. The
collection of coefficients $u_{i}$ and $t_{j}$ are the components of the
unknown algebraic arrays, named $\mathbf{u}$ and $\mathbf{t}$ of sizes
$3n_{n}$ and $3n_{f}$, with $n_{n}$ the number of nodes in the hexahedral mesh
and $n_{f}$ the number of interface elements. This approach has the main
advantage of being naturally coupled with a finite volume pressure solution
computed on the same domain discretization at no additional cost, since both
the traction and the pressure are represented using the same space on the same
grid. On the other side, in order to ensure the LBB-stability of the proposed
mixed finite element spaces, a tailored jump stabilization has been proposed
in [53]. For a complete analysis of the LBB-stability of a general pair of
mixed finite element spaces, see Elman et al. [80]. An implementation of the
presented algorithm can be found in [81].

We emphasize that even if a linear elastic constitutive relation is used for
the porous medium, the set of equations reported in Eqs. (16a)-(17) represents
a non linear problem, because a consistent partitioning of the fracture
surface is unknown and has to be computed depending on the solution vectors.
To be more specific, constraints in Eqs. (2)-(3) are the Karush-Kuhn-Tucker
(KKT) conditions and we are dealing with a non-linear optimization problem
[77]. As already mentioned, the solution strategy is based on an active-set
strategy, with each non-linear problem addressed by a classical exact Newton
algorithm.

###  2.4 Linearization and linear system solution

The use of a mixed finite element approximation produces a linearized step
with a generalized saddle-point Jacobian matrix [82]. The Jacobian is
generally non-symmetric because of the contribution related to the friction
component of the traction when the fracture slides. In particular, at each
Newton iteration, the linear system that has to be solved is:

$J\delta\mathbf{x}=-\mathbf{r},$ (18)

with the $2\times 2$ block matrix $J$, the residual vector $\mathbf{r}$, and
the solution vector $\delta\mathbf{x}$ given by:

$J=\begin{bmatrix}\frac{\partial\mathcal{R}_{u}}{\partial\boldsymbol{u}}&\frac{\partial\mathcal{R}_{u}}{\partial\boldsymbol{t}}\\\
\frac{\partial\mathcal{R}_{t}}{\partial\boldsymbol{u}}&\frac{\partial\mathcal{R}_{t}}{\partial\boldsymbol{t}}\end{bmatrix},\quad\mathbf{r}=\begin{bmatrix}\mathcal{R}_{u}\\\
\mathcal{R}_{t}\end{bmatrix},\quad\text{and}\quad\delta\mathbf{x}=\begin{bmatrix}\delta\mathbf{u}\\\
\delta\mathbf{t}\end{bmatrix},$ (19)

where $\mathcal{R}_{u}$ and $\mathcal{R}_{t}$ are computed at the current
counter level $l$. As usual, the updated solution vector at level $l+1$ is:

$\begin{bmatrix}\mathbf{u}\\\ \mathbf{t}\\\
\end{bmatrix}^{l+1}=\begin{bmatrix}\mathbf{u}\\\ \mathbf{t}\\\
\end{bmatrix}^{l}+\begin{bmatrix}\delta\mathbf{u}\\\ \delta\mathbf{t}\\\
\end{bmatrix}.$ (20)

For the detailed expression of the Jacobian, see [52]. The resulting linear
system is characterized by a large and sparse matrix and it is necessary to
use a preconditioned iterative method for its efficient solution. To achieve
satisfactory results, the use of a suitable preconditioner is mandatory. Since
the properties of the linear system may change significantly as the simulation
proceeds and the fault elements change status, also the preconditioner must
evolve. Among others, an idea is to exploit the scalability intrinsically
present in the multigrid approach and combine it with the known physics-based
partitioning of the blocks to be able to solve the saddle-point matrix. For
details on robust and efficient techniques used for the solution of this
peculiar linear system, the reader may refer to [42, 53, 83].

###  2.5 Constitutive model for fracture: slip weakening

In this work we use both the classical Coulomb criterion with a constant
friction coefficient, and a slip-weakening friction law with a variable
friction coefficient. Originally used by Andrews [84] to take into account the
change from static to dynamic friction, slip-weakening friction laws [85, 86]
are based on the concept that the shear stiffness of the fracture decreases as
sliding occurs. From a mathematical viewpoint, the standard Coulomb criterion
reads:

$\left\|\boldsymbol{t}_{T}\right\|_{2}\leq c-t_{N}\mu,\quad\text{with
}\mu=\tan\varphi,$ (21)

while a more general slip-weakening friction law reads:

$\left\|\boldsymbol{t}_{T}\right\|_{2}\leq
c-t_{N}\mu(\|\boldsymbol{g}_{T}\|_{2}).$ (22)

A simple expression to account for the friction reduction with fault motion is
provided by a piecewise linear function, as shown in Fig. 3, where the
friction coefficient linearly decreases from the static value $\mu_{s}$ down
to the dynamic value $\mu_{d}$ at a sliding value equal to $D_{c}$. For larger
sliding values, the friction coefficient remains constantly equal to
$\mu_{d}$. Other analytical expression can be also used to simulate the
friction coefficient reduction with the sliding, such as an exponential law
(see Fig. 3), which has the advantage to allow for a smooth variation that can
be differentiated everywhere. It provides:

$\mu=\mu_{d}+\left(\mu_{s}-\mu_{d}\right)\exp\left(-\frac{\left\|\boldsymbol{g}_{T}\right\|_{2}}{D_{c}}\right).$
(23)

A similar smooth behavior can be formulated based on inverse trigonometric
functions (see Fig. 3):

$\mu=\mu_{d}+\left(\mu_{s}-\mu_{d}\right)\left(1-\frac{2}{\pi}\arctan\frac{\left\|\boldsymbol{g}_{T}\right\|_{2}}{D_{c}}\right).$
(24)

Figure 3: Slip weakening constitutive laws. From left to right: piecewise
linear friction law, exponential friction law and inverse trigonometric
friction law.

In order to compare the three different expressions, we use the simple 1D
problem sketched in Fig. 4. The selected physical parameter set is:
$\mu_{s}=\tan(30^{\circ})$, $\mu_{d}=\tan(10^{\circ})$, $D_{c}=2$ mm, with the
spring stiffness $K=11\times 10^{9}$ N/m and a compression load $N=3\times
10^{7}$ N. The first three values are representative of the conditions
typically found in the seismogenic gas fields within the Rotliegend
stratigraphic units in the Netherlands [45, 87]. The physical quantities of
interest are shown in Fig. 4. The primary variable is always the displacement
of the point where the external load $N$ is applied, while the outcomes are:
(i) the friction strength $F$, (ii) the relative displacement $u_{r}$ between
the body connected to the spring and the fixed basement, (iii) the global
system stiffness $\overline{K}$, and (iv) the internal energy $U$. Though the
response in terms of friction strength are different, both relative
displacement and global energy are comparable. By distinction, the global
stiffness behaves differently and for two cases out of three it reaches
negative amounts that are greater in absolute value than the original spring
stiffness $K$. The finite element approach used in the present modeling
analysis is based on the global equilibrium of the system and not on a local
(elemental) balance. This is the reason why a comparison of the global energy
associated with the different laws is meaningful. At the elemental level, it
is desirable to avoid negative stiffness, which could potentially lead to
friction instabilities. Hence, we chose to work with the law based on the
inverse trigonometric function, i.e., the only providing a minimum stiffness
smaller in absolute value than the original one. In such a way, we can ensure,
at least for conditions similar to the ones used in this example, a positive
global stiffness.

Figure 4: Mono-dimensional friction system with 1 degree of freedom. On the
left: sketch of the model used in this example. On the right, from the top
left, (i) friction strength, (ii) relative displacement, (iii) global system
stiffness and (iv) internal energy. Continuous, dotted and dashed line
represent linear, exponential and inverse trigonometric slip-weakening
formulation, respectively

##  3 Model set-up

###  3.1 Conceptual model

In order to test the mathematical model and identify the main mechanisms
governing fault reactivation in UGS fields, we use a simplified geological
model representative of the typical features of the Rotliegend UGS reservoirs,
such as Norg and Grikpskerk [18], e.g., see Fig. 5. These reservoirs are
bounded by normal faults with a significant throw (up to a 250 m) and consist
of a few compartments separated by internal faults. The gas fields are located
between 2000 and 3000 m of depth, with the Rotliegend reservoir rock
characterized by an average net thickness of 150-200 m. Detailed information
about the geological setting and typical geometric features of UGS reservoirs
in the Netherlands can be found, among others, in [16, 37, 88, 89] and
published reports [17, 18].

Figure 5: On the left: base Zechstein semblance map of the Norg UGS (in blue)
and surrounding area with traces of the bounding faults and localization of
the recorded seismic events. On the right: conceptual map of the Norg field
with major and minor faults highlighted in blue and red, respectively.

Based on those features, we define a conceptual model composed of two adjacent
compartments, 2000$\times$2000 m wide, 200 m thick, and 2000 m deep, where UGS
activities are carried out. The reservoir compartments are laterally confined
by two families of orthogonal faults, denoted as F1-F2 (parallel to $y$-axis)
and F4-F5 (parallel to $x$-axis). Another fault, denoted as F3, separates the
two reservoir blocks (Fig. 6). The two compartments have only a partial
hydraulic connection depending on the sealing properties of fault F3, so the
pore pressure distribution in space and time may be different. Faults F1 and
F2 are inclined with respect to the vertical $z$-axis by a dip angle equal to
$\pm 10^{\circ}$, while F3, F4 and F5 are vertical faults, as shown in Fig. 6.
The faults extend from -3000 m to -1600 m depth, i.e., they terminate within
the caprock sealing the reservoir, called Zeichestein formation. Notice that
the blocks have a 200-m offset along the vertical direction, corresponding to
the entire thickness of the reservoir, relative to the Rotliegend formation
located in the sideburden.

Figure 6: On the left: plain view of the model. On the right: vertical
sections of the conceptual model along the trace A-A and B-B shown on the
left.

###  3.2 Finite element-interface element discretization

The reservoir is embedded in a 30-km wide square domain. The overall model
size is much larger than the reservoir dimension to minimize the effects of
the (arbitrary) boundary conditions on the solution in the area of interest
(Fig. 6). The bottom of the model is 5000-m deep and the land surface is
located at the elevation of 0 m. Standard conditions with zero displacement
and zero pore pressure variation on the outer and bottom boundaries are
prescribed, whereas the land surface is a traction-free boundary.

A 3D finite element mesh of the selected domain is built by using hexahedral
elements, which are particularly suitable for the symmetric configuration with
the faults parallel to the Cartesian axes. Fig. 7 shows an axonometric view of
the full computational grid used in the geomechanical model. The mesh consists
of 253,165 nodes and 236,208 hexahedral elements with a finer discretization
in the reservoir layers, i.e., at depth between 2000 and 2200 m. The element
size within the reservoir is 100$\times$100$\times$20 m. Fig. 8 shows the
fault system embedded in the continuous 3D grid as discretized by 5,215
interface elements. The state of each element of the faults is synthetically
evaluated with the aid of the criticality index defined as:

$\chi=\frac{\|\boldsymbol{t}_{T}\|_{2}}{\tau_{\max}}=\frac{\|\boldsymbol{t}_{T}\|_{2}}{c-t_{N}\tan\left(\varphi\left(\|\boldsymbol{g}_{T}\|_{2}\right)\right)}.$
(25)

From Eq. (25), it is easy to see that $\chi\in[0,1]$, where 0 is associated
with the safest condition and 1 to plastic sliding.

Figure 7: Axonometric view of the computational domain used for the
geomechanical simulations: full 3D finite element grid (left) and interface
element grid (blue) embedded in a portion of the full 3D grid (right). Figure
8: Interfece element discretization of the fault discontinuities. The planar
trace of F1, F2 and F3 is parallel to the y-axis, whereas that of F4 and F5 is
parallel to the $x$-axis. F3 is the central fault separating the two reservoir
compartments.

###  3.3 Simulated scenarios

To evaluate the capabilities of the presented numerical model and understand
the possible mechanisms causing fault reactivation during CGI and UGS, a few
scenarios are simulated in the typical setting of the Rotliegend reservoirs in
the Netherlands. The main geological and geomechanical parameters are reported
in Tab. 1. For the sake of simplicity, a linear elastic behavior is assumed in
the reservoir during the UGS activities.

The pressure history prescribed in an active well located in each compartment
is sketched in the leftmost frame of Fig. 9. We assume a 10-y duration for the
PP phase, where the pressure drops linearly by up to 20 MPa. After this
period, a 2-year CGI phase follows, where the pressure recovers to the initial
(undisturbed) value $p_{i}$, and then UGS cycles start. They are characterized
by a 6-month extraction period, during which the pressure drops by 10 MPa, and
a 6-month injection period, when the pressure returns to $p_{i}$.

Table 1: Formation-dependent geomechanical parameters. See Fig. 6 for a detail on the depths. layer | density [kg/m3] | Young modulus [GPa] | Poisson ratio
---|---|---|---
Overburden | 2200 | 10.0 | 0.25
Upper Zechstein Salt (-1500 to -1800 m) | 2100 | 35.0 | 0.30
Lower Zechstein Salt (below -1800 m) | 2100 | 20.0 | 0.30
Reservoir (Upper Rotliegend) | 2400 | 11.0 | 0.15
Underburden | 2600 | 30.0 | 0.20

Figure 9: On the left: sketch of the pore pressure variation in time
prescribed in an active well of the reservoir compartments. On the right:
initial normal stress with respect to the fault orientation. The principal
stresses $\sigma_{h}$, $\sigma_{H}$ and $\sigma_{v}$ are parallel to the
Cartesian axes. Faults F4 and F5 are more loaded because of their
orthogonality to $\sigma_{H}$.

To initialize the simulation, the undisturbed stress regime must be
prescribed. We assume it has the principal effective stress tensor directions
aligned with the Cartesian axes, in particular,
$\sigma_{1}=\sigma_{v}=\sigma_{z}$, $\sigma_{2}=\sigma_{H}=\sigma_{y}$, and
$\sigma_{3}=\sigma_{h}=\sigma_{x}$, where $\sigma_{v}$ denotes the vertical
compressive stress, and $\sigma_{H}$ and $\sigma_{h}$ the largest and smallest
compressive horizontal principal stress, respectively. At the reservoir
average depth, i.e., $z=-2100$ m, we have $\sigma_{v}=-25.4$ MPa,
$\sigma_{h}=M_{1}\sigma_{v}=-18.8$ MPa, $\sigma_{H}=M_{2}\sigma_{v}=-21.1$
MPa, with $M_{1}=0.40$ and $M_{2}=0.47$. The initial normal stress acting on
the faults is shown in Fig. 9.

Two scenarios have been simulated based on the parameters describing the
Coulomb frictional criterion. In the reference scenario (scenario 1)
$\varphi_{s}=30^{\circ}$ and fault weakening is not accounted for. The effect
of slip-weakening behavior is investigated in scenario 2, where the friction
angle reduces from $\varphi_{s}=30^{\circ}$ to $\varphi_{d}=10^{\circ}$ in a
slip distance of $D_{c}=2$ mm. Cohesion $c=2$ MPa in both scenarios. Finally,
the time step is 1 year during the PP phase, then during CGI and UGS phases it
is reduced to 2 months.

##  4 Numerical results

The objective of the representative simulations reported herein is to evaluate
the fault reactivation risk during the different stages of the UGS activities
in the conceptual reservoir. For this reason, we mainly focus on the
criticality index $\chi$ defined in Eq. (25). For the sake of clarity and ease
of readability, $\chi$ is represented for each fault as a function of depth
only, i.e., for each $z$-value we compute the $\chi$ average for the stripe of
interface elements located at the same depth. Another significant quantity is
the maximum sliding, i.e., the maximum value of $\|\boldsymbol{g}_{T}\|_{2}$
simulated along each fault. These two quantities are closely related each
other, since a single element can slide only when $\chi=1$. However, we prefer
to propose an averaged version of $\chi$, so as to obtain information on the
criticality state of the entire fracture at a given depth.

The last quantity used to interpret the results and analyze the fault behavior
is the tangential component of the traction. In particular, we use $t_{T,z}$,
i.e., the vertical component of $\boldsymbol{t}_{T}$. Usually, the 2-norm of
the tangential traction is analyzed, i.e., $\|\boldsymbol{t}_{T}\|_{2}$, but
this does not provide information on the shear direction. However, thanks to
the symmetric geometry of the conceptual model, in some locations there is no
horizontal component of $\boldsymbol{t}_{T}$, thus,
$\|\boldsymbol{t}_{T}\|_{2}=|t_{T,z}|$. The two quantities share the same
modulus, but the vertical component carries additional information on the
sliding direction.

Figure 10: Time behavior of the production/injection rate used in the OPM
simulation.

###  4.1 Pore pressure variation

As previously mentioned, in this work we adopt a one-way coupled approach,
thus the multiphase flow prediction is computed first. The simulation is
performed through the open-source reservoir simulator Open Porous Media [74,
75]. As a reference scenario, a typical year-long cycle of UGS activity has
been considered, with the injection-production history represented in Fig. 10.
Fig. 11 shows the location of the injection/production wells with respect to
the fault system. Note that, to avoid any interpolation among computational
grids, the OPM finite volume mesh exactly corresponds to the finite element
grid of a single block within the 3D geomechanical model. The characteristic
horizontal and vertical permeability of a reservoir in the study area are
$k_{h}=600$ mD and $k_{v}=300$ mD, respectively. The “working gas” volume
amounts to 6.5$\times$109 Sm3 per compartment.

Figure 11: Location of the injection/production wells in the two reservoir
blocks (left) and axonometric view of the 3D computational grid used in OPM to
simulate the injection/production phase in each reservoir compartment (right).
The OPM mesh exactly corresponds to the finite element grid of a single block
within the 3D geomechanical model.

The numerical results in terms of pressure variation are summarized in Fig.
12. The figure shows the depth-averaged pressure behavior along a section
passing through the production/injection wells every 3 months. After 3 months
the maximum production rate is achieved, after 6 months the production phase
ends, after 9 months the maximum injection rate is met, and, finally, after 12
months the simulation ends. Notice that the pressure perturbation during the
entire production (or injection) phase is almost uniform in space and varies
approximately within the interval between 0 and -10 MPa with respect to the
initial value $p_{i}$. This outcome shows that, for the setting defined in
these representative simulations, the spatial gradient of the pore pressure
variation into each compartment is expected to be quite limited. Hence,
considering a constant pressure variation value for each reservoir block
appears also to be a reasonable assumption.

Figure 12: Depth-averaged value of the pore pressure variation during a
production/injection cycle as obtained by the OPM flow simulator.

###  4.2 Analysis of the fault reactivation risk

The mechanisms for the possible fault reactivation have been investigated in
scenario 1. The value $\chi_{\max}=1$ is reached on faults F1 and F2 at
loading step 9, with $\chi_{\max}$ up to 0.8 at the end of CG and UGS
injection phases (Fig. 13). Conversely, $\chi=0$ on fault F3 irrespective of
the loading step due to the symmetry of the geometry and loading
configurations. A comparison between the behavior versus depth of the
criticality index along fault F1 and F4 and the distribution of $\chi$ on the
whole fault system at the end of PP are shown in Fig. 14. Notice that the most
critical condition develops along the top and bottom of the reservoir in
agreement with previous modeling study [35]. Moreover, faults F4 and F5
exhibit smaller values of $\chi$ with respect to F1 and F2, showing that a
sub-vertical orientation is usually more likely to reactivate.

Figure 13: Behavior of $\chi_{\max}$ from all the loading steps for each
fault. Note that due to symmetry F1 and F2 behave identically, as well as F4
and F5.

Figure 14: On the left: behavior of the criticality index $\chi$ vs depth at
loading step 10 (end of PP) for faults F1 and F4. On the right: $\chi$ factor
on all the fracture surfaces at loading step 10.

Figure 15: On the left: location of the selected element. On the right: stress
path $\|\boldsymbol{t}_{T}\|_{2}$ vs $-t_{N}$ for the element highlighted on
the left sketch by a red dot. The red line is the yield bound. Numbers along
the path denote the loading steps. It can be easily recognized the primary
production (loading steps 1 to 10), the cushion gas injection (loading steps
10 to 12) and the underground gas storage (loading steps 12 to 12.5 –
production – and 12.5 to 13 – injection).

Fig. 15 shows the stress path in the $t_{N}-\|\boldsymbol{t}_{T}\|_{2}$ plane
experienced by a representative element located on fault F1 at the top of the
reservoir. The actual stress state touches the yield bound at loading step
number 9 and remains on the yield surface till the end of PP (loading step
10). During CGI, the stress state initially departs from the yield condition
but returns close to it during the last part of the injection when the
pressure recovers to the initial value. UGS behaves elastically over a new
path with respect to what experienced during the last part of the CGI phase,
with an almost constant $t_{N}$ value. Again, the stress state approaches the
critical condition at the end of the UGS injection phase when $p$ rises back
to $p_{i}$.

Figure 16: On the top left: distribution of vertical component of the shear
stress $t_{T,z}$ for the loading steps (l.s.) 0 (initial condition), 10, 11,
12, 12.5, 13 on fault F1 (dip = 10∘). On the top right: time behavior of
$\|\boldsymbol{t}_{T}\|_{2}$ for the points denoted by the thick black dots in
the previous frame located at the top, bottom, and center of the reservoir.
Positive values mean that the shear stress is directed upward. On the bottom:
sketches representing the reservoir deformation, shear stress direction, and
inactive/active portions of the fault at the same loading steps.

A deeper explanation for this behavior can be found by analyzing the actual
direction of the shear stress. Fig. 16 shows the vertical component of
tangential traction $t_{T,z}$ on fault F1 at loading steps 0 (initial
condition), 10, 11, 12, 12.5, and 13. This component is meaningful because of
the symmetry of the model, indeed, we have that
$\|\boldsymbol{t}_{T}\|_{2}=|t_{T,z}|$. Sketches of the reservoir-fault-
sideburden conditions are provided for the same loading steps. The initial
shear stress differs from the null value because of the fault dip. The largest
value of $t_{T,z}$ are observed at the end of PP (loading step 10). Reservoir
compaction induced by the pressure depletion is accompanied by fault
reactivation. Note that a positive and negative shear stress characterizes the
reservoir bottom and top, respectively. As physically expected because of the
compaction mechanism, the direction of the shear stress is oriented toward the
center of the reservoir. When CGI starts, the shear stress orientation changes
and the reactivated part of the fault returns stick. At loading step 11, half
of the pore pressure change has been recovered. As the reservoir expands due
to pressure recovery, $t_{T,z}$ decreases on the reservoir top and bottom (the
orientation remains the same but the absolute value decreases) and an almost
null $t_{T,z}$ on the previously sliding IEs is obtained at this step.
Differently, $t_{T,z}$ does not significantly change for the elements
surrounding the activated stripes of the fault. The reservoir continues to
recover pressure and re-expand until loading step 12. During this second part
of CGI shear stress increases, with a sign opposite to that experienced during
PP (Fig. 16). A mirror behavior occurs for the IEs at the reservoir bottom.
Therefore, expansion during CGI increases the criticality condition of the
fault (mainly at the reservoir top and bottom) due the stress re-distribution
after the sliding developed over the PP. Fig. 13 shows that faults F1 and F2
approach the criticality state ($\chi_{\max}>0.8$) when the pressure recover
the initial value at the end of CGI and UGS injection phase, i.e., in a
pressure state close to the initial undisturbed one, which is not generally
expected to be associated to fault reactivation.

Figure 17: Maximum sliding versus time for the investigated scenarios. On the
left: reference case (scenario 1). On the right: using slip-weakening
constitutive law (scenario 2).

###  4.3 Slip-weakening effect

The adopted Coulomb frictional criterion can handle slip-weakening effects.
Here the outcome of a slip-weakening constitutive law for the fault behavior
is compared to that previously obtained using a static friction coefficient
equal to $\varphi_{s}=30^{\circ}$. The two parameters defining the new
constitutive law are $\varphi_{d}$ and $D_{c}$, i.e., the dynamic friction
angle and the slip weakening distance, respectively. In the simulated
scenario, the friction angle reduces from $\varphi_{s}=30^{\circ}$ to
$\varphi_{d}=10^{\circ}$ in a slip distance of $D_{c}=2$ mm.

Figure 18: Effect of the Coulomb parameters on $\chi_{\max}$ at increasing
loading steps for each fault. As usual, the pairs F1-F2 and F4-F5 behave
identically due to symmetry. The proposed scenario corresponds to
$\varphi_{d}=10^{\circ}$ and $D_{c}=2$ mm. Figure 19: Zoom of Fig. 18 over the
cushion gas injection and UGS phases for faults F1 = F2 and F4 = F5.

Fig. 17 provides the time behavior of the fault maximum sliding for the
proposed scenario. It can be seen that the current sliding is more than twice
that obtained using a static friction angle. Fig. 18 shows a comparison
between the criticality index during the entire simulation for scenario 1 and
2. It can be noticed that the new constitutive law causes F1, F2, F4 and F5 to
slip as well at the end of the cushion gas and UGS injection phases, but not
at loading step 12.5, i.e., at the end of the 6-month UGS production phase
(see the zoom in Fig. 19).

Finally, Fig. 20 shows the stress path for the same location as in Fig. 15.
Because of the reduced friction angle, the yield surface is reached more
easily during PP, at the end of CGI, and at the end of UGS phases. As observed
for the reference scenario, the elastic phases develop with an almost constant
normal stress because of the selected ratios between the reservoir and
overburden stiffness and between the pressure change in the reservoir and
within the fault. The stress path and the yield bound are quite complex due to
weakening. Moreover, due to the very small friction angle
($\varphi=10^{\circ}$), a large part of UGS is characterized by a the stress
state that develops either on the yield surface or very close to it.

Figure 20: Stress path $\|\boldsymbol{t}_{T}\|_{2}$ vs $-t_{N}$ for the F1
element highlighted in Fig. 15. The dashed and the continuous red lines are
the yield bound corresponding to the static condition ($\varphi_{s}$) and
after the slip distance $D_{c}$ is overcome, respectively. The numbers along
the path denote the loading steps. The primary production (loading steps 1 to
10), the cushion gas injection (loading steps 11 to 12) and the underground
gas storage (loading steps 12 to 12.5 – production – and 12.5 to 13 –
injection) can be easily recognized.

##  5 Conclusions

The first underground gas storage site became operational as early as 1915
[90]. Since then, this technology has spread to all continents, reaching
nowadays more than 600 facilities worldwide. Despite this large use, there are
risks related to the possible reactivation of existing geological fractures.
Although it is a “rare” event from a statistical viewpoint [28], it deserves
the proper attention due to its strong social and economic effects. Most of
recorded human-induced seismic events can be explained by a pressure increase
until it exceeds the initial value, triggering the shear stress on the fault
surface to reach the limit strength. However, there are recorded events that
cannot be explained by this mechanism. They are the so-called “unexpected”
seismic events, which occur when the pressure is in the range already
experienced during the primary production. The main scope of this work is to
identify these phenomena, explain their basic processes, and define some safe
operational bandwidth for the UGS reservoir management with reference to the
gas fields located in the Roetligend formation, the Netherlands. To accomplish
these aims, we use a computational modeling approach for the accurate and
robust simulation of the mechanics of faulted porous rocks.

The overall work is split into two parts. This paper deals with Part I, which
concerns the development, implementation and test of the mathematical and
numerical model used for computational simulations. A one-way coupled strategy
is adopted to deal with the poro-mechanical interaction. First, the set of
governing relationships for frictional contact mechanics are introduced, then
the weak variational formulation is derived and discretized. We use Lagrange
multipliers to prescribe the normal and frictional constraints on faults, a
mixed-dimensional approach, and a mixed finite element discretization with
displacement in the 3D porous body and traction on the fault surfaces are the
main unknowns. In order to be consistent with classical finite volume
discretizations for the multiphase flow, we focus on low-order hexehedral
elements for the 3D continuum and a piecewise constant representation of the
traction on the contact surfaces, thus requiring a proper stabilization to
ensure the regularity of the resulting generalized saddle-point problem. An
active-set algorithm and an exact Newton method are implemented for the
solution of the overall non-linear problem, while ad hoc preconditioning
strategies are used to allow for and accelerate the convergence of the inner
linear Krylov solver. A discussion on the slip-weakening constitutive law for
fault frictional behavior is also provided. Finally, the model is applied to
two realistic scenarios, carried out on a conceptual model built from an
idealization of real UGS fields located in the formation of interest. Modeling
simulations allow to identify the main mechanisms potentially inducing a fault
reactivation during UGS activities, even in “unexpected” situations where the
current stress state appears to be less demanding than what the porous medium
had already experienced in the past. The use of a slip-weakening rheological
model for the frictional behavior can increase the chance of producing a fault
reactivation during CGI and UGS activities.

Part II of this work will focus on the model application in a real-world
scenario, with an extensive sensitivity analysis on the factors that can
mostly impact on the reactivation chances. Further development concerns
widening the feasible parameter ranges, e.g., testing different constitutive
laws for the continuous medium, using different parameter values, changing the
fault positions and orientation. The analysis will be extended to other kinds
of storage activities, such as CO2 geological sequestration or underground H2
and N2 storage. The final objective is to draw some guidelines to define a
safe operational bandwidth for the management of storage reservoirs in the
Netherlands, and, at the same time, build a methodological example that can be
successfully extended to other real-world experiences.

## CRediT authorship contribution statement

Andrea Franceschini: Methodology, Software, Writing - Original Draft,
Visualization. Claudia Zoccarato: Conceptualization, Methodology, Formal
analysis, Investigation. Selena Baldan: Investigation, Visualization. Matteo
Frigo: Software, Investigation. Massimiliano Ferronato: Conceptualization,
Methodology, Writing - Review and Editing, Supervision. Carlo Janna: Software.
Giovanni Isotton: Software. Pietro Teatini: Conceptualization, Methodology,
Formal analysis, Writing - Review and Editing, Supervision, Funding
acquisition.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgements

This research was supported by the State Supervision of Mines (SodM), Ministry
of Economic Affairs (The Netherlands), Project KEM01 “Safe Operational
Bandwidth of Gas Storage Reservoirs” grant. Portions of this work were
performed within the 2020 INdAM-GNCS project “Optimization and advanced linear
algebra for PDE-governed problems”. Computational resources were provided by
University of Padova Strategic Research Infrastructure Grant 2017: “CAPRI:
Calcolo ad Alte Prestazioni per la Ricerca e l’Innovazione”.

## References

  * Doglioni [2018] C. Doglioni,  A classification of induced seismicity,  Geoscience Frontiers 9 (2018) 1903–1909. doi:10.1016/j.gsf.2017.11.015.
  * Foulger et al. [2018] G. R. Foulger, M. P. Wilson, J. G. Gluyas, B. R. Julian, R. J. Davies,  Global review of human-induced earthquakes,  Earth-Science Reviews 178 (2018) 438–514. doi:10.1016/j.earscirev.2017.07.008.
  * Keranen and Weingarten [2018] K. M. Keranen, M. Weingarten,  Induced Seismicity,  Annual Review of Earth and Planetary Sciences 46 (2018) 149–174. doi:10.1146/annurev-earth-082517-010054.
  * Segall et al. [1994] P. Segall, J.-R. Grasso, A. Mossop,  Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France,  Journal of Geophysical Research: Solid Earth 99 (1994) 15423–15438. doi:10.1029/94jb00989.
  * Hettema et al. [2000] M. H. H. Hettema, P. M. T. M. Schutjens, B. J. M. Verboom, H. J. Gussinklo,  Production-Induced Compaction of a Sandstone Reservoir: The Strong Influence of Stress Path,  SPE Reservoir Evaluation & Engineering 3 (2000) 342–347. doi:10.2118/65410-pa.
  * Candela et al. [2019] T. Candela, S. Osinga, J.-P. Ampuero, B. Wassing, M. Pluymaekers, P. A. Fokker, J.-D. van Wees, H. A. de Waal, A. G. Muntendam-Bos,  Depletion-Induced Seismicity at the Groningen Gas Field: Coulomb Rate-and-State Models Including Differential Compaction Effect,  Journal of Geophysical Research: Solid Earth 124 (2019) 7081–7104. doi:10.1029/2018jb016670.
  * Buijze et al. [2017] L. Buijze, P. A. van den Bogert, B. B. Wassing, B. Orlic, J. ten Veen,  Fault reactivation mechanisms and dynamic rupture modelling of depletion-induced seismic events in a Rotliegend gas reservoir,  Netherlands Journal of Geosciences 96 (2017) s131–s148. doi:10.1017/njg.2017.27.
  * Chang and Segall [2016] K. W. Chang, P. Segall,  Injection-induced seismicity on basement faults including poroelastic stressing,  Journal of Geophysical Research: Solid Earth 121 (2016) 2708–2726. doi:10.1002/2015jb012561.
  * Rutqvist et al. [2016] J. Rutqvist, A. P. Rinaldi, F. Cappa, P. Jeanne, A. Mazzoldi, L. Urpi, Y. Guglielmi, V. Vilarrasa,  Fault activation and induced seismicity in geological carbon storage - Lessons learned from recent modeling studies,  Journal of Rock Mechanics and Geotechnical Engineering 8 (2016) 789–804. doi:10.1016/j.jrmge.2016.09.001.
  * Eyre et al. [2019] T. S. Eyre, D. W. Eaton, D. I. Garagash, M. Zecevic, M. Venieri, R. Weir, D. C. Lawton,  The role of aseismic slip in hydraulic fracturing-induced seismicity,  Science Advances 5 (2019). doi:10.1126/sciadv.aav7172.
  * Tan et al. [2020] Y. Tan, J. Hu, H. Zhang, Y. Chen, J. Qian, Q. Wang, H. Zha, P. Tang, Z. Nie,  Hydraulic Fracturing Induced Seismicity in the Southern Sichuan Basin Due to Fluid Diffusion Inferred From Seismic and Injection Data Analysis,  Geophysical Research Letters 47 (2020) e2019GL084885. doi:10.1029/2019gl084885.
  * Deflandre et al. [2018] J. P. Deflandre, J. Laurent, D. Michon, E. Blondin,  Microseismic surveying and repeated VSPs for monitoring an underground gas storage reservoir using permanent geophones,  First Break 13 (2018) 129–138. doi:10.3997/1365-2397.1995008.
  * Cesca et al. [2014] S. Cesca, F. Grigoli, S. Heimann, A. Gonzalez, E. Buforn, S. Maghsoudi, E. Blanch, T. Dahm,  The 2013 September-October seismic sequence offshore Spain: a case of seismicity triggered by gas injection?,  Geophysical Journal International 198 (2014) 941–953. doi:10.1093/gji/ggu172.
  * Zhou et al. [2017] P. Zhou, H. Yang, B. Wang, J. Zhuang,  Seismological Investigations of Induced Earthquakes Near the Hutubi Underground Gas Storage Facility,  Journal of Geophysical Research: Solid Earth 124 (2017) 8753–8770. doi:10.1029/2019jb017360.
  * Hager and Toksoz [2009] B. H. Hager, M. N. Toksoz, Technical Review of Bergermeer Seismicity Study TNO Report 2008-U-R1071/B, Technical Report, Massachusetts Institute of Technology, 2009.
  * Kraaijpoel et al. [2013] D. Kraaijpoel, D. Nieuwland, B. Dost,  Microseismic Monitoring and Subseismic Fault Detection in an Underground Gas Storage,  in: Proc. 4th EAGE Passive Seismic Workshop, Amsterdam, Netherlands, European Association of Geoscientists & Engineers, 2013, pp. 1–3. doi:10.3997/2214-4609.20142354.
  * TNO [2015] TNO, Injection-Related Induced Seismicity and its Relevance to Nitrogen Injection: Description of Dutch Field Cases, Technical Report, TNO R10906, 2015.
  * Nederlandse Aardolie Maatschappij BV [2016] Nederlandse Aardolie Maatschappij BV, Norg UGS fault reactivation study and implications for seismic threat, Technical Report, Nam Ep201610208045, 2016.
  * Plaat [2009] H. Plaat,  Underground gas storage: Why and how,  Geological Society, London, Special Publications 313 (2009) 25–37. doi:10.1144/sp313.4.
  * Cornot-Gandolphe [2019] S. Cornot-Gandolphe, Underground Gas Storage in the World - 2019 Status, Technical Report, CEDIGAZ Insights # 35, 2019.
  * Zakeri and Syri [2015] B. Zakeri, S. Syri,  Electrical energy storage systems: A comparative life cycle cost analysis,  Renewable and Sustainable Energy Reviews 42 (2015) 569–596. doi:10.1016/j.rser.2014.10.011.
  * Mouli-Castillo et al. [2019] J. Mouli-Castillo, M. Wilkinson, D. Mignard, C. McDermott, R. S. Haszeldine, Z. K. Shipton,  Inter-seasonal compressed-air energy storage using saline aquifers,  Nature Energy 4 (2019) 131–139. doi:10.1038/s41560-018-0311-0.
  * Sopher et al. [2019] D. Sopher, C. Juhlin, T. Levendal, M. Erlström, K. Nilsson, J. P. Da Silva Soares,  Evaluation of the subsurface compressed air energy storage (CAES) potential on Gotland, Sweden,  Environmental Earth Sciences 78 (2019). doi:10.1007/s12665-019-8196-1.
  * Mosleh et al. [2019] M. H. Mosleh, M. Sedighi, M. Babaei, M. Turner, Geological sequestration of carbon dioxide, Academic Press, 2019, pp. 487–500. doi:10.1016/B978-0-12-814104-5.00016-8.
  * Teatini et al. [2011] P. Teatini, N. Castelletto, M. Ferronato, G. Gambolati, C. Janna, E. Cairo, D. Marzorati, D. Colombo, A. Ferretti, A. Bagliani, F. Bottazzi,  Geomechanical response to seasonal gas storage in depleted reservoirs: A case study in the Po River basin, Italy,  Journal of Geophysical Research: Earth Surface 116 (2011). doi:10.1029/2010jf001793.
  * European Commission [2023] European Commission, Questions and Answers on the new EU rules on gas storage, https://ec.europa.eu/commission/presscorner/detail/en/qanda_22_1937, 2023. Accessed: 2023-06-28.
  * Li et al. [2022] H. Z. Li, P. M. Saint-Vincent, M. Mundia-Howe, N. J. Pekney,  A national estimate of US underground natural gas storage incident emissions,  Environmental Research Letters 17 (2022) 084013. doi:10.1088/1748-9326/ac8069.
  * HiQuake [2023] HiQuake, The Human-Induced Earthquake Database, https://inducedearthquakes.org, 2023\. Accessed: 2023-06-28.
  * Gautier [2003] D. L. Gautier, Carboniferous-Rotliegend Total Petroleum System Description and Assessment Results Summary, Technical Report, U. S. Geological Survey Bulletin 2211, 2003\.
  * van Wees et al. [2014] J. van Wees, L. Buijze, K. Thienen-Visser, M. Nepveu, B. Wassing, B. Orlic, P. Fokker,  Geomechanics response and induced seismicity during gas field depletion in the Netherlands,  Geothermics 52 (2014) 206–219. doi:10.1016/j.geothermics.2014.05.004.
  * Uta [2017] P. Uta, Recent Intraplate Earthquakes in Northwest Germany - Glacial Isostatic Adjustment and/or a Consequence of Hydrocarbon Production, Ph.D. thesis, Leibniz University Hannover, 2017\. doi:10.15488/9088.
  * van Thienen-Visser and Breunese [2015] K. van Thienen-Visser, J. N. Breunese,  Induced seismicity of the Groningen gas field: History and recent developments,  The Leading Edge 34 (2015) 664–671. doi:10.1190/tle34060664.1.
  * van Wees et al. [2017] J. van Wees, P. A. Fokker, K. Van Thienen-Visser, B. B. Wassing, S. Osinga, B. Orlic, S. A. Ghouri, L. Buijze, M. Pluymaekers,  Geomechanical models for induced seismicity in the Netherlands: inferences from simplified analytical, finite element and rupture model approaches,  Netherlands Journal of Geosciences 96 (2017) s183–s202. doi:10.1017/njg.2017.38.
  * Zbinden et al. [2017] D. Zbinden, A. P. Rinaldi, L. Urpi, S. Wiemer,  On the physics-based processes behind production-induced seismicity in natural gas fields,  Journal of Geophysical Research: Solid Earth 122 (2017) 3792–3812. doi:10.1002/2017jb014003.
  * Haug et al. [2018] C. Haug, J.-A. Nüchter, A. Henk,  Assessment of geological factors potentially affecting production-induced seismicity in North German gas fields,  Geomechanics for Energy and the Environment 16 (2018) 15–31. doi:10.1016/j.gete.2018.04.002.
  * Nagelhout and Roest [1997] A. Nagelhout, J. Roest,  Investigating fault slip in a model of an underground gas storage facility,  International Journal of Rock Mechanics and Mining Sciences 34 (1997) 212.e1–212.e14. doi:10.1016/s1365-1609(97)00051-8.
  * Orlic et al. [2013] B. Orlic, B. B. T. Wassing, C. R. Geel,  Field scale geomechanical modeling for prediction of fault stability during underground gas storage operations in a depleted gas field in the Netherlands,  in: 47th US Rock Mechanics / Geomechanics Symposium, Paper #ARMA 13-300, American Rock Mechanics Association, 2013, pp. 1–11.
  * Teatini et al. [2019] P. Teatini, M. Ferronato, A. Franceschini, M. Frigo, C. Janna, C. Zoccarato,  Gas storage in compartmentalized reservoirs: a numerical investigation on possible “unexpected” fault activation,  in: 53rd US Rock Mechanics / Geomechanics Symposium, Paper #ARMA 19-1991, American Rock Mechanics Association, 2019, pp. 1–9.
  * Teatini et al. [2020] P. Teatini, C. Zoccarato, M. Ferronato, A. Franceschini, M. Frigo, C. Janna, G. Isotton,  About geomechanical safety for UGS activities in faulted reservoirs,  Proceedings of the International Association of Hydrological Sciences 382 (2020) 539–545. doi:10.5194/piahs-382-539-2020.
  * Teatini et al. [2023] P. Teatini, S. Baldan, C. Zoccarato, M. Ferronato, A. Franceschini, M. Frigo, C. Janna, G. Isotton,  Unexpected fault activation in underground gas storage. Part II: Safe operational bandwidth for reservoirs in the Roetligend formation, the Netherlands.,  Geomechanics for Energy and the Environment (2023) submitted.
  * Franceschini et al. [2016] A. Franceschini, M. Ferronato, C. Janna, P. Teatini,  A novel Lagrangian approach for the stable numerical simulation of fault and fracture mechanics,  Journal of Computational Physics 314 (2016) 503–521. doi:10.1016/j.jcp.2016.03.032.
  * Franceschini et al. [2019] A. Franceschini, N. Castelletto, M. Ferronato,  Block preconditioning for fault/fracture mechanics saddle-point problems,  Computer Methods in Applied Mechanics and Engineering 344 (2019) 376–401. doi:10.1016/j.cma.2018.09.039.
  * Isotton et al. [2019] G. Isotton, P. Teatini, M. Ferronato, C. Janna, N. Spiezia, S. Mantica, G. Volonte,  Robust numerical implementation of a 3D rate-dependent model for reservoir geomechanical simulations,  International Journal for Numerical and Analytical Methods in Geomechanics 43 (2019) 2752–2771. doi:10.1002/nag.3000.
  * Gambolati et al. [2000] G. Gambolati, P. Teatini, D. Baú, M. Ferronato,  Importance of poroelastic coupling in dynamically active aquifers of the Po river basin, Italy,  Water Resources Research 36 (2000) 2443–2459. doi:10.1029/2000wr900127.
  * Buijze et al. [2019] L. Buijze, P. A. J. van den Bogert, B. B. T. Wassing, B. Orlic,  Nucleation and Arrest of Dynamic Rupture Induced by Reservoir Depletion,  Journal of Geophysical Research: Solid Earth 124 (2019) 3620–3645. doi:10.1029/2018JB016941.
  * Bizzarri [2010] A. Bizzarri,  How to promote earthquake ruptures: different nucleation strategies in a dynamic model with slip-weakening function,  Bulletin of the Seismological Society of America 100 (2010) 923–940. doi:10.1785/0120090179.
  * Kikuchi and Oden [1988] N. Kikuchi, J. T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1988. doi:10.1137/1.9781611970845.
  * Laursen [2003] T. A. Laursen, Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis, Springer-Verlag Berlin Heidelberg, 2003. doi:10.1007/978-3-662-04864-1.
  * Wriggers [2006] P. Wriggers, Computational Contact Mechanics, 2nd ed., Springer-Verlag Berlin Heidelberg, 2006. doi:10.1007/978-3-540-32609-0.
  * Wohlmuth [2011] B. Wohlmuth,  Variationally consistent discretization schemes and numerical algorithms for contact problems,  Acta Numerica 20 (2011) 569–734. doi:10.1017/s0962492911000079.
  * Simo and Hughes [1998] J. C. Simo, T. J. R. Hughes, Computational Inelasticity, Springer-Verlag New York, 1998. doi:10.1007/b98904.
  * Franceschini et al. [2020] A. Franceschini, N. Castelletto, J. A. White, H. A. Tchelepi,  Algebraically stabilized Lagrange multiplier method for frictional contact mechanics with hydraulically active fractures,  Computer Methods in Applied Mechanics and Engineering 368 (2020) 113161\. doi:10.1016/j.cma.2020.113161.
  * Franceschini et al. [2022] A. Franceschini, N. Castelletto, J. A. White, H. A. Tchelepi,  Scalable preconditioning for the stabilized contact mechanics problem,  Journal of Computational Physics 459 (2022) 111150. doi:10.1016/j.jcp.2022.111150.
  * Coussy [2004] O. Coussy, Poromechanics, John Wiley & Sons Ltd. Chichester, England, 2004\. doi:10.1002/0470092718.
  * De Souza Neto et al. [2011] E. A. De Souza Neto, D. Peric, D. R. Owen, Computational methods for plasticity: theory and applications, John Wiley & Sons, 2011. doi:10.1002/9780470694626.
  * Xikui and Zienkiewicz [1992] L. Xikui, O. Zienkiewicz,  Multiphase flow in deforming porous media and finite element solutions,  Computers & structures 45 (1992) 211–227. doi:10.1016/0045-7949(92)90405-O.
  * Rutqvist et al. [2002] J. Rutqvist, Y.-S. Wu, C.-F. 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 (2002) 429–442. doi:10.1016/S1365-1609(02)00022-9.
  * Yadigaroglu and Hewitt [2018] G. Yadigaroglu, G. F. Hewitt, Introduction to Multiphase Flow. Basic Concepts, Applications and Modelling, Springer, 2018\. doi:10.1007/978-3-319-58718-9.
  * Wilkins et al. [2021] A. Wilkins, C. P. Green, J. Ennis-King,  An open-source multiphysics simulation code for coupled problems in porous media,  Computers & Geosciences 154 (2021) 104820. doi:10.1016/j.cageo.2021.104820.
  * Garipov et al. [2016] T. T. Garipov, M. Karimi-Fard, H. A. Tchelepi,  Discrete fracture model for coupled flow and geomechanics,  Computational Geosciences 20 (2016) 149–160. doi:10.1007/s10596-015-9554-z.
  * Zoccarato et al. [2016] C. Zoccarato, D. Baù, M. Ferronato, G. Gambolati, A. Alzraiee, P. Teatini,  Data assimilation of surface displacements to improve geomechanical parameters of gas storage reservoirs,  Journal of Geophysical Research: Solid Earth 121 (2016) 1441–1461. doi:10.1002/2015jb012090.
  * Teatini et al. [2014] P. Teatini, N. Castelletto, G. Gambolati,  3D geomechanical modeling for CO2 geological storage in faulted formations. A case study in an offshore northern Adriatic reservoir, Italy,  International Journal of Greenhouse Gas Control 22 (2014) 63–76. doi:10.1016/j.ijggc.2013.12.021.
  * Castelletto et al. [2013a] N. Castelletto, G. Gambolati, P. Teatini,  Geological CO2 sequestration in multi-compartment reservoirs: Geomechanical challenges,  Journal of Geophysical Research: Solid Earth 118 (2013a) 2417–2428. doi:10.1002/jgrb.50180.
  * Castelletto et al. [2013b] N. Castelletto, M. Ferronato, G. Gambolati, C. Janna, P. Teatini,  Compartmentalization Effects in Geologic CO2 Sequestration. A case Study in an Offshore Reservoir in Italy,  in: Poromechanics V: Proceedings of the Fifth Biot Conference on Poromechanics, 2013b, pp. 1265–1270. doi:10.1061/9780784412992.151.
  * Castelletto et al. [2013c] N. Castelletto, P. Teatini, G. Gambolati, D. Bossie-Codreanu, O. Vincké, J.-M. Daniel, A. Battistelli, M. Marcolini, F. Donda, V. Volpi,  Multiphysics modeling of CO2 sequestration in a faulted saline formation in Italy,  Advances in water resources 62 (2013c) 570–587. doi:10.1016/j.advwatres.2013.04.006.
  * Eymard et al. [2000] R. Eymard, T. Gallouët, R. Herbin,  Finite volume methods,  Handbook of numerical analysis 7 (2000) 713–1018. doi:10.1016/S1570-8659(00)07005-8.
  * Manzini and Russo [2008] G. Manzini, A. Russo,  A finite volume method for advection–diffusion problems in convection-dominated regimes,  Computer Methods in Applied Mechanics and Engineering 197 (2008) 1242–1261. doi:10.1016/j.cma.2007.11.014.
  * Droniou [2014] J. Droniou,  Finite volume schemes for diffusion equations: introduction to and review of modern methods,  Mathematical Models and Methods in Applied Sciences 24 (2014) 1575–1619. doi:10.1142/S0218202514400041.
  * Wan et al. [2003] J. Wan, L. Durlofsky, T. Hughes, K. Aziz,  Stabilized finite element methods for coupled geomechanics-reservoir flow simulations,  in: SPE Reservoir Simulation Symposium, 2003, pp. 1–11. doi:10.2118/79694-MS.
  * Jha and Juanes [2007] B. Jha, R. Juanes,  A locally conservative finite element framework for the simulation of coupled flow and reservoir geomechanics,  Acta Geotechnica 2 (2007) 139–153. doi:10.1007/s11440-007-0033-0.
  * Prevost and Sukumar [2016] J. H. Prevost, N. Sukumar,  Faults simulations for three-dimensional reservoir-geomechanical models with the extended finite element method,  Journal of the Mechanics and Physics of Solids 86 (2016) 1–18. doi:10.1016/j.jmps.2015.09.014.
  * Abushaikha et al. [2017] A. S. Abushaikha, D. V. Voskov, H. A. Tchelepi,  Fully implicit mixed-hybrid finite-element discretization for general purpose subsurface reservoir simulation,  J. Comput. Phys. 346 (2017) 514–538. doi:10.1016/j.jcp.2017.06.034.
  * Nardean et al. [2023] S. Nardean, M. Ferronato, A. Abushaikha,  Block constrained pressure residual preconditioning for two-phase flow in porous media by mixed hybrid finite elements,  Computational Geosciences (2023). doi:10.1007/s10596-023-10238-x.
  * The Open Porous Media Initiative [2023] The Open Porous Media Initiative, Open Porous Media, https://opm-project.org, 2023\. Accessed: 2023-07-21.
  * Rasmussen et al. [2021] A. F. Rasmussen, T. H. Sandve, K. Bao, A. Lauser, J. Hove, B. Skaflestad, R. Klöfkorn, M. Blatt, A. B. Rustad, O. Sævareid, K.-A. Lie, A. Thune,  The open porous media flow reservoir simulator,  Computers & Mathematics with Applications 81 (2021) 159–185. doi:10.1016/j.camwa.2020.05.014.
  * Ferronato et al. [2012] M. Ferronato, C. Janna, G. Pini,  Parallel solution to ill-conditioned FE geomechanical problems,  International Journal for Numerical and Analytical Methods in Geomechanics 36 (2012) 422–437. doi:10.1002/nag.1012.
  * Nocedal and Wright [2006] J. Nocedal, S. Wright, Numerical optimization, Springer Science & Business Media, 2006. doi:10.1007/b98874.
  * Antil et al. [2018] H. Antil, D. P. Kouri, M.-D. Lacasse, D. Ridzal, Frontiers in PDE-constrained Optimization, volume 163, Springer, 2018. doi:10.1007/978-1-4939-8636-1.
  * Goodman et al. [1968] R. E. Goodman, R. L. Taylor, T. L. Brekke,  A model for the mechanics of jointed rock,  Journal of Soil Mechanics & Foundations Div. 94 (1968) 637–659.
  * Elman et al. [2014] H. C. Elman, D. J. Silvester, A. J. Wathen, Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics, Oxford University Press, USA, 2014. doi:10.1093/acprof:oso/9780199678792.001.0001.
  * The Atlas Project Team [2023] The Atlas Project Team, Atlas web page, https://www.m3eweb.it/atlas, 2023\. URL: https://www.m3eweb.it/atlas.
  * Benzi et al. [2005] M. Benzi, G. H. Golub, J. Liesen,  Numerical solution of saddle point problems,  Acta Numerica 14 (2005) 1–137. doi:10.1017/s0962492904000212.
  * Franceschini et al. [2022] A. Franceschini, M. Ferronato, M. Frigo, C. Janna,  A reverse augmented constraint preconditioner for Lagrange multiplier methods in contact mechanics,  Computer Methods in Applied Mechanics and Engineering 392 (2022) 114632\. doi:10.1016/j.cma.2022.114632.
  * Andrews [1976] D. J. Andrews,  Rupture velocity of plane strain shear cracks,  Journal of Geophysical Research 81 (1976) 5679–5687. doi:10.1029/JB081i032p05679.
  * Ida [1972] Y. Ida,  Cohesive force across the tip of a longitudinal-shear crack and Griffith’s specific surface energy,  Journal of Geophysical Research 77 (1972) 3796–3805. doi:10.1029/JB077i020p03796.
  * Palmer and Rice [1973] A. C. Palmer, J. R. Rice,  The growth of slip surfaces in the progressive failure of over-consolidated clay,  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 332 (1973) 527–548. doi:10.1098/rspa.1973.0040.
  * Hunfeld et al. [2020] L. B. Hunfeld, J. Chen, S. Hol, A. R. Niemeijer, C. J. Spiers,  Healing Behavior of Simulated Fault Gouges From the Groningen Gas Field and Implications for Induced Fault Reactivation,  Journal of Geophysical Research: Solid Earth 125 (2020) e2019JB018790. doi:10.1029/2019JB018790.
  * Fokker et al. [2013] P. Fokker, B. Wassing, F. Van Leijen, R. Hanssen, D. Nieuwland,  Data assimilation of PS-InSAR movement measurements applied to the Bergermeer gas field,  in: International EAGE Workshop on Geomechanics and Energy, 2013, pp. cp–369. doi:10.3997/2214-4609.20131971.
  * Wassing et al. [2017] B. Wassing, L. Buijze, J. Ter Heege, B. Orlic, S. Osinga,  The impact of viscoelastic caprock on fault reactivation and fault rupture in producing gas fields,  in: 51st US Rock Mechanics/Geomechanics Symposium, 2017, pp. 1–13.
  * Plaat [2009] H. Plaat,  Underground gas storage: Why and how,  Geological Society, London, Special Publications 313 (2009) 25–37. doi:10.1144/SP313.4.