ARMA 16-658
Modelling of fault reactivation and fault slip in producing
gas fields using a slip-weakening friction law
Wassing, B.B.T., Buijze, L. and Orlic, B.
TNO, Earth, Environmental and Life Sciences, Utrecht, The Netherlands
This paper was prepared for presentation at the 50th US Rock Mechanics / Geomechanics Symposium held in Houston, Texas, USA, 26-29 June
2016. This paper was selected for presentation at the symposium by an ARMA Technical Program Committee based on a technical and critical
review of the paper by a minimum of two technical reviewers. The material, as presented, does not necessarily reflect any position of ARMA, its
officers, or members.
ABSTRACT: Geomechanical numerical simulations were conducted to analyze the stability of faults during gas production. A
FLAC3D model of a fault intersecting a producing gas reservoir was developed which incorporates the fully dynamic behavior of
the fault and surrounding rock mass, and a fault frictional behavior based on a slip-weakening law. The simplified reservoir and
fault geometry of the model are representative for the Dutch gas fields. Simulated fault slip displacements and fault slip lengths
were used to calculate moment magnitudes of induced seismic events. In addition, results of the fully dynamic model with slip-
weakening frictional behavior were compared to results of a static geomechanical model with slip-weakening. Comparison of
model results shows that, for a first-order assessment of induced seismicity, the static models can be used as a simplified and
computationally less expensive alternative to the fully dynamic fault rupture models.
between in-situ stress conditions, pore pressure
1. INTRODUCTION
depletion, (differential) reservoir compaction, reservoir
and fault stress evolution, and the onset of fault
reactivation. In many cases poro-elastic stress changes
and an ideal-plastic failure law for the faults are used to
explain the temporal and spatial distribution of induced
seismicity caused by reservoir depletion (e.g. Mulders,
2003, Orlic et al., 2013, Van Den Bogert, 2015). Such
models give valuable insights into the potential for fault
reactivation and are generally suitable to identify faults
which are more or less prone to reactivation and to
assess the onset of fault reactivation in a non-critical
stress regime. However, these models give little
information on the post-failure dynamic behavior once
the fault slip has been initiated. Insight into the post-
failure behavior of faults is of paramount importance to
understand the dynamic rupture behavior of the faults
and associated seismicity, and to estimate source
parameters such as the area of fault slip, stress drops and
(a)seismic slip displacements with ongoing production.
These parameters can be used to assess (trends in)
magnitudes of induced earthquakes due to ongoing gas
production, and are considered an important input to
seismic hazard analysis.
The operation of a gas field causes dynamic changes of
the pore pressure and therefore changes in the stress state
of the reservoir and surrounding rocks. Production-
induced stress changes can destabilize faults which
transect or bound the reservoir, or are located in the
vicinity of the reservoir, causing associated seismicity.
Whether or not faults are (seismically) reactivated during
reservoir depletion depends on a complex interplay of
many factors, such as the combination of initial stress
state of the reservoir and faults, the magnitude of pore
pressure changes, the geometry of the reservoir and
faults, such as orientation, dip and fault offset, and the
geomechanical properties of the rocks and faults.
In many onshore producing gas fields in the Netherlands
seismicity has been observed during the depletion of the
gas reservoirs. Induced seismicity has been recorded in
26 out of 150 producing gas fields, with maximum
magnitudes of the seismic events up to Mw 3.6. In all of
the gas fields pressure depletion was significant prior to
the onset of recorded seismicity (Van Wees et al., 2014).
The marked delay between the start of reservoir
depletion and the onset of seismicity has been
interpreted as an indication that the in-situ stress
conditions in the northern part of the Netherlands are
non-critical.
The main objective of the study presented in this paper is
to obtain a better understanding of the process of fault
reactivation, rupture, related stress drops and (a)seismic
slip on faults within a gas reservoir which compacts due
to the ongoing extraction of gas and pore pressure
depletion. The study focusses on the evolution of
Various publications describe the use of analytical and
numerical geomechanical models to explain the relation
stresses, displacements and reactivated fault area after
the onset of fault reactivation, taking into account the
post-failure behavior of the fault. A geomechanical finite
difference model of a fault intersecting a producing gas
reservoir was built in FLAC3D (Itasca, 2013). The
model was used to analyze fault reactivation and rupture
caused by the decline of pore pressures. The model is
overlain by a 50 m thick caprock of Zechstein rocksalt.
Above the Zechstein caprock and below the reservoir,
rocks consist of undifferentiated overburden,
respectively underburden. Reservoir, overburden and
underburden are cut by a continuous fault plane with a
dip of 60⁰ and 0 m offset.
based on
a
simplified, schematic cross-section
The size of the FLAC3D model is selected in such a way
as to minimize the effects of boundaries on the stress and
deformation. The model covers a domain of 5000 m *5
m*6000 m. The vertical boundaries are fixed laterally to
constrain displacements perpendicular to the sides of the
model. The lower horizontal boundary of the model is
fixed vertically. The resolution of the finite difference
mesh increases towards the fault surface and the
boundaries of the reservoir. A high model resolution of 1
m in the vertical direction and 5 m in the horizontal
direction is used around the fault segments intersecting
the reservoir (for a discussion of choice of vertical
resolution, see also section 3). The element sizes
increase up to a maximum size of a few hundred meters
towards the outer boundaries of the model. At the
location of the fault surface an interface is defined. As
only one, 5 m thick element is defined in the out-of-
plane y-direction of the model, deformation conditions
in the model are plane strain (i.e. no displacements are
allowed in the y-direction, see also Figure 1).
representing a reservoir and fault geometry characteristic
for a typical Rotliegend gas field in the northern part of
the Netherlands. A slip-weakening law for frictional
failure has been implemented in FLAC3D to describe
the post-failure behavior of the fault.
In this study two types of numerical analyses were
performed: a static and a fully dynamic analysis. Results
of the fully dynamic models, in terms of area of fault
reactivation within and outside the reservoir, distribution
of fault slip, stress drop and related maximum
magnitude, were compared to results obtained for the
static model runs. This way the contribution of mass and
inertia forces to the stress drop, total fault slip, slip area
and magnitude forecasts was analyzed. The rationale
behind this comparison is that, if effects of inertia forces
are limited, static models can be applied as a simplified
and computationally less expensive alternative to the
fully dynamic calculations.
In the next section 2 a brief description of the geometry,
rock parameters, initial stresses and loading conditions
for the reference geomechanical model is given. Section
2 also describes the implementation of the constitutive
law for post-failure fault behavior and the modelling
approach for static and dynamic fault rupture analysis. In
section 3 the static and dynamic modelling results for the
reference scenario are summarized and compared.
Subsequently, in section 4, the effect of stress drop on
fault slip, reactivated fault area and difference between
static and dynamic models is analyzed. In the last section
5 the implications of this work for modelling of fault
reactivation are addressed and main conclusions from
this work are summarized.
2. MODEL DESCRIPTION AND APPROACH
2.1. Model geometry, rock parameters, initial
stress and pore pressure conditions
Model geometry and rock parameters of the reference
scenario are chosen such as to capture the main aspects
of a typical depleting Rotliegend reservoir in the
northern parts of the Netherlands. Figure 1 gives a
schematic representation of the geometry of the model
and depths of individual layers. The model comprises a
200 m thick Rotliegend sandstone reservoir at a depth of
2800 - 3000 m below surface level (bsl). The reservoir is
Figure 1. Schematic representation of the FLAC3D model
geometry. Fault dip in the reference scenario is 60⁰.
Table 1. Geomechanical properties, initial pressure and stress gradients for the base scenario.
Lithology
E
v
Pgrad,vert (initial)
MPa/m
Pgrad,hor (initial)
MPa/m
Pgrad,water (initial)
MPa/m
Overpressure
(MPa)
ρ
(GPa) (-)
(kg/m3)
2200
2200
Overburden
Rocksalt
Sandstone reservoir 2200
Underburden 2200
15
15
15
15
0.17 0.022
0.17 0.022
0.17 0.022
0.17 0.022
0.016
0.022
0.016
0.016
0.011
0.011
0.002
0.011
-
-
4.8
-
The overburden layers, Zechstein rocksalt caprock,
reservoir rocks and the underburden layers are modelled
as elastic materials. Table 1 summarizes rock densities
and elastic rock parameters, i.e. Young’s modulus (E)
and Poisson’s ratio (v) of the reservoir rock, Zechstein
caprock and burden. The elastic properties of the
reservoir rocks are based on experimental tests on
Rotliegend sandstones cored in the Dutch gas fields. For
this schematic model representative of a typical
Rotliegend field, average rock properties are used for the
reservoir rocks. Though large variations in lithology of
burden and reservoir and related rock properties can
exist, an extensive sensitivity analysis of parameter
uncertainty is not within the scope of this study.
Whereas stress paths in a depleting reservoir depend on
the Poisson’s ratio of the reservoir rocks (see also
section 2.2), in a horizontally layered model geometry
like the base scenario, with mainly uniaxial reservoir
compaction, the effect of elasticity parameters of the
caprock and burden on stress and deformation is
expected to be limited. In order to facilitate the
interpretation of the dynamic rupture models,
differentiation of over- and underburden is kept to a
minimum. Densities and elasticity parameters of all
layers are kept uniform and chosen equal to reservoir
parameters.
σv=σhmin=σHmax). The high initial horizontal stresses
impede fault rupture through the visco-elastic Zechstein
caprock layers. Figure 2a presents the pore pressure,
horizontal and vertical effective stress in the rocks and
the resulting normal effective stress and shear stress on
the fault. Stress path gradients are summarized in Table
1.
The reservoir rocks are assumed to be filled with gas,
with pore pressure gradients in the reservoir of 0.002
MPa/m. Pore pressure gradients in the Zechstein
caprock, overburden and underburden are hydrostatic,
i.e. 0.011 MPa/m. The reservoir is assumed to be
initially overpressured, i.e. with a pore pressure of 35
MPa at the top of the reservoir at 2800 m bsl, since
many Rotliegend reservoirs in the northern Netherlands
are overpressured (Verweij et al., 2014). The average
vertical lithostatic stress gradient for σv in reservoir,
rocksalt and burden is 0.022 MPa/m, consistent with a
rock density of 2200 kg/m3 for all layers. Stress
gradients for total minimum (σhmin) and maximum
horizontal stress (σHmax) conditions in the reservoir rocks
and burden are 0.016 MPa/m, which is regarded as a
representative average value for the northern
Netherlands (Van Wees et al., 2014). Stress conditions
in the Zechstein rocksalt are assumed to be isotropic due
to visco-elastic behavior of the rocksalt (K0=1 and
Figure 2. a) Vertical effective (σ’v) and horizontal effective
(σ’h) stress, and pore pressure (P) in the reservoir, caprock and
burden and resulting normal effective stress (σ’n), shear stress
(τ) and b) shear capacity utilization (SCU) on the fault. Beige-
shaded area presents sandstone reservoir, blue-shaded area
represents Zechstein rocksalt. Black circle indicates region of
highest initial fault criticality at the reservoir top.
The initial strength of the fault is determined by the
Mohr Coulomb constitutive model:
휏푚푎푥 = 퐶 + 휇푠푡푎푡휎푛′
(1)
where τmax is the shear strength of the fault or the
maximum shear stress the fault can resist before it fails,
C is cohesion, µstat is friction coefficient (i.e. static
friction coefficient before the onset of failure) and σ’n is
effective normal stress on the fault. The fault is assumed
to be cohesionless. The friction coefficient is μstat = 0.6,
which is in agreement with typical values for sandstone.
The graph in Figure 2b presents the Shear Capacity
Utilization (SCU) at the fault for initial stress conditions.
SCU gives the proximity to failure and is defined as
(Mulders, 2003):
friction during fault slip. Various publications have
addressed the use of velocity-dependent rate-and-state
friction laws (McClure and Horne, 2011, Dieterich et al.,
2015) and slip-weakening friction laws for modelling
fault friction and induced seismic fault response (Cappa
and Rutqvist, 2011, Cappa and Rutqvist, 2012, Buijze et
al., 2015). In this study we chose a slip-weakening
constitutive law to model the evolution of friction during
fault slip. A schematic representation of the slip
weakening law is given in Figure 3 (see also e.g.
Ohnaka, 2013). The slip weakening law is described by
the parameters of critical slip distance Dc and static and
dynamic friction coefficient. When static friction μstat is
exceeded, friction decays linearly with slip displacement
to dynamic friction μdyn. Dc is the critical displacement
for friction to degrade to µdyn. The dynamic friction
coefficient μdyn and the critical slip distance Dc in
combination with the initial stress conditions have a
large influence on the stress drop and the propagation of
slip along the fault. The dynamic friction coefficient µdyn
in the slip-weakening model is 0.5 and critical distance
Dc=0.01 m (Buijze et al., 2015). The slip-weakening law
does not take into account the effects of velocity-
dependency and healing.
휏
휏
푆퐶푈 = 휏
= (휇
)
(2)
′
휎
푚푎푥
푠푡푎푡
푛
where τ is the shear stress on the fault. Figure 2b shows
that fault criticality is highest at the top of the reservoir
rocks.
2.2. Pore pressure depletion
Pressure depletion of the reservoir rocks is equal on both
sides of the fault. It is assumed that faults are permeable
and the pore pressures in the fault are similar to pore
pressures in the reservoir rocks. Then, for a uniformly
depleting horizontally layered reservoir with a large
lateral extent, like the base scenario, the stress path
coefficients and stress path during pore pressure
depletion up to the onset of fault slip can be derived
from poro-elastic theory (Mulders, 2003):
1−2푣
훾ℎ = 푑휎 = 훼 ( 1−푣
)
(3)
(4)
ℎ
푑푃
훾푣 = 푑휎푣 = 0
푑푃
where σv is vertical total stress, σh is horizontal total
stress, dP is pore pressure change, ν is Poisson’s ratio, α
is Biot’s coefficient, γh and γv are the horizontal and
vertical stress path coefficients, which give the changes
in total horizontal and total vertical stress per unit of
pressure change.
Figure 3. a) Evolution of shear stress during fault rupture.
b) Simplified constitutive law for slip-dependent fault
frictional weakening as used in FLAC3D (after: Ohnaka,
2013). Parameter µstat denoting static friction coefficient, µdyn
dynamic friction coefficient and Dc the critical displacement
for friction to degrade to µdyn. c) Evolution of shear stress and
stress drop related to fault slip, with τ0 denoting initial shear
stress on the verge of fault slip, τd dynamic or residual shear
stress and τpeak peak shear stress.
Based on the above equations, for the base scenario with
a fault dip of 60⁰, μstat = 0.6 and initial stress and pore
pressure conditions in the reservoir as described in
section 2.1, a total pore pressure depletion of 13.4 MPa
is computed for the onset of failure at the fault segment
intersecting the top of the reservoir. In the FLAC3D
model, pore pressures are lowered by 13.4 MPa, up to
the onset of fault slip at the top of the reservoir. From
In contrast to an ideal-plastic Mohr Coulomb failure law,
where friction coefficient is constant during fault slip,
advanced constitutive laws for fault frictional behavior
can be used to model the post-failure evolution of
13.4 MPa depletion onwards, the reservoir sandstone is
further depleted with small load steps of 0.1 MPa.
and mass of the grid points is derived from the real
densities of the surrounding elements. Results from the
static calculations can be compared to results from the
fully dynamic rupture calculations, which take into
account mass and inertia forces. In this way the
contribution of inertia forces to the evolution of fault
slip, reactivated fault area and stress drop can be
quantified.
The time step in dynamic calculation is automatically
controlled by FLAC3D. The time step depends on the
element size and compressional wave speed. In the base
scenario, a constant time step of 10-6 s is used. In the
dynamic calculations, the no-displacement boundary
conditions are replaced by the so-called quiet, viscous
boundary conditions, to minimize reflections of the
outward propagating waves back into the model.
Figure 4. Stress path (in blue) in a reservoir during depletion
and release of shear stresses during fault rupture. The smaller
Mohr circle represents the initial stress state before depletion,
the larger Mohr circle represents the stress state at the onset of
fault slip.
3. RESULTS FOR THE BASE SCENARIO
3.1. Sensitivity to mesh resolution
In the static and dynamic models with slip weakening,
shear stresses and slip rates are expected to vary
significantly in the breakdown zone behind the rupture
front (Figure 3). The along-fault mesh resolution was
adjusted to the size of the breakdown zone. Here the
recommendations of Day et al. (2005) were followed to
resolve the breakdown zone with at least five grid points.
An upper limit for the size of the breakdown zone, just
after initiation of fault rupture, can be derived from the
equations given in Day et al. (2005):
The evolution of fault slip length, relative shear
displacements, fault shear stress and fault friction
coefficient is closely monitored during simulations. The
FLAC3D model is run in static mode during the first
stage of slip, in which the growth of the slip area and the
increase in slip displacements is still controlled by the
applied loading. At a certain pressure drop, due to the
increase of fault slip length and associated decrease of
apparent fault stiffness, the fault weakening rate exceeds
the rate of unloading, leading to instability. At this stage,
an acceleration in the growth of the fault slip length and
rapid increase in the slip displacements is observed, with
associated rapid drop in friction coefficient and shear
stress on the fault, which is indicative of the nucleation
of a ‘seismic event’. The growth in fault slip length and
fault slip displacements is no longer controlled by the
loading rate and fault slip extends spontaneously by the
release of elastic strain energy stored in the rocks
surrounding the fault (Ohnaka, 2013). A ‘seismic event’
in FLAC3D is defined as a sharp increase in the amount
of fault slip displacement and a simultaneous sharp
decrease in the fault shear stress and fault friction
coefficient. The pressure drop and slip length at which
nucleation occurs are denoted dPnucl, and Lnucl
respectively.
∗
퐶 휇 퐷
1
푐
퐿0 = (휏 −휏 )
(5)
푠
푑
where L0 is the size of the breakdown zone, C1 is a
constant, and μ*=G/(1-v), with G the shear modulus of
the rocks, τs is shear stress on the verge of fault slip and
τd is dynamic or residual shear stress. Based on
parameters for the reference scenario, an L0 of
approximately 35 m is computed, which gives a
minimum along-fault resolution of 7 m for the zones and
interface elements within, above and below the reservoir
layer. As the size of the breakdown zone is expected to
reduce with the increase in slip rates during the rupture
process (Day et al., 2005), the sensitivity of model
outcomes to varying along-fault mesh resolutions was
tested. For the reference scenario as defined in section
2.1 three model runs were performed with vertical mesh
resolutions between a depth of 2750 m and 3100 m bsl
of 1 m, 2 m and 5 m. Aspect ratios of the zones near the
fault do not exceed 1:5. For all 3 models, a reservoir
pore pressure depletion of 13.4 MPa was imposed,
which initiated fault slip on the top of the reservoir.
Additional pore pressure depletion with load steps of 0.1
MPa was applied, up to nucleation of a seismic event.
2.3. Static versus dynamic analysis
At dPnucl a static and dynamic model calculation is
performed. In both static and dynamic FLAC3D
calculations, the explicit finite difference approach is
used to solve the equations of motion. In a static
analysis, mass and inertia are fictitious terms, as they are
used as a means to reach the static equilibrium in a
numerically stable way. In a fully dynamic rupture
analysis, the configuration is switched to dynamic mode,