1.
Introduction
Earthquakes induced by fluid injection into the subsurface pose a major challenge for the geoenergy in
-
dustry and society in general (Ellsworth,
2013
; Grigoli et al.,
2017
). Tectonically quiescent regions, where
dormant faults could be reactivated are particularly challenging, as their infrastructure is often not designed
for large-magnitude induced earthquakes (McGarr et al.,
2015
). At the same time, some faults have been
observed to slip stably at aseismic speeds of 10
−7
–10
−2
m/s, in response to fluid injection (Cornet et al.,
1997
;
Duboeuf et al.,
2017
; Guglielmi et al.,
2015
; Scotti & Cornet,
1994
; Wei et al.,
2015
). While induced earth
-
quakes have been located anywhere from a few meters to tens of kilometers from injection wells (Goebel
& Brodsky,
2018
), the spatial extent of fluid-induced aseismic slip is not as well characterized due to the
paucity of direct observations. Understanding what conditions lead to seismic versus aseismic and localized
versus widespread fault reactivation is central to physics-based hazard forecasting.
An outstanding opportunity to investigate these questions is offered by a decametric-scale fluid injection
experiment recently conducted in an underground tunnel intercepting a dormant fault in a carbonate
Abstract
While the notion that injecting fluids into the subsurface can reactivate faults by reducing
frictional resistance is well established, the ensuing evolution of the slip is still poorly understood. What
controls whether the induced slip remains stable and confined to the fluid-affected zone or accelerates
into a runaway earthquake? Are there observable indicators of the propensity to earthquakes before
they happen? Here, we investigate these questions by modeling a unique fluid-injection experiment
on a natural fault with laboratory-derived friction laws. We show that a range of fault models with
diverging stability with sustained injection reproduce the slip measured during pressurization. Upon
depressurization, however, the most unstable scenario departs from the observations, suggesting that the
fault is relatively stable. The models could be further distinguished with optimized depressurization tests
or spatially distributed monitoring. Our findings indicate that avoiding injection near low-residual-friction
faults and depressurizing during slip acceleration could help prevent large-scale earthquakes.
Plain Language Summary
Fluid injections into the Earth's crust are common practice in
the exploitation of subsurface energy resources such as geothermal energy, shale gas, and conventional
hydrocarbons. These injections can perturb nearby fault structures and hence induce earthquakes and
transient slow slip. Understanding what controls the stability (i.e., the propensity to generate earthquakes)
and spatial extent of the fault response as well as identifying precarious faults is crucial to minimize
the seismic hazard associated with these industrial practices. Here, we take a step toward this goal by
modeling a unique experiment, in which water was injected into a natural fault and the resulting slip
measured directly at depth. We first show that multiple models can explain the observations equally well,
while pressure is increased in the experiment. In these models, how stable the fault response is with
further injection and how large of a zone is reactivated compared to the fluid-affected region depends
on frictional properties. We then demonstrate that the slow slip response to a decrease in injection
pressure further constrains the range of admissible models. Our work suggests that it may be possible to
identify potentially hazardous faults with optimally designed injection tests without inducing damaging
earthquakes.
LAROCHELLE ET AL.
© 2021. American Geophysical Union.
All Rights Reserved.
Constraining Fault Friction and Stability With Fluid-
Injection Field Experiments
Stacy Larochelle
1
, Nadia Lapusta
1,2
, Jean-Paul Ampuero
3
, and Frédéric Cappa
3,4
1
Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USA,
2
Division of
Engineering and Applied Science, California Institute of Technology, Pasadena, CA, USA,
3
IRD, CNRS, Observatoire de
la Côte d'Azur, Université Côte d'Azur, Géoazur, Sophia Antipolis, Nice, France,
4
Institut Universitaire de France, Paris,
France
Key Points:
•
Multiple frictional models with
different stability reproduce the slip
observed during the pressurization
stage of a field experiment
•
The depressurization phase
provides additional constraints on
hydromechanical parameters and
hence fault stability
•
Fault stability and the spatial extent
of slip relative to the pressurized
region depend on residual friction
versus initial stress levels
Supporting Information:
Supporting Information may be found
in the online version of this article.
Correspondence to:
S. Larochelle,
stacy.larochelle@caltech.edu
Citation:
Larochelle, S., Lapusta, N., Ampuero,
J.-P., & Cappa, F. (2021). Constraining
fault friction and stability with fluid-
injection field experiments.
Geophysical
Research Letters
,
48
, e2020GL091188.
https://doi.org/10.1029/2020GL091188
Received 7 OCT 2020
Accepted 5 MAY 2021
10.1029/2020GL091188
RESEARCH LETTER
1 of 11
Geophysical Research Letters
formation (Guglielmi et al.,
2015
) (Figure
1a
). During the experiment, the fluid pressure and fault slip were
recorded at the injection site. Although the observed slip was mostly aseismic, it is important to understand
if the observations contained sufficient information to determine whether slip would have accelerated into
an earthquake rupture if injection had continued. Previous efforts to model the field experiment with a
slip-weakening friction law concluded that aseismic slip outgrew the pressurized zone, potentially leading
to a runaway earthquake with continued injection (Bhattacharya & Viesca,
2019
).
Here, we use the data from the field experiment to examine the issue of slow and confined versus fast
and runaway slip in models with more realistic, laboratory-derived rate-and-state friction laws (Dieter
-
ich,
1979
,
2007
; Ruina,
1983
) consistent with laboratory results on materials from this specific fault zone
(Cappa et al.,
2019
). Furthermore, we use the modeling to identify promising avenues to quantify the fault
properties and control injection-induced seismicity hazard. We adopt a fully dynamic computational frame
-
work that resolves both aseismic and seismic slip on faults. We keep other model ingredients relatively
simple to better understand frictional effects in the presence of a diffusing fluid. For example, we do not
explicitly model the change in fault permeability induced by slip as in previous studies (Bhattacharya &
Viesca,
2019
; Cappa et al.,
2019
; Guglielmi et al.,
2015
). Nonetheless, we find that multiple frictional scenar
-
ios of varying spatial behavior and proneness to large earthquakes match the slip observations of the field
LAROCHELLE ET AL.
10.1029/2020GL091188
2 of 11
Figure 1.
In situ measurement and modeling of fault slip induced by fluid injection. (a) Schematic of the field
experiment presented in Guglielmi et al. (
2015
) in which fluid injected into a borehole crossing a natural but inactive
fault caused its reactivation. A special borehole probe (SIMFIP) was used to measure the fault displacements directly
at the injection site. (b) Pressure, flow rate, and fault slip measured during the field experiment. The colored lines
and associated parameters correspond to the three different hydrological models considered in this study. The gray
area indicates the depressurization stage that has not been shown nor modeled in prior studies. (c) Schematic of
the model used to simulate slip on a fault plane embedded in an elastic bulk medium. Snapshots of a sample fluid
pressure diffusion scenario and its resulting fault slip are shown for illustration (the darker colors indicate later times).
Schematics (a) and (c) are not to scale.
Geophysical Research Letters
experiment equally well during fault pressurization. We also find that depressurization provides further
constraints that could help identify potentially hazardous faults.
2.
Data and Methods
2.1.
A Unique Fluid-Injection Experiment on a Natural Fault
The unprecedented field experiment involved injecting water directly into the fault zone and measuring the
resulting fault slip at a depth of 280 m with a specially designed borehole probe (Guglielmi et al.,
2015
) (Fig
-
ure
1a
). Prior to the experiment, the shear and normal stress acting on the fault were estimated at 1.65 ± 0.5
and 4.25 ± 0.5 MPa, and the permeability and bulk modulus of the initially dry fault at 7 × 10
−12
m
2
and
13.5 ± 3.5 GPa, respectively. Figure
1b
summarizes the main observations of the experiment, including the
deceleration of slip associated with depressurization not discussed in previous works. The slip measured
during the pressurization phase displays three distinct slip stages. At first, the fault is inactive and no signif
-
icant slip is recorded. The second stage initiates between 300 and 400 s, when slip rates attain
∼
10
−7
m/s and
the accumulated slip becomes measurable within the timeframe of the experiment. Stage 3 corresponds to
the sharp acceleration to slip velocities of
∼
10
−6
m/s without any significant increase in injection pressure
at
∼
1200 s. Hydromechanical modeling suggests that 70% of the 20-fold increase in permeability during
the experiment occurred prior to this acceleration (Guglielmi et al.,
2015
). Laboratory experiments were
also performed on grinded materials from the fault zone to further constrain the rate-and-state frictional
properties (Cappa et al.,
2019
).
2.2.
Diffusion of Pore Fluid Pressure Into the Fault Zone
We model the field experiment as a fluid injection into a planar fault embedded in an elastic medium (Fig
-
ures
1a
and
1c
). We simulate the fluid injection by prescribing an evolution of pore pressure at the center of
the fault that approximates the pressure history of the field experiment (Figure
1b
, top). Simulations with
a smooth pressure evolution result in similar but easier to interpret simulation results than those with the
exact pressure history (Figures
S1–S2
).
The imposed pressure diffuses axisymmetrically into the fault plane as follows:
2
2
, ,,
1
prtprtprt
t
rr
r
(1)
where
p
is the pore pressure,
r
is radial distance,
t
is time, and
is the hydraulic diffusivity. The diffusion
is numerically implemented using a forward finite difference scheme. Injection pressure is prescribed at a
distance of
inj
r
= 0.05 m from the center of the fault to mimic the experimental procedure. Although we
prescribe zero pressure boundary conditions far from the injection point to emulate the initially dry fault,
the choice of boundary condition is not essential here because the size of the simulated fault (250 m) is
larger than that of the pressure diffusion. Models with larger fault domains produce nearly identical results
(Figure
S3
).
Although both pressure and flow rate are reported as part of the field experiment, the exact value of the
hydraulic diffusivity
is still uncertain because the spatial extent of the pressurized zone and the fault
thickness over which the diffusion occurs,
b
, are poorly constrained. The volumetric flow rate,
Q
, depends
on
and
b
as:
inj
inj
2
,
s
rbS
p
Qt
rt
gr
(2)
where
S
S
is the specific storage and
the density of water. Hence, for a given flow rate, there is a trade-off
between the fault thickness
b
over which the fluid diffusion occurs, the hydraulic diffusivity
and the
specific storage
S
S
of the fault zone (and hence permeability
S
S
k
g
where
is the dynamic viscosity of
water). Note that
affects
Qt
via both the prefactor and the
inj
,
p
rt
r
term in Equation
2
. In Section
3
,
LAROCHELLE ET AL.
10.1029/2020GL091188
3 of 11
Geophysical Research Letters
we use hydraulic diffusivities of 0.04, 0.20, and 0.85 m
2
/s to match field experimental measurements of
slip for different friction regimes. Assuming a specific storage of
S
S
2 × 10
−4
m
−1
as in Bhattacharya and
Viesca (
2019
), for example, these hydraulic diffusivities correspond to permeability values of 0.8, 4, and
17 × 10
−12
m
2
that are within the ranges presented in previous studies that considered permeability en
-
hancement: 0.8−1.3 × 10
−12
m
2
(Bhattacharya & Viesca,
2019
) and 7–100 × 10
−12
m
2
(Guglielmi et al.,
2015
).
These permeability values are also consistent with the flow rates measured in the field experiment, for
reasonable values of the fault thickness
b
of 29, 6.7, and 1.8 cm, respectively (Figure
1b
). While consider
-
ing permeability enhancement may be necessary to match the finer features of the pressure and flow rate
histories (unless the fault thickness
b
affected by fluid flow varies with time or with space), all three com
-
binations of the parameters we use reproduce the hydrologic observations to the first order. We, therefore,
consider a range of constant hydraulic diffusivity (and hence permeability) values in our search for models
that reproduce the main features of the experimental observations.
2.3.
Numerical Modeling of Fluid-Induced Fault Slip
As fluid pressure increases and diffuses into the fault plane, fault friction eventually decreases and meas
-
urable slip ensues (Figure
1c
). We model this induced fault slip using a fully dynamic 2D antiplane bound
-
ary integral method capable of simulating the complete seismic cycle including both aseismic and seismic
deformation (Lapusta et al.,
2000
; Noda & Lapusta,
2013
). Fault slip is governed by the following elastody
-
namic equation:
ini
, ,
,,
2
s
xtfpxtFxtVxt
c
(3)
where
is the shear stress,
f
the friction coefficient,
the normal stress,
ini
the initial (i.e., background)
shear stress,
F
a linear functional which depends on the slip history,
,
the shear modulus of the elastic
medium,
s
c
the shear wave speed, and
V
the slip rate. The friction coefficient in Equation
3
follows an em
-
pirical rate-and-state formulation derived from laboratory experiments which describes the dependence of
f
on the slip rate and a state variable
(Dieterich,
1979
,
2007
; Ruina,
1983
):
RS
,
ln
ln
VV
fVfab
D
V
(4)
where
a
and
b
are the direct and evolutionary rate-and-state parameters,
RS
D
is the critical slip distance
and
f
is a reference coefficient of friction at reference slip rate
V
. The state variable is assumed to evolve
according to the aging law (Marone,
1998
; Ruina,
1983
).
As the fault in the experiment is inactive prior to the fluid stimulation, the modeled fault is not loaded
tectonically. Fault slip is thus purely fluid induced, that is, no significant slip would occur without the in
-
jection within the time scales considered in the simulations. To initialize the models, we impose shear and
normal stresses in agreement with the values reported at the field site prior to the experiment (Guglielmi
et al.,
2015
) and initial state variable values consistent with a dormant, highly healed fault (Text
S1
; Fig
-
ures
S4–S7
). The corresponding initial slip rate is then computed from Equation
4
.
3.
Results
3.1.
Models in Agreement With the Slip Observations During Pressurization
By first limiting our analysis to the pressurization stage of the experiment (up to 1400 s), we find that the
observations are equally well reproduced by a family of models. Three representative cases, which we de
-
note lower, intermediate and higher friction models, are shown in Figures
2a–2c
and
S8–S11
and Table
S1
.
Below, we explain how we constrained these models by examining how the various parameters govern
the transitions between the different slip stages and considering the trade-off between friction and fluid
pressure.
LAROCHELLE ET AL.
10.1029/2020GL091188
4 of 11
Geophysical Research Letters
At the beginning of all simulations, slip rates are low and both inertial effects and elastic stress transfers are
negligible. Equation
3
is then reduced to:
ini
,,
fVpxt
(5)
As
p
increases and
ini
remains constant over time,
f
must increase via growing slip rates in order for Equa
-
tion
5
to remain true, resulting in a balance between the direct frictional effect and changes in pore pressure
(Dublanchet,
2019
). Slip rate and friction continue increasing until slip becomes significant at
V
∼
10
−7
m/s.
The onset of significant slip thus approximately coincides with the maximum friction reached during the
simulations (Figures
2a
,
2b
, and
S8
). The peak friction,
p
f
, can be approximated as:
LAROCHELLE ET AL.
10.1029/2020GL091188
5 of 11
Figure 2.
Multiple simulated scenarios match the pressurization stage of the experiment but respond differently
to depressurization. (a) Temporal evolution of pore fluid pressure, slip, and slip rate for three model scenarios (solid
curves) that reproduce the observations (black dots) during the field-experiment pressurization. (b) Simulated evolution
of friction with slip at the injection site; the three scenarios correspond to lower (red), intermediate (green), and higher
(blue) residual friction in comparison to the fault prestress (black dashed line). Note that only the intermediate and
higher friction faults result in slip consistent with the depressurization part. (c) Key frictional and hydraulic properties
of the three scenarios. (d) Similar to (a) but for an improved depressurization: Reducing injection pressure once slip
starts to accelerate would allow to distinguish between all three cases, helping to constrain the fault friction properties.
Geophysical Research Letters
ini
RS
ln
ln
p
s
VV
ffab
D
V
(6)
where
s
V
= 10
−7
m/s. The state variable remains at its initial value,
ini
, as it has not evolved significantly
yet due to negligible slip and short healing time compared to its large initial value. Moreover, because the
fluid pressure at the injection site is known at all times, we can relate
p
f
to the timing of slip initiation,
s
t
:
ini
0,
p
s
f
pt
(7)
It is thus possible to control
s
t
by computing the corresponding
p
f
with Equation
7
, and selecting
f
,
a
,
b
,
ini
, and
RS
D
such that Equation
6
is satisfied. The three example models have
s
t
between 300 and 400 s and
p
f
between 0.84 and 0.99 (Figures
2b
and
S8
).
Once significant slip starts accumulating, the fault begins weakening, until it reaches steady state and fric
-
tion reaches its quasi-static residual value of
ln /
r
ffabVV
at the latest stage of the fault pres
-
surization experiment (Figures
2b
and
S8
). As in Dublanchet (
2019
)'s rate-strengthening models, we find
that this transition to steady state is accompanied with a marked acceleration in slip rate (Phase II in Dub
-
lanchet,
2019
) which we assume to explain the acceleration observed at 1200 s.
The critical slip distance,
c
, over which friction weakens from
p
f
to
r
f
can be approximated as:
RS
/
pr
c
ff
bD
(8)
since
RS
fb
D
. Furthermore, from elasticity, slip is related to stress drop by:
Δ
Δ
h
(9)
where
h
is the length of the slipping zone. By equating Equations
8
and
9
at the center of the fault, we can
estimate the slipping zone size,
ac
h
, at which steady state is reached and Stage 3 initiates:
RS
ac
Δ
pr
D
ff
h
b
(10)
Moreover, by choosing
V
to be on the same order of magnitude as the fastest slip rate measured during the
field experiment (
V
= 10
−6
m/s), we can approximate
r
f
with
f
, since the contribution of
ln /
abVV
becomes small compared to that of
f
. Equation
10
can then be rewritten in terms of known parameters as:
ini
RS
RS
ac
ini
ac
ln
ln
0,
s
VV
ab
D
D
V
h
b
fpt
(11)
where
ac
t
denotes the onset of Stage 3. For all the simulations presented in this work, we find that adding a
prefactor of 3 to Equation
11
provides a good estimate of the slipping zone size at
ac
t
(Text
S2
). Remarkably,
ac
h
only depends on quantities at the injection site. We can thus control the initiation of Stage 3 in our sim
-
ulations by tuning the model parameters such that the slipping zone reaches length
ac
h
at
∼
1200 s, as is the
case for our three representative models in Figure
3
.
Another critical aspect in these simulations is the balance between friction and the pore pressure forcing.
Figures
S20–S23
illustrate how the aseismic slip zone grows with decreasing
f
and increasing
, respec
-
tively. In particular, during Stage 3, the spatial extent of the slipping zone with respect to the pressurized
LAROCHELLE ET AL.
10.1029/2020GL091188
6 of 11
Geophysical Research Letters
zone as well as the slip rate at the injection site depend on the difference between the residual and initial
friction,
ini
r
ff
, which controls the elastic energy available to drive fault rupture once initiated (Bhattacha
-
rya & Viesca,
2019
; Dublanchet,
2019
; Galis et al.,
2017
; Garagash & Germanovich,
2012
) (Figures
S19a–
S19c
). Note that this is distinct from the difference between peak and initial friction,
ini
p
ff
(e.g., Gis
-
chig,
2015
), which controls the timing of fault reactivation as discussed above.
Given all these considerations, for each diffusion scenario presented in Figure
1b
, we find a corresponding
frictional model by adjusting
f
such that the simulated slip matches the observations during the first two
slip stages and produces a sufficiently large slip transient during Stage 3. To be able to use
r
f
(and hence
f
) values in agreement with the range 0.55–0.65 inferred from laboratory experiments on the grinded fault
zone material (Cappa et al.,
2019
), we set
ini
f
to 0.54 (
ini
= 2.15 MPa,
= 4.00 MPa), which is within the
uncertainty range of the initial stress measurements. A smaller
ini
f
would require smaller values of
f
out
-
side of this range to obtain the same slip at the injection site. Moreover, the selected values of
f
restrict the
range of possible values for the term
ini
RS
ln
/
bVD
in Equation
6
in order for slip to initiate between 300
and 400 s, which in turn restricts factor
RS
/
Db
in Equation
11
in order for Stage 3 to initiate at 1200 s.
The factor
RS
D
, which appears in estimates of critical nucleation lengths also needs to be large enough to
avoid nucleation of dynamic events within the experimental time (e.g., Rice & Ruina,
1983
; Rubin & Am
-
puero,
2005
). Finally, we fine tune parameters
a
and
ini
to adjust the slope and timing of the acceleration,
respectively. Note that decreasing
a
, while keeping
b
constant increases the slope of the slip acceleration,
due to the (weak) dependence of
r
f
on
ab
, and eventually leads to the nucleation of a dynamic event
right at
ac
t
(Figure
S16
and
S19d–S19f
). This procedure results in a family of models with
f
= 0.48–0.60,
ab
= −0.001–−0.005 (
b
= 0.016),
ini
= 1.2 × 10
12
–7.0 × 10
12
s and
= 0.04–0.85 m
2
/s that match the slip
observations equally well during pressurization.
Although the three models exhibit comparable slip histories at the injection site, they differ in features that
were not directly accessible to field observation. In particular, their spatial behaviors differ qualitatively
(Figure
3
,
S9–S11
). Defining the pressurized zone with 0.5 MPa pressure contours as in previous works, the
lower friction scenario produces an aseismic front that outruns the pressurized region, within 1400 s, as in
slip-weakening models (Bhattacharya & Viesca,
2019
) (Figure
3d
). By contrast, in the higher friction model,
which reproduces the observations equally well, aseismic slip remains confined well within the pressur
-
ized area (Figure
3f
). Our models demonstrate that slip did not necessarily extend beyond the pressure
perturbation during the experiment; that explaining a slip history at a single point in space is a nonunique
problem; and that further hydromechanical complexity is not required to explain the observed slip to first
LAROCHELLE ET AL.
10.1029/2020GL091188
7 of 11
Figure 3.
Whether the slipping zone is contained within or outruns the pressurized zone depends on fault friction.
Spatial and temporal evolution of (a–c) pore fluid pressure and (d–f) slip rate for the three scenarios of Figure
2
during
pressurization. The purple line shows the estimate
ac
h
of the slipping zone for the acceleration stage. Black dashed
lines indicate the extent of the pressurized zone defined by 0.5 MPa fluid pressure contours. During the pressurization
stage, the slipping zone of the lower-friction case outruns the pressurized zone while the intermediate and higher
friction cases remain confined to the pressurized zone.