of 22
Supporting Information for:
Unraveling scaling properties of slow-slip events
Luca Dal Zilio
1
,
2
, Nadia Lapusta
1
,
2
, Jean-Philippe Avouac
1
,
2
1
Division of Geological and Planetary Sciences, California Institute of Technology
2
Division of Engineering and Applied Sciences, California Institute of Technology
Contents of this file
1. Methods
2. Tables S1 and S2
3. Figures S1 to S5
4. Supporting References
Numerical methodology for simulating slow-slip events
For numerical simulations of long-term fault behavior, we use the numerical code BICyclE (Boundary
Integral Cycles of Earthquakes) (Lapusta et al., 2000; Lapusta & Liu, 2009; Noda & Lapusta, 2010;
Noda et al., 2013), which solves the coupled problem of elasto-dynamic equations of motion, fault
friction, and off-fault diffusion of heat and pore fluids. At each fault point during slip, dynamic
shear stress (or forcing) must be equal to friction (or fault shear strength). The resulting spontaneous
evolution of fault slip is a complex non-linear process governed by multiple feedbacks. Shear stress
at each fault point is affected by motion of all other fault points through wave-mediated dynamic
and eventually static stress changes. At the same time, frictional resistance itself is coupled to slip
rate and, through the state variable, to prior history of slip. Our computational approach employs an
efficient boundary-integral method in the Fourier domain, which allows us to resolve the slow tectonic
loading, slow-slip event nucleation and propagation — including possible radiation of seismic waves
— and the afterslip transient that can follow large slip events. An adaptive time stepping is used to
numerically resolve slow tectonic loading, slow-slip nucleation and propagation, and aseismic slip
transients.
We consider a megathrust fault segment embedded into an elastic medium, loaded by deep-seated
slip at the long term slip rate, and governed by rate-and-state friction (Dieterich, 1979; Dieterich et
al., 1981; Ruina, 1983; Marone, 1998; Beeler et al., 1996, 2008):
τ
=
f
σ
=
[
f
+
aln
V
V
+
bln
θV
D
RS
]
σ ,
(1)
dt
= 1
V θ
D
RS
,
(2)
where
τ
is the shear resistance,
V
is the instantaneous slip rate,
σ
=
σ
p
is the effective normal stress
accounting for pore fluid pressure
p
,
f
is the steady-state friction at the reference slip rate
V
,
a
,
b
, and
D
RS
are rate-and-state parameters, and
θ
is a state variable that describes the evolutionary effects in
response to changing slip rates. This law, commonly referred to as the Dieterich-Ruina or aging rate-
and-state friction, has experimental support (Beeler et al., 1996) and partial theoretical justification
(Rice et al., 2001). Other formulations for the evolution of the state variable exist, such as the slip
law (Ruina, 1983) as well as various composite laws, and the formulation that best describes various
laboratory experiments remains a topic of ongoing research (Bhattacharya et al., 2017). For a constant
slip velocity
V
, the state variable
θ
evolves to its steady-state value
θ
ss
=
D
RS
/V
, transforming the
shear resistance
τ
into its steady-state form
τ
ss
=
[
f
+ (
a
b
)
ln
V
V
]
σ ,
(3)
which highlights the importance of the non-dimensional parameter (
a
b
) of the rate dependence of
friction: stable slip for velocity-strengthening (VS) friction, (
a
b
)
>
0, and conditionally unstable
slip for velocity-weakening (VW) friction, (
a
b
)
<
0 (Rice & Ruina, 1983; Scholz, 1998).
Some of our models, including the reference model discussed in the main text, include the evolution
of pore pressure
p
due to inelastic shear dilatancy
φ
according to the following equation (Segall &
Rice, 1995; Segall et al., 2010; Segall & Rice, 2006):
∂p
∂t
=
α
hy
2
p
∂y
2
F
(
y
)
β
dt
,
(4)
where
α
hy
is the hydraulic diffusivity,
β
is the specific storage,
dφ/dt
is the rate of change in the
porosity due to inelastic deformation, and
F
(
y
)
is the Gaussian function representing the distribution
of the inelastic porosity over the thickness of the shear zone (
W
sz
). Following Segall and Rice (1995),
we associate the effects of dilatancy and compaction with the evolution of state variable within the
fault gouge, as motivated by experiments of Marone et al. (1990), so that:
δφ
=
 ln
(
V
0
θ
D
RS
)
,
(5)
dt
=

d
dt
ln
(
V
0
θ
D
RS
)
=

θ
dt
,
(6)
where

is a constant derived empirically (Table S1) (Marone et al., 1990), and
D
RS
is the character-
istic slip. Above steady state conditions, i.e. for
θ > D
RS
/V
,
θ
decreases (gouge dilatancy), while
below steady state,
θ
increases (gouge compaction). With this approach, slow-slip events nucleate
in areas of unstable friction under drained conditions but, as slip accelerates, dilatancy reduces pore
pressure quenching instability (Segall et al., 2010).
The thrust fault between a subducting oceanic plate and overlying continental plate is modeled as a
2-D planar interface embedded in a 3-D elastic half-space. The fault extends 320 km along the strike
direction and 60 km along the downdip direction. The planar fault is discretized by 1360 x 280 grid
cells along the strike and downdip directions, respectively. The cell sizes are thus 250 m along strike
and along downdip. A uniform tectonic loading at rate of 40 mm/yr is applied further downdip to
represent the convergence rate between the Juan de Fuca and North American plates in northern Cas-
cadia (Michel et al., 2018). We assign a velocity-weakening (VW) friction patch extending 300 km
along the strike direction (length
L
) and 25 km along the downdip direction (width
W
); elsewhere,
fault friction is assumed to be velocity strengthening (VS). The model setup with a narrow width is
designed to favor SSEs. The effective compressive normal stress
σ
is taken to be constant at 10 MPa;
this distribution of
σ
is appropriate for an over-pressurized crust at depth (Suppe, 2014; Audet et al.,
2009; Leeman et al., 2016). The corresponding elevated pore fluid pressures are required by all fric-
tion models to reproduce realistic properties of slow-slip events (Shibazaki & Iio, 2003; Liu & Rice,
2005, 2007; Shibazaki & Shimamoto, 2007; Rubin, 2008; Liu & Rubin, 2010; Matsuzawa et al., 2010;
Hawthorne & Rubin, 2013; D. Li & Liu, 2016; Luo & Ampuero, 2018; Viesca & Dublanchet, 2018);
high pore-fluid pressures also promote slow slip in models with bulk effects (Gao & Wang, 2017). As
SSEs are now found nearly universally at subduction zones (Bartlow et al., 2011; Kato & Nakagawa,
2014; Tsang et al., 2015; Fu et al., 2015; Michel et al., 2018; S. Li et al., 2016; Villegas-Lanza et al.,
2016; Houston et al., 2011), such high pore fluid pressures — and the associated low effective nor-
mal stresses — may be prevalent at most subduction zones, at least at the locations of slow-slip events.
Theoretical estimates of nucleation size and cohesive zone size
For 2D problems with velocity-weakening rate-and-state friction and no dilatancy effects, two theo-
retical estimates of the earthquake nucleation size are
h
rr
(Rice & Ruina, 1983) and
h
ra
(Rubin &
Ampuero, 2005):
h
rr
=
π
4
μ
D
RS
σ
(
b
a
)
;
(7)
h
ra
=
2
π
μ
D
RS
b
σ
(
b
a
)
2
;
(8)
where
μ
=
μ
for mode
III
,
μ
=
μ/
(1
ν
)
for mode II,
μ
is the shear modulus,
ν
is the Poisson’s
ratio, and the estimate (7) is valid for
0
.
5
< a/b <
1
, as explained in Rubin and Ampuero (2005). In
the parameter regime
0
< a/b <
0
.
3781
, that study advocates
b
instead of
(
b
a
)
2
/b
in the above
formula; however,
(
b
a
)
2
/b
b
as
a
0
, so that eq. 8 should approximately capture both cases.
Indeed, in our exploration of parameters for the occurrence of slow-slip events, as discussed in the
next section, we see no break in behavior between the two regimes of
a/b
.
In 3D problems, the nucleation size is given by Chen and Lapusta (2009):
h
= (
π
2
/
4)
h
ra
.
(9)
The theoretically estimated nucleation sizes in our models are
8.5–11 km based on eq. 7, and
29–31 km based on eq. 8. Table S1 summarizes the values of parameters used in the models. The
chosen values of
D
RS
are larger than those obtained in laboratory, which facilitate our numerical
computation (Lapusta & Liu, 2009).
Two important physical scales in the problem are the nucleation size
h
and cohesive zone size
Λ
.
The nucleation size is a crucial length scale during interseismic periods and
h
/
x
is an important
criterion to assess spatial resolution. The cohesive zone size is an important resolution criterion in
dynamic rupture because it gives the spatial length scale over which the shear stress drops from its
peak to residual values at the propagating rupture front (Palmer & Rice, 1973; Day et al., 2005). For
rate-and-state friction law,
Λ
0
, the size of
Λ
at zero rupture speed
v
r
, is given as:
Λ
0
=
C
1
μ
D
RS
.
(10)
where
C
1
is a constant and equal to
9
π/
32
if the stress traction distribution within the cohesive zone
is linear in space. Day et al. (2005) established that
Λ
0
/
x
of 3 to 5 is required to resolve dynamic
rupture. The ratio of nucleation zone size and cohesive zone size is given by
Λ
0
/
h
(
b
a
)
2
/
b
2
.
With the chosen values of
a
and
b
, this ratio in our models is between 0.3 and 0.5. Therefore, resolving
the cohesive zone is the more stringent numerical criterion here. In our models, we choose the cell
size small enough to resolve the cohesive zone with at least 30 cells. Hence the spatial discretization
in the simulations, with the cell size of
x
= 250 m, is small enough to resolve the evolution of stress
and slip rate.
The addition of dilatancy stabilizes the problem (Segall & Rice, 1995; Segall et al., 2010; Rubin,
2008; Liu & Rubin, 2010). For the parameters used (see Table S1), our SSEs nucleate in areas of
unstable friction under drained conditions and hence the nucleation size estimate
h
rr
is valid. This
occurs because the off-fault diffusivity is high enough that the pore pressure inside the shearing zone
equilibrates with the pore pressure outside rapidly enough for the dilatancy not to affect pore fluid
pressure.
Parameter regime for SSEs in our simulations
We extensively explore the effect of frictional parameters on the resulting behavior of our model (Fig.
S1). By increasing the characteristic slip
D
RS
(from 4 to 90 mm; Table S1), our models show four
different slip patterns: 1) fast (seismic) events, 2) slow-slip events, 3) decay oscillations, and 4) stable
creep (Fig. S1A). We also analyze the effects of the rate-and-state properties by performing experi-
ments with different
b
a
(from 0.004 to 0.01; Fig. S1B), and
(
b
a
)
/b
(from 0.28 to 0.71), which
give results in agreement with previous studies (Veedu & Barbot, 2016; Barbot, 2019). According
to our models (Fig. S1E), SSEs appear in a parameter regime between steady slip and regular (fast)
earthquakes.
In particular, the conditions that lead to SSEs for rate-and-state faults with no dilatancy effect are
given by:
1
3
W < h
rr
<
1
2
W
;
h
ra
> W
;
(11)
That is, the SSEs appear when the fault width
W
is smaller than the nucleation estimate
h
ra
but
larger, by a factor of 2 to 3, then the nucleation estimate
h
rr
. This condition makes qualitative sense
given the nature of the estimates. On the one hand,
h
rr
was obtained as the critical wavelength from
a linearized quasi-static stability analysis of steady sliding to perturbations of different wavelengths,
with the smaller wavelengths leading to decaying perturbations and larger wavelengths leading to
growing perturbations (Rice & Ruina, 1983). Hence
h
rr
is an estimate of the wavelengths that would
result in some slip instability, explaining why the SSEs spontaneously arise when
h
rr
is smaller than
W
. On the other hand,
h
ra
was obtained from the energy balance of the nucleation process that
takes the form of an expanding crack, giving the upper bound on the crack dimension (Rubin &
Ampuero, 2005). Previous studies have shown that
h
ra
is a good estimate of the nucleation zone at
the transition from quasi-static to seismic (weave-producing) slip (Chen & Lapusta, 2009). Hence
making
W
smaller than
h
ra
ensures that accelerating quasi-static slip does not transition into faster,
seismic earthquakes.
Dilatancy effects are then explored for models M1–M4 listed in Table S2. The models show that
dilatancy slows down the rupture velocity of SSEs, shifting the results towards relatively longer event
durations and thus enabling a wide time spectrum ranging from a few days to several weeks and
months (Ide et al., 2007).
Measurements of SSE duration and moment release
For each numerical experiment, we analyze all slow-slip events to determine their duration, along-
strike length, area, rupture velocity, and magnitude. The events are identified using a slip-velocity
threshold. With this approach, each SSE lasts while the slip velocity exceeds a threshold equal to
1
.
26
x
10
8
m/s (i.e., 10 times faster than the long-term slip rate). In order to assess possible un-
certainties in the duration and moment release of each event, we also include results for a standard
deviation of 20% (
1
.
26
x
10
8
m/s
±
2
.
5
x
10
9
m/s). The sensitivity of our results to the slip-
velocity threshold is tested by assuming a different definition for SSEs, with the velocity threshold
which is only 2 times faster than the long-term slip rate and a minimum moment release rate of 1
x
10
14
Nm/day. The resulting properties of SSEs have minimal variations (Fig. S4). The choice
of the slip-velocity threshold is to ensure that (1) selected SSEs have velocities exceeding the plate
convergence rate at the minimum and (2) it is around the lower resolution limit of GPS slip inversion
models in order to compare modeled SSE source parameters to those inferred from geodetic observa-
tions. For SSEs in Cascadia, Michel et al. (2018) found that the minimum resolved slip rate is 0.21
mm/day (i.e., 2 times the long-term slip rate) at the lower resolution limits for their network inversion
filter models, whereas Wech and Bartlow (2014) found a minimum resolved slip rate of
0.25 to 0.5
mm/day (2.2–4.5 times the long-term slip rate).
Measurements of stress drop
We consider the average energy-based stress drops (
σ
E
) (Noda et al., 2013), which has been re-
cently used to analyze earthquakes (Ye et al., 2016).
σ
E
corresponds to averaging the stress drop
distribution
σ
(
x
1
,x
3
)
the actual final slip
u
at each point as the weighting function:
σ
E
=
σ
u dS
u dS
,
(12)
where
S
is the planar fault area,
σ
is the initial shear stress minus the final shear stress at each rup-
tured location,
u
is the final slip distribution produced by the event. This measure is more convenient
than the more commonly used moment-based stress drop, as the moment-based stress drop requires
the determination of the shape of the ruptured area and computation of the appropriate weighting
function (47). Defining the ruptured area of a slow-slip event would require an additional criterion of
what constitutes significant slip, which is typically defined as a percentage of the maximum slip.
Parameter
Symbol
Value
Background characteristic slip distance
D
RS
5
-
90
mm
Rate-and-state properties in VW region
a
0
.
002
-
0
.
010
b
0
.
006
-
0
.
016
Rate-and-state properties in VS region
a
0
.
019
b
0
.
014
Effective normal stress
̄
σ
=
σ
n
p
10
MPa
Dilatancy coefficient

0;
3
·
10
4
Hydraulic diffusivity
α
hy
1
·
10
6
m
2
/s
Specific storage
β
5
·
10
11
1/Pa
Thickness of the shear zone
W
sz
1
·
10
3
m
Reference friction coefficient
f
0
.
6
Reference slip velocity
V
10
6
m/s
Poisson’s ratio
ν
0
.
25
Shear modulus
μ
30
GPa
Loading rate
V
pl
40
mm/yr
Spatial resolution
x
250
m
Width of VW region
W
25
km
Length of VW region
L
300
km
Table S1:
Model parameters used in this study.
Model
a
b
D
RS
(mm)
Dilatancy coefficient (

)
M1
0.004
0.014
45
/
M1-D*
0.004
0.014
45
3
·
10
4
M2
0.002
0.006
19
/
M3
0.0028
0.0088
27
/
M4
0.0032
0.0112
36
/
* The reference model shown in Fig. 2–4.
Table S2:
Parameters of the selected models, for which the scaling is shown in Fig. S4. M1 and M1-D have the same properties other
than dilatancy coefficient.