manuscript submitted to
JGR: Solid Earth
A spectral boundary-integral method for faults and
1
fractures in a poroelastic solid: Simulations of a
2
rate-and-state fault with dilatancy, compaction, and
3
fluid injection
4
El ́ıas Rafn Heimisson
1
,
2
, Shengduo Liu
3
, Nadia Lapusta
2
,
3
, John Rudnicki
4
5
1
Swiss Seismological Service, ETH Zurich, Zurich, Switzerland
6
2
Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USA
7
3
Division of Engineering and Applied Science, California Institute of Technology, Pasadena, California,
8
USA
9
4
Department of Civil and Environmental Engineering and Department of Mechanical Engineering,
10
Northwestern University, Evanston, IL, USA
11
Key Points:
12
•
We present a novel spectral boundary-integral method (SBIM) for a 2D poroe-
13
lastic solid
14
•
We solve for fault slip with fully-coupled dilatancy and injection on a rate-and-
15
state fault
16
•
The method is applied to study the influence of several parameters on fault sta-
17
bility
18
Corresponding author: El ́ıas Rafn Heimisson,
elias.heimisson@sed.ethz.ch
–1–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
Abstract
19
Fluid-fault interactions result in many two-way coupled processes across a range of length
20
scales, from the micron scale of the shear zone to the kilometer scale of the slip patch.
21
The scale separation and complex coupling render fluid-fault interactions challenging to
22
simulate and may ultimately limit our understanding of experimental data and induced
23
seismicity. Here we present spectral boundary-integral solutions for in-plane interface
24
sliding and opening in a poroelastic solid. We solve for fault slip in the presence of rate-
25
and-state frictional properties, inelastic dilatancy, injection, and the coupling of a shear
26
zone and a diffusive poroelastic bulk. The shear localization zone is treated as having
27
a finite-width and non-constant pore pressure, albeit with a simplified mathematical rep-
28
resentation. The dimension of the 2D plane strain problem is reduced to a 1D problem
29
resulting in increased computational efficiency and incorporation of small-scale shear-
30
zone physics into the boundary conditions. We apply the method to data from a fault
31
injection experiment that has been previously studied with modeling. We explore the
32
influence of inelastic dilatancy, bulk poroelastic response, and bulk diffusivity on the sim-
33
ulated fault slip due to the injection. Dilatancy not only alters drastically the stability
34
of fault slip but also the nature of pore pressure evolution on the fault, causing signif-
35
icant deviation from the standard square-root-of-time diffusion. More surprisingly, vary-
36
ing the bulk’s poroelastic response (by using different values of the undrained Poisson’s
37
ratio) and bulk hydraulic diffusivity can be as critical in determining rupture stability
38
as the inelastic dilatancy.
39
Plain Language Summary
40
Earthquakes occur on faults deep in the Earth’s crust. At this depth, the faults are
41
surrounded by rock and water that fills up pores and fractures in the rock. This water
42
affects how the surrounding crust responds to earthquakes or slip on the faults. Water
43
also plays an important role within the faults since it will decrease or increase the fric-
44
tional resistance if it causes pressurization or depressurization, respectively. A common
45
cause of pressurization in faults is by an injection of fluid, which is done for many dif-
46
ferent purposes ranging from geothermal exploitation, carbon sequestration, or waste-
47
water disposal. Here we develop a new efficient method to simulate fault slip and earth-
48
quakes in a porous and fluid-filled medium. This allows us to better understand the role
49
of water in earthquake processes, either in the medium surrounding the fault or within
50
–2–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
the fault. We compare our method to a previously studied experiment where water was
51
injected directly into a fault and slip measured. In addition, we investigate certain phys-
52
ical properties of the porous rock that have not received much attention in the litera-
53
ture. We find that they significantly influence if earthquakes occur due to injection.
54
1 Introduction
55
The role of fluids in seismic and aseismic faulting processes has been of significant
56
interest in the last few years. Mounting evidence indicates that fluids may play an im-
57
portant role in a diverse set of mechanisms that alter fault slip behavior ranging from
58
earthquake triggering to slow slip events.
59
The most prominent example of fluid and fault interactions is the clear link between
60
fluid injection and induced seismicity, as originally pointed out by Raleigh et al. (1976);
61
Hsieh and Bredehoeft (1981) and remains a critical issue (e.g. Ellsworth, 2013). This phe-
62
nomenon has a straightforward mechanical explanation: higher pore pressures, due to
63
injection, reduce the effective normal stress and thus the frictional resistance of the fault.
64
The fault then slips faster and may accelerate the generation of seismic instabilities. This
65
problem has been frequently modeled with a straightforward implementation of one-way
66
coupling of pore pressure and frictional strength where pore pressure perturbations are
67
imposed and slip or number of seismic events are computed. Injection into faults may
68
lead to sustained aseismic transients (e.g. Viesca & Dublanchet, 2019; Bhattacharya &
69
Viesca, 2019), which may later become seismic events depending on the frictional prop-
70
erties of the fault (Larochelle et al., 2021a). A more detailed investigation of this prob-
71
lem reveals considerable complexity in pore pressure evolution if heterogeneous perme-
72
ability structures and poroelasticity are considered (e.g. Yehya et al., 2018)
73
The poroelastic properties of the crust have lately been receiving more interest, most
74
prominently as a long-ranging and fast-acting mechanism in which faults can be stressed
75
due to injection or extraction (Segall & Lu, 2015). However, there is also significant lit-
76
erature on the role of poroelasticity in influencing the nucleation or propagation of seis-
77
mic and aseismic ruptures (Rudnicki & Koutsibelas, 1991; Dunham & Rice, 2008; Jha
78
& Juanes, 2014; Heimisson et al., 2019, 2021). An effect of particular importance in re-
79
gard to the influence of poroelasticity is that, during in-plane sliding, compression and
80
dilation of the host rock induces pore pressure change in the shear zone (Heimisson et
81
–3–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
al., 2019, 2021); this effect is discussed further in section 1.1. Thus the poroelastic re-
82
sponse of the bulk, induced by an ongoing rupture, may influence the effective normal
83
stress and hence shear resistance to the rupture, creating a feedback loop. Poroelastic-
84
ity also influences and introduces a diffusion-dependent time-evolving shear stress on the
85
fault plane with significant implications for the stability of sliding (Heimisson et al., 2021).
86
Processes other than poroelasticity may change pore pressure in an active shear
87
zone and affect rupture and instability formation on faults. The generation of aseismic
88
slip transients on faults is believed to be related to pore fluids. For example, transient
89
slow slip events (SSEs) in subduction zones are thought to be related to high pore pres-
90
sure conditions (e.g., Liu & Rice, 2007; B ̈urgmann, 2018). A primary challenge in ex-
91
plaining the mechanics of transient slow slip is to understand why it starts, but does not
92
become an earthquake. One potential mechanism is a geometric restriction, in which the
93
high-pore-pressure region is large enough to cause slip acceleration, for example, due to
94
rate-and-state velocity-weakening friction properties, but too small for that slip to be-
95
come seismic (Liu & Rice, 2005, 2007). Another potential explanation is the change from
96
velocity-weakening to velocity-strengthening friction with increasing slip rates (Shibazaki
97
& Shimamoto, 2007; Hawthorne & Rubin, 2013; Leeman et al., 2016). Rate-and-state
98
faults with velocity-strengthening friction and additional destabilizing effects can also
99
produce SSEs in models with poroelasticity (Heimisson et al., 2019) and viscoplasticity
100
(Tong & Lavier, 2018). Inelastic dilatancy of granular fault gouge, which can lead to a
101
reduction in pore pressure and stabilize fault slip, has been highlighted as a naturally
102
present fluid-related mechanism that can explain how slow slip transients do not evolve
103
into seismic events (e.g. Segall & Rice, 1995; Segall et al., 2010). Modeling of fault slip
104
with inelastic dilatancy can explain many properties of slow slip events, including their
105
scaling (Dal Zilio et al., 2020).
106
Multiple mechanisms may act at a time. Recently, numerical simulations have started
107
exploring the simultaneous injection and inelastic dilatancy in a diffusive shear zone (Ciardo
108
& Lecampion, 2019; Yang & Dunham, 2021). However, these efforts have been limited
109
to a non-diffusive and elastic bulk. Coupling with a poroelastic bulk introduces another
110
degree of complexity, where elastic dilation and compression of the bulk generate pore
111
pressure transients. Further complexity is introduced by field observations indicating that
112
permeability of the shear zone in a fault core may be very different from the surround-
113
ing damage zone and host rock (e.g. Wibberley & Shimamoto, 2003). Further, the shear-
114
–4–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
ing of gouge material can dramatically reduce the permeability perpendicular to the shear-
115
ing direction and thus result in the shear zone having a significantly anisotropic perme-
116
ability (Zhang et al., 1999).
117
Here we present a spectral boundary-integral method that allows us to simulate
118
quasi-dynamic slow and fast slip on a rate-and-state fault with dilatancy/compaction
119
and fluid flow in a plane-strain poroelastic medium. We take a boundary layer approach
120
where the outer solution, which is the spectral representation of the poroelastic bulk, treats
121
the fault as a zero-thickness interface with suitable boundary conditions. However, the
122
inner solution considers the fault to be a finite-width shear zone. We consider the fric-
123
tional properties of the shear zone to be determined by its width-averaged properties.
124
The bulk is an isotropic standard quasi-static Biot poroelastic solid with a hydraulic dif-
125
fusivity
c
. The shear zone has frictional strength described by rate-and-state friction, with
126
inelastic state-dependent dilatancy and compaction and anisotropic permeability: the
127
permeability across the shear zone is different than the permeability along the shear zone.
128
The inelastic state-dependent dilatancy and compaction are implemented using the Segall
129
and Rice (1995) approach, as explained later. We frequently refer to this process only
130
as ”dilatancy” for the sake of brevity, and that is also how it is commonly referred to
131
in the fault mechanics community. However, we remind the reader that the ”dilatancy”
132
law also predicts compaction under certain conditions. The pore pressure in the layer
133
is simplified and assumed to be bi-linear where the two linear profiles are continuous at
134
the center of the shear zone (as in Heimisson et al., 2021, see also section 1.1). The spec-
135
tral representation uses analytical convolution kernels, which are truncated for efficiency
136
similar to Lapusta et al. (2000), but at time scales relevant for the bulk diffusion at the
137
specific wavenumber.
138
When slip speed becomes high enough in a narrow enough shear layer with small
139
enough permeability, then thermal pressurization of pore fluids due to shear heating may
140
also become important (e.g. Rice, 2006; Bizzarri & Cocco, 2006). While such effects may
141
be critical for seismic rupture evolution (e.g. Noda & Lapusta, 2013), they may be neg-
142
ligible or at least much less pronounced in the nucleation phases of the seismic cycle (Segall
143
& Rice, 2006; Segall, 2010), which are primarily the focus of this study. Consequently,
144
we do not account for thermal pressurization.
145
–5–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
The paper first discusses the general problem setup (section 1.1). For complete-
146
ness, there is a quick review of governing equations and boundary conditions (section
147
2). However, we highlight that a more complete description is found in Heimisson et al.
148
(2021) with the exception of added complexity introduced into the fluid mass balance
149
(section 2.3.1) not included in previous work. In section 3, we provide the analytical spec-
150
tral boundary-integral solutions for sliding and opening of an interface in a plane-strain
151
poroelastic solid. The numerical approach taken to solve the coupled problem - with di-
152
latancy, compaction, and injection in a poroelastic solid - is described in section 4. Fi-
153
nally, we show an application of the method (section 5), where we use constraints from
154
a field experiment (Guglielmi et al., 2015) and a recent numerical study that modeled
155
the field experiment data (Larochelle et al., 2021a). Finally, we discuss the role of poroe-
156
lasticity, and other fluid-based mechanisms, in the dynamics of injection-induced seis-
157
mic and aseismic slip.
158
1.1 Problem description
159
The general problem setup can be divided into three domains. Two are isotropic
160
poroelastic half-spaces, which we call the bulk, one in
y >
0 region and the other in
y <
161
0 region. The third is a shear zone made from fault gouge, which separates the two half-
162
spaces (Figure 1a). The two poroelastic half-spaces are assumed to have the same ma-
163
terial properties, which we characterize through the shear modulus
G
, Skempton’s co-
164
efficient
B
, drained Poisson’s ratio
ν
, undrained Poisson’s ratio
ν
u
, and hydraulic dif-
165
fusivity
c
(e.g., Cheng, 2016; Detournay & Cheng, 1995; Rice & Cleary, 1976). In some
166
cases, other poroelastic parameters may be displayed for compactness, legibility, and in-
167
tuition. However, the implementation of the method we present uses the aforementioned
168
five.
169
The shear zone is a thin layer of half-width
. Here thin indicates that
should be
170
much smaller than any significant variation in fields, such as slip or pressure, along the
171
x
−
axis
, which is fundamental for accuracy of the boundary-layer treatment of the shear
172
zone. The properties of the shear zone or fault gouge are characterized by reference poros-
173
ity
φ
0
, inelastic dilatancy coefficient
γ
(Segall & Rice, 1995), and pore-pressure and normal-
174
stress dependent void-volume compressibilities
β
p
n
and
β
σ
n
. In addition, the intact gouge
175
material compressibilities are
β
p
g
and
β
σ
g
, and the fluid compressibilities are
β
p
f
and
β
σ
f
.
176
The frictional strength of the shear zone is determined by the reference coefficient of fric-
177
–6–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
tion
f
0
, the characteristic state evolution distance
L
, the constitutive parameter
a
that
178
scales the direct rate dependence of friction, and the constitutive parameter
b
that scales
179
the state dependence of friction. These parameters and properties of the shear zone are
180
the same as in Heimisson et al. (2021) where a more detailed discussion is offered. We
181
also note that their meaning is presented in the context of the governing equations in
182
section 2. The hydraulic properties of the layer are somewhat different here compared
183
to Heimisson et al. (2021). First, we consider that there may be a source of fluid mass
184
in the layer, for example by injection, indicated by
Q
. Second, we include an anisotropic
185
mobility (permeability over dynamic fluid viscosity). In particular, the mobility in the
186
y
direction,
κ
cx
can be different from the mobility in the
x
direction
κ
cy
. Thus, fluids
187
injected into the fault have multiple migration paths, along the shear zone, perpendic-
188
ular to the shear zone, and in both
x
and
y
directions in the bulk. Furthermore, an in-
189
crease in pore pressure in the bulk can migrate into the shear zone and also into the bulk
190
on the other side. (Figure 1a)
191
–7–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
fl
ow directions
in shear zone
y
x
diffusion in bulk in
both x and y directions
shear zone /
fault gouge
Poroelastic bulk
Poroelastic bulk
slip across shear zone
0
y
y=0
Distance from shear zone center
compression
dilation
Distance along layer
Pore pressure change
x
0
Pore pressure change
a
b
panel b
0
Pore pressure change
dilatancy dominated
pressure pro
fi
le
slip dominated
pressure pro
fi
le
injection dominated
pressure pro
fi
le
isolated
Figure 1.
Schematic overview of the problems setup and possible pore pressure profiles
scenarios in the shear zone.
a
Injection occurs in a thin shear zone embedded between two poroe-
lastic solids of the same properties. This injection causes fluid migration along the shear zone,
across the shear zone, and into the bulk. The evolving pore fluid pressure leads to slip across the
shear zone.
b
Pore pressure profiles that can occur during the propagation of a single rupture
induced by injection. If the pore pressure diffusion is ahead of the rupture, then the shear zone
has increased pressure at the center
p
c
(right-most profile). Once considerable slip has occurred
the inelastic dilatancy may have reduced the pressure, even to the point of being less than hy-
drostatic, which we may call a dilatancy dominated pore pressure (left-most profile). Between
the two cases of an injection dominated regime and a dilatancy dominated regime we expect at
or near the rupture tip the two effects may cancel. However, the compression and dilation of the
host rock induced by the inhomogeneous slip can significantly change the pore pressures on either
side of the shear zone (
p
+
and
p
−
).
–8–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
A key question in induced seismicity is to understand when so-called runaway rup-
192
tures happen, that is ruptures that propagate well outside a pressurized region. This is
193
a useful focal point to explain some of the general dynamics that we expect from the de-
194
scribed problem above. When injection into a fault occurs, there are two important length
195
scales along the
x
dimension (Figure 1) that can interact and explain the dynamics of
196
the slip. First, how far the pressure front from the injection site has diffused, which we
197
can define as the region of significantly elevated pore pressure. Second, how far the rup-
198
ture tip has propagated, which can be understood as the region of significant fault slip.
199
If a fault has relatively low shear stress, i.e., its shear stress over initial effective normal
200
stress is significantly below its reference friction coefficient, or is well-healed, which may
201
be common in injection experiments, the pore pressure front controls how far the rup-
202
ture tip can move since the frictional resistance is too great outside the pressure front
203
(e.g., Larochelle et al., 2021a). However, if a fault is relatively well-stressed, or if the slip-
204
ping region enters a more well-stressed portion of the fault or a portion of the fault with
205
lower friction, then the rupture may become self-sustained and rupture outside the pres-
206
sure front. Thus the rupture may initially be contained by the pressure front, but evolve
207
to become a runaway rupture.
208
The interplay of the rupture tip and pressure front provides a useful qualitative ex-
209
planation of the transition from a confined to runaway rupture. However, additional com-
210
plexity, which is related to the pressure profile across the fault, plays an important role
211
in determining the if, when or how such a rupture can happen. If a rupture is initiated
212
in a shear zone by injection, the pressure profile across the shear zone (i.e. pressure change
213
with
y
, Figure 1
b
) can be dominated by different mechanisms depending on whether ob-
214
serving the profile at a
x
coordinate that is ahead of the rupture, at the tip or behind
215
the tip (Figure 1b ). This will be particularly prominent for an in-plane rupture direc-
216
tion due to the volumetric straining of the bulk. If the pressurized zone is ahead of the
217
rupture the shear zone central pressure (
p
c
) will be elevated. The pore pressures adja-
218
cent to the shear zone (
p
+
and
p
−
) will also be elevated due to the leak-off into the bulk.
219
Near the tip region, the influence of dilatancy has started to lower the pore pressure
p
c
,
220
but furthermore volumetric straining of the bulk has caused an increase in pore pressure
221
on the compressive side (
p
+
) and decrease on the dilating side (
p
−
) due to poroelastic
222
coupling. Finally, behind the tip dilatancy may have further reduced the pressure
p
c
and
223
possibly reversed the sign compared to the background equilibrium pressure and caused
224
–9–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
flow back into the shear zone. We thus suggest that in order to model rupture propa-
225
gation, earthquake nucleation, and understand runaway ruptures in a fluid-saturated medium
226
due to injection, we must consider coupling that arises from the interplay of several mech-
227
anisms that alter the pore pressure.
228
2 Governing equations
229
This section will describe the conservation laws, friction laws, and boundary con-
230
ditions. All the governing equations and boundary conditions with the exception of Sec-
231
tion 2.3.1, which describes the fluid mass balance, are the same as in Heimisson et al.
232
(2021). We state the equation with brief explanations for completeness, but refer the reader
233
to Heimisson et al. (2021) for more elaborate discussion and derivations.
234
2.1 Poroelastic Bulk
235
The quasi-static theory of poroelasticity can be described as four coupled partial
236
differential equations written in terms of displacements
u
i
and pressure changes
p
rel-
237
ative to an equilibrium pressure state (e.g., Detournay & Cheng, 1995; Cheng, 2016)
238
Gu
i,kk
+
G
1
−
2
ν
u
k,ki
=
αp
,i
(1)
and
239
1
M
p
,t
−
κp
,kk
=
−
αu
k,kt
,
(2)
where the material parameters are as follows:
G
: shear modulus,
ν
: drained Poisson’s
240
ratio,
α
: Biot-Willis parameter,
M
: Biot modulus, and
κ
is the mobility (the ratio be-
241
tween the permeability and fluid viscosity). In later expressions a different set of poroe-
242
lastic material parameter may be used for compactness and increased intuition.
243
In this work, we assume plane strain deformation, in which case the governing equa-
244
tions can be reduced to three. Further simplification and decoupling of the governing equa-
245
tions is possible by using the McNamee-Gibson displacement functions (McNamee & Gib-
246
son, 1960; Verruijt, 1971). In obtaining solutions to equations (1) and (2) we follow the
247
strategy explained in the Appendix of Heimisson et al. (2019) using the McNamee-Gibson
248
displacement functions but using the boundary conditions listed in the next section.
249
–10–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
2.1.1 Boundary conditions
250
Here we apply the same boundary conditions as in (Heimisson et al., 2021) at the
251
interface, i.e. the shear zone and at infinity.
252
lim
y
→
0
±
u
+
x
−
u
−
x
=
δ
x
,
(3)
lim
y
→
0
±
u
+
y
−
u
−
y
=
δ
y
,
(4)
lim
y
→±∞
u
±
x
= 0 and
u
±
y
= 0
,
(5)
lim
y
→±∞
p
±
= 0
,
(6)
lim
y
→
0
±
σ
+
xy
−
σ
−
xy
= 0
,
(7)
lim
y
→
0
±
σ
+
yy
−
σ
−
yy
= 0
,
(8)
where we have dropped the index notation and used
x
and
y
(as represented in Figure
253
1a).
254
The pore pressure in the shear zone is assumed to be bi-linear as in Heimisson et
255
al. (2021) This is a generalization of the leaky interface used in the plane strain dislo-
256
cation solution of Song and Rudnicki (2017). The pore pressure across the shear zone
257
is parameterized in terms of pressure at the center
p
c
at
y
= 0 and the pressure at the
258
shear zone boundaries where the poroelastic bulk meets the shear zone, that is,
p
±
at
259
y
=
±
. We can explicitly write out the assumed pore pressure profile as:
260
p
(
y
) =
y
(
p
+
−
p
c
) +
p
c
if 0
< y <
p
(
y
) =
y
(
p
c
−
p
−
) +
p
c
if
−
< y <
0
.
(9)
Thus equating the fluid mass flux into the shear zone and in the the bulk, and vice versa,
261
gives rise to a pressure gradient boundary condition:
262
dp
±
dy
∣
∣
∣
∣
y
=0
±
=
±
κ
cy
κ
(
p
±
−
p
c
)
,
(10)
where
κ
cy
is the shear zone mobility in the
y
direction and
κ
is the poroelastic bulk mo-
263
bility which is rated to the bulk hydraulic diffusivity by
c
=
Mκ
. We note that bound-
264
ary conditions for the bulk are applied at
y
= 0
±
but in the description of the shear
265
–11–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
zone we treat it as a finite layer with thickness between
y
=
±
. This is because we take
266
a boundary layer approach (similar to Appedix B of Rudnicki & Rice, 2006) where the
267
inner solution, the shear zone, is assumed to have a finite thickness. However, the outer
268
solution, the bulk, approximates the layer as having an infinitesimal thickness. Thus the
269
assumption that any variation along the length of the shear zone occurs over a length
270
scale much smaller than
is implicit. In other words, we always require that
k
1,
271
with
k
representing the wavenumber (inverse of a wavelength) of any field that varies
272
along the
x
-dimension.
273
2.2 Frictional properties
274
As in Heimisson et al. (2021) we represent the frictional strength of the layer in an
275
averaged sense.
276
Let us assume that the frictional strength of every point in the layer can be rep-
277
resented as follows:
278
τ
(
x,t
)
σ
(
x,t
)
−
p
(
x,y,t
)
=
f
(
x,y,t
)
for
−
< y < ,
(11)
where
τ
(
x,t
) is the sum of all contributions to the shear stress, both initial background
279
value and slip contributions. We note that the shear stress is assumed to be spatially con-
280
stant across the layer. Similarly,
σ
(
x,t
) represents background initial effective normal
281
stress (normal stress minus the ambient pore pressure) in addition to the slip induced
282
changes in normal stress and we assume is spatially constant across the layer. However,
283
we have separated from the description the perturbations in pore pressure
p
(
x,y,t
) since,
284
as previously discussed, they cannot be assumed to be constant in
y
. Using equation (9)
285
and averaging over the layer, we obtain:
286
τ
(
p
c
−
p
+
) log
(
σ
−
p
−
σ
−
p
c
)
+ (
p
c
−
p
−
) log
(
σ
−
p
+
σ
−
p
c
)
2(
p
c
−
p
−
)(
p
c
−
p
+
)
=
〈
f
〉
,
(12)
with the
〈
f
〉
representing the frictional coefficient of the layer. We have explored using
287
the equation above for modeling the interface frictional strength, but we find that it ren-
288
ders very similar results as an linearized approximation valid in the limit of the pore pres-
289
sure changes being small compared to the background normal stress:
290
–12–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
τ
= (
σ
−〈
p
(
t
)
〉
)
〈
f
〉
,
(13)
where
〈
p
(
t
)
〉
is the average pressure across the layers and can be computed directly
291
〈
p
〉
=
1
2
∫
−
p
(
y
)
dy
=
1
2
(
p
c
+
p
+
+
p
−
2
)
.
(14)
Equation (13) further offers a simpler interpretation of the role of the pore pressure in
292
the effective normal stress compared to equation (12), which helps in understanding the
293
simulation results.
294
We interpret the averaged friction coefficient
〈
f
〉
of the shear zone as being rep-
295
resented by the rate-and-state friction law (e.g., Dieterich, 1979; Ruina, 1983; Marone,
296
1998):
297
〈
f
〉
=
1
2
∫
−
f
(
x,y,t
)
dy
=
a
arcsinh
[
V
2
V
0
exp
(
f
0
+
b
log(
V
0
θ/L
)
a
)]
,
(15)
where we use the regularized form of the friction law that is also valid for slip speeds
V
298
much smaller than the reference slip speed
V
0
(Rice & Ben-Zion, 1996; Ben-Zion & Rice,
299
1997; Lapusta et al., 2000). Here
a
and
b
are constitutive parameters that describe the
300
rate dependence and state dependence of friction, respectively. Further,
f
0
is the refer-
301
ence coefficient and
L
is the characteristic slip distance over which the state evolves. The
302
state variable is described by the aging law (Ruina, 1983):
303
dθ
dt
= 1
−
θV
L
(16)
We note that here we have introduced a minor difference compared to (Heimisson
304
et al., 2021). We represent friction using the regularized friction law whereas the non-
305
regularized version was discussed by (Heimisson et al., 2021). In the linearized analy-
306
sis treated by Heimisson et al. (2021), there is no difference between the two versions.
307
2.3 Shear Zone
308
Here we analyze the fluid and solid constituent mass balance of the shear zone gouge.
309
This analysis is largely based on Heimisson et al. (2021) although here we introduce new
310
physical processes into fluid mass balance, which are detailed below. Heimisson et al. (2021)
311
–13–
ESSOAr | https://doi.org/10.1002/essoar.10510458.2 | CC_BY_NC_ND_4.0 | First posted online: Mon, 8 Aug 2022 05:12:30 | This content has not been peer reviewed.