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 dP_{nucl, }and L_{nucl }
respectively.
∗
퐶 휇 퐷
1
푐
퐿_{0 }= _{(휏 −휏 ) }
(5)
푠
푑
where L_{0 }is the size of the breakdown zone, C_{1 }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 L_{0 }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 dP_{nucl }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,