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.