Proceedings Volume Cover
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  
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  
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 60and 0 m offset.  
based on  
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.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.  
Pgrad,vert (initial)  
Pgrad,hor (initial)  
Pgrad,water (initial)  
(GPa) (-)  
Sandstone reservoir 2200  
Underburden 2200  
0.17 0.022  
0.17 0.022  
0.17 0.022  
0.17 0.022  
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  
σ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  
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:  
푚푎푥 = 퐶 + 휇푠푡푎′  
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.  
푆퐶푈 = 휏  
= (휇  
where τ is the shear stress on the fault. Figure 2b shows  
that fault criticality is highest at the top of the reservoir  
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−푣  
= 푑휎= 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  
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.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  
퐶 휇 퐷  
0 = (휏 −휏 )  
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,