Plain Language Summary
Subduction zones host the world’s largest earthquakes, which exhibit
depth-dependent characteristics such as depletion of high-frequency
seismic radiation at shallow depth. Some recent studies propose that
elastic properties of wall rocks determine these features despite the
fact that fault friction has long been considered to be a dominant
factor. Here we design a suite of dynamic earthquake rupture models to
explore the roles of wall-rock properties and fault friction. We find
that fault friction plays more important roles than wall rock properties
in rupture propagation, slip distribution, and high-frequency depletion
at shallow depth. On the other hand, even though the low-velocity rock
layes decrease rupture velocity to some extent, they enhance and do not
reduce high-frequency radiation at shallow depth. We conclude that fault
friction plays more important roles than wall-rock properties in
depth-dependent rupture characteristics of subduction zone earthquakes.
1 Introduction
Subduction zones host the world’s largest earthquakes (Kanamori, 1986).
Subduction zone earthquakes exhibit depth-dependent seismic
characteristics, such as decrease in normalized source duration with an
increase in depth (e.g., Bilek & Lay, 1999), enhanced high frequency
radiation as depth increases (e.g., Rushing & Lay, 2012), slower
rupture velocity toward the trench for the shallow tsunami earthquakes
(Bilek & Lay, 2002), and heterogeneous coseismic slip distribution over
different depths (e.g., Ammon et al., 2005; Ide et al., 2011). In
particular, recent great earthquakes with wide along-dip rupture
extents, such as the 2004 Mw 9.1 Sumatra, 2010 Mw 8.8 Chile, and 2011 Mw
9.0 Tohoku earthquakes, show clearly depth-dependent variations in
frequency contents of seismic radiation and slip distribution.
High-frequency seismic radiations are imaged in the downdip portions of
the megathrusts by large seismic network back-projection methods (e.g.,
Lay et al., 2010; Kiser and Ishii, 2011; Koper et al., 2012; Ishii,
2011; Lay et al., 2012). Inversions of seismic, geodetic, and tsunami
data show that large slip with weak high-frequency seismic radiation
occurs in the updip portions of the megathrusts (e.g., Lay et al., 2010;
Tong et al., 2010; Ammon et al., 2011; Hayes, 2011; Ide et al., 2011).
These features in depth-varying rupture characteristics motivate a
four-domain conceptual model for megathrusts (Lay et al., 2012). These
four domains are A) near-trench domain where tsunami earthquakes with
very weak high-frequency radiation or aseismic slip occur, B) central
megathrust domain with large seismic slip and modest high-frequency
radiation, C) downdip domain with modest seismic slip and significant
coherent high-frequency seismic radiation, and D) transition domain
further downdip featuring slow-slip events, low frequency earthquakes
and seismic tremor. Although the increase in seismic velocities with
depth is recognized to likely cause increasing rupture velocity with
depth, the four domains are largely controlled by frictional properties
(including seismic, aseismic, and conditionally stable) in the
conceptual megathrust model (Lay et al., 2012).
Frictional properties on the plate interface control the wide spectrum
of slip behaviors (Scholz, 1998), with diverse observations of ordinary
earthquakes, low-frequency earthquakes, and tectonic tremor (Lay, 2015).
In a conceptually generic model (Bilek & Lay, 2002; Kodaira et al.,
2004; Scholz, 1998) of slip instability at subduction zones, the top
several kilometers are in a stable regime, where velocity-strengthening
fault conditions dominate. In the seismogenic zone, ranging from the
upper limit of ~4 km depth to the lower limit of
~35 km depth, unstable slip and velocity-weakening fault
conditions dominate. The downdip stable regime (> 35 km
depth) is mainly controlled by velocity-strengthening behaviors. This
conceptual model is further supported by experimental evidence, in which
frictional properties are depth-varying and temperature dependent
(Blanpied et al., 1995; den Hartog & Spiers, 2013), and has been widely
used in subduction earthquake simulations (e.g., Im et al., 2020; Liu &
Rice, 2005; Liu & Rice, 2007; Meng et al., 2022).
Recently, heterogeneous upper-plate properties are proposed to determine
depth-varying rupture behavior of megathrust earthquakes (Sallares and
Ranero, 2019; Sallares et al., 2021; Prada et al., 2021). Using 48
P-wave velocity (Vp) models obtained from wide-angel
reflection and refraction surveys across circum-Pacific and Indian Ocean
subduction zones, Sallares and Ranero (2019) develop a global model of
Vp(z). They average Vp at the lower part
of the upper plate as a function of interplate boundary depth below
seafloor (z) and calculate depth profiles of density ρ(z), S-wave
velocity Vs(z), and rigidity µ(z) with experiment-determined empirical
relationships of ρ (Vp) and Vs (Vp) (Brocher, 2005). They find that Vp
increases by a factor of 2.0-2.5 from ~3.0 km/s at 1 km
depth to ~6.5 km/s at 25 km depth, with decreasing
gradient downwards. They derive a depth profile of slip based on µ(z),
assuming the same rupture area for the same size of earthquakes at
different depths. Similarly, they obtain a depth profile of rupture
duration based on Vs(z), assuming rupture velocity being 70-90% of Vs.
Essentially, slip and rupture duration are inversely proportional to
rigidity and Vs, respectively, in their results. Sallares et al. (2021)
perform a site-specific study of the 1992 Mw 7.7 Nicaragua tsunami
earthquake. They obtain the upper-plate elastic properties from
wide-angel reflection and refraction seismic data and multichannel
seismic reflection across the rupture area. They also calculate the
moment release, slip and stress drop distributions of the earthquake
from a finite fault inversion. Consistent with Sallares and Ranero
(2019), they emphasize the dominant role of upper-plate elastic
properties in controlling large slip and long duration of the event at
shallow depth in this tsunami earthquake. Prada et al. (2021) perform 3D
dynamic rupture and tsunami simulations to explore the influence of
depth-varying upper-plate elastic properties (with a global model
developed by Sallares and Raneo, 2019) on rupture characteristics and
tsunamigenesis. They compare slip, rupture duration and frequency
content from different scenarios with different velocity structures.
They use a linear slip-weakening friction law with constant values of
friction parameters along the fault, including static and dynamic
frictional coefficients and the critical slip distance, for dynamic
rupture simulations. Therefore, their dynamic rupture models can be
considered as essentially having uniform friction properties along the
fault. Their models reproduce depth-varying rupture features in terms of
slip, rupture duration and frequency content that agree with Sallares
and Ranero (2019).
Both fault frictional properties and wall rock properties are deemed to
affect slip instability and dynamic rupture propagation. The
heterogeneous coseismic slip of great earthquakes (Mw\(\geq\) 8.0) along subduction zones appears to be more complicated than
that can be explained by a purely frictional or rigidity effect. For
instance, the largest coseismic slip can be concentrated near the trench
such as in the 2011 Mw 9.0 Tohoku-Oki earthquake (Ide et
al., 2011), while some subduction zones exhibit a rupture propagation
barrier near the trench such as the 2010 Mw 8.8 Maule
earthquake (Lin et al., 2013). This phenomenon implies combined effects
of spatial-varying frictional properties and wall rock rigidity. In the
classic spring-slider model, slip instability is controlled by the
sliding interface’s properties (i.e., fault frictional properties) and
stiffness of the spring (i.e., surrounding material properties)
(Dieterich, 1979; Rice & Ruina, 1983; Scholz, 1998). In dynamic rupture
simulations, many studies assume simple velocity structure such as a
homogeneous material, partly because they primarily explore effects of
heterogeneous friction and/or complex fault geometry on dynamic rupture
propagation. Effects of fault-bounding material properties are explored
by many other dynamic rupture modeling studies, most notably for
bimaterial problems (e.g., Harris and Day, 1997; Andrews and Ben-Zion,
1997; Duan, 2008a; Ampuro and Ben-Zion, 2008) and for low-velocity fault
zone problems (e.g., Harris and Day, 1997; Duan, 2008b; Huang and
Ampuero, 2011). Nevertheless, effects of heterogeneous velocity
structure at subduction zones, in particular shallow low-velocity
materials in the upper plate, need to be better understood and
incorporated into numerical models of shallow subduction zones. The
recent studies by Sallares and Ranero (2019), Sallares et al. (2021),
and Prada et al. (2021) make a significant contribution to this
endeavor. However, without contrasting and quantifying effects of the
two factors, namely fault friction and surrounding material property, in
one framework of physics-based models, it is difficult to ascertain
which of the two factors plays a more important role, among many other
factors such as nonplanar fault geometry and heterogeneous stress state.
There are some studies that incorporate both depth-varying frictional
properties and rigidity in dynamic rupture simulations for subduction
zone earthquakes. For example, Kozdon and Dunham (2013) perform 2D
dynamic rupture simulations of the 2011 Mw 9.0 Tohoku earthquake to
understand why and how the rupture could reach the trench. In their
models, heterogeneous velocity structure from seismic surveys of the
Japan trench (Miura et al., 2001, 2005) is included. They use a rate-
and state- dependent friction law (RSF) with velocity-weakening
frictional properties for the central portion of the subduction
interface. At the shallow portion beneath the accretionary prism, they
test several different frictional properties, including
velocity-weakening, neutrally stable, and velocity-strengthening. Their
preferred model that is validated against seafloor deformation and GPS
data shows that the shallow portion of the subduction is velocity
strengthening. They find that waves radiated from deep slip reflect off
the seafloor, causing large stress changes to the shallow portion of the
subduction interface that drive the rupture through the
velocity-strengthening region to the trench. Lotto et al. (2017) perform
a series of numerical simulations that couple dynamic rupture and
tsunami propagation in 2D models to explore compliant prisms’ effects on
tsunamigenesis. A compliant prism with reduced rigidity is embedded in
an otherwise homogeneous material and the shallow portion of the
subduction plane beneath the prism has variable frictional properties
with velocity-weakening properties at its down-dip extension. They find
that increasing prism compliance enhances shallow slip, and that a more
velocity-weakening behavior leads to increased slip both beneath the
prism and further downdip along the plate boundary fault. Although these
studies include both depth-varying friction and rigidity in their
dynamic rupture models and shed some lights onto effects of the two
factors, they have their specific objectives other than comparing roles
of the two factors in rupture dynamics. With the recent series of
studies (Sallares and Ranero, 2019; Sallares et al., 2021; Prada et al.,
2021) emphasizing a dominant role of depth-varying upper-plate rigidity,
it is imperative for the scientific community to better understand roles
of the two factors in rupture characteristics of subduction zone
earthquakes.
In this study, we compare and quantify the roles of depth-varying
frictional properties and rigidity in depth-dependent rupture
characteristics of subduction zone earthquakes, using physics-based
dynamic rupture models. We include both depth-varying upper-plate
rigidity of Sallarès and Ranero (2019) and depth-varying frictional
properties along the subduction interface (e.g., Lay, 2015; Lay et al.,
2012; Scholz, 1998) in the target model. We also perform dynamic rupture
simulations on other comparative models. By contrasting rupture
characteristics from these models, we quantify roles of the two factors
in determining depth-dependent rupture characteristics observed in
recent large subduction zone earthquakes.
2 Methods and Model
2.1 Dynamic rupture simulator
We use the three-dimensional finite element code EQdyna (Duan &
Oglesby, 2006; Duan, 2010; Duan, 2012; Luo & Duan, 2018; Liu & Duan,
2018) to simulate a suite of dynamic rupture scenarios. EQdyna is an
explicit finite element (FEM) dynamic rupture simulator that has been
verified in the Southern California Earthquake Center/U.S. Geological
Survey (SCEC/USGS) Spontaneous Rupture Code Verification Project (Harris
et al., 2009; Harris et al., 2011; Harris et al., 2018). In this
research, seismic waves propagate in an elastic medium and rupture on
the fault is governed by the rate-and-state friction law (RSF) with
aging law (Dieterich, 1979) implemented in EQdyna (Luo & Duan, 2018) to
explore major features of earthquake ruptures in the dynamic phase,
following equation (1):
\begin{equation}
\tau=\sigma(f_{0}+a\ln\frac{V}{V_{0}}+b\ln\frac{V_{0}\theta}{D_{c}})\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1),\nonumber \\
\end{equation}where a and b are constitutive frictional parameters
determined in laboratory experiments, Dc is the
critical slip distance for the exponential healing process after a
velocity stepping, and f0 (set to be 0.6) is a
reference friction coefficient associated with a reference steady state
slip rate V0 (set to be 10-6m/s). The state variable, θ , is a description of sliding history
and evolves as a function of V , θ , andDc according to the aging law:
\begin{equation}
\frac{\text{dθ}}{\text{dt}}=1-\frac{\text{Vθ}}{D_{c}}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2).\nonumber \\
\end{equation}2.2 Fault geometry and boundary conditions
Our model geometry has length of 190 km along strike, width of 200 km
perpendicular to strike, and depth of 50 km. The fault plane dips
(\(\phi\)) 15˚ and has both length and width of 150 km (Figure 1). The
top boundary of the model (Z = 0) intersects with the free surface, and
the side and bottom boundaries are perfectly matched layer (PML) that
absorbs seismic waves (Liu & Duan, 2018). The left (X =
Xmin = \(-\) 95 km) and right (X = Xmax= 95 km) boundaries are fixed along X-axis (i.e., zero displacement). We
create the finite element (FE) mesh of the model largely using
hexahedral elements for computational efficiency, with
fault-node-spacing of 200 m. To conform the shallow-dipping
(\(\phi=15\)) fault geometry, we cut a hexahedral element into two
wedge elements along the fault plane based on the degeneration technique
(e.g., Duan, 2010; Duan, 2012; Hughes, 2000; Luo & Duan, 2018). The
element sizes around the fault along the x-axis, y-axis, and z-axis are\(\Delta x\) = 200 m, \(\Delta y=\Delta x\cos\phi\) = 193 m, and\(\Delta z=\Delta x\sin\phi\) = 52 m, respectively.