Figure 3. Fault friction and hanging wall velocity structure setup in each scenario. Left panels: friction setup. Middle panels: hanging wall velocity setup. Right panels: footwall velocity setup. Color coding in the middle and right panels: the red curve indicatingVp , blue curve indicatingVs , and black dashed curve indicating density.
2.4 Rupture nucleation and resolution
In the dynamic rupture simulation, a region on the fault with velocity-weakening property is a necessary condition for nucleation. To initiate an instability, this velocity-weakening zone must be larger than a critical nucleation patch size h* , which is determined by the energy balance of a quasi-statically expanding crack (Lapusta et al., 2000; Rice, 1993; Rubin & Ampuero, 2005). Here, we use an estimation for 3D modeling according to Chen and Lapusta (2009) and Lapusta and Liu (2009):
\begin{equation} h^{*}=\frac{\pi}{2}\frac{\mu^{*}bD_{c}}{(b-a)^{2}\sigma}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3),\nonumber \\ \end{equation}
where \(\mu^{*}\) is \(\mu\) for a mode III crack and\(\mu/\left(1-\nu\right)\) for a mode II crack, \(\mu\) is the shear modulus, and \(\nu\) is the Poisson’s ratio. In our simulations, we assign a nucleation patch at 110 km along dip with radius of 4 km and slip rate of 0.01 m/s to artificially initiate a rupture event.
During dynamic rupture process, shear stress and slip rate change dramatically in the cohesive zone at the rupture front, which requires a certain number of elements to resolve these features (Day et al., 2005). The spatial resolution of the cohesive zone is thus critical for simulating dynamic rupture propagation (Day et al., 2005), which constrains the element size of the model (e.g., Duan & Day, 2008). The size of the cohesive zone, Λ0, at rupture speed\(v_{R}\) = 0+ under the RSF law follows
\begin{equation} \Lambda_{0}=C_{1}\frac{\mu^{*}L}{\text{bσ}}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (4),\nonumber \\ \end{equation}
where C1 is a constant of 9\(\pi\)/32 (Lapusta & Liu, 2009). For our FEM scheme, it is found that\(\Lambda_{0}/\Delta x\) of 2.4 with an element size \(\Delta x\) of 200 m can well resolve the cohesive zone (Meng et al., 2022). Taken the parameters choices in Section 2.3 , we set the model parameters considering equations (3) and (4).
Another consideration for resolution is time step. For dynamic rupture and seismic wave propagation, the time step (dt ) isα d/Vp , where α is a constant between 0 and 1 and d is the minimum element size (e.g., Liu et al., 2021). Given d = Δz = 52m, Vp = 6.7 km/ s, and α = 0.26, we set dt = 0.002 s.
2.5 Rupture scenarios
Figure 3 shows the fault friction and hanging wall material property setup for five dynamic rupture scenarios. The dynamic rupture scenarios are all nucleated at 110 km. All scenarios incorporate a downdip transition from velocity-weakening behavior at 120-km downdip (depth of ~30 km) to velocity strengthening behavior at 150-km downdip (depth of ~40 km) (white area between the fault patch and the dashed line in Figure 1). To account for the effects of updip transition from velocity-strengthening behavior near the trench to velocity-weakening behavior downdip, we set the along-dip transition distance L , together with a conditionally stable layer with an along-dip distance d (Figure 2a; left panels in Figure 3). For the depth-varying rigidity, we incorporate a multi-layered non-uniform velocity structure of the upper plate following Sallares & Ranero (2019) (middle and right panels in Figure 3). Scenario 1 is a reference scenario that assumes homogeneous velocity structure (Figures 3b and 3c) and homogeneous friction with L of 0 and d of 0 (Figure 3a). Scenario 2, on the other hand, is a most realistic setup that includes depth-varying velocity structure in the hanging wall (Figure 3e) and a two-layer velocity structure footwall (Figure 3f) and depth-varying friction with L of 40 km and d of 30 km (Figure 3d) among the scenarios. Scenarios 3 and 4 aim to quantify the effects of friction, both with homogeneous velocity structure in the both walls (Figures 3h, 3i, 3k, 3l) and depth varying friction withL of 40 km, but different d values of 0 km and 30 km (Figures 3g and 3j), respectively. Finally, Scenario 5 aims to quantify the effects of depth-varying rigidity, with depth-varying velocity structure in the hanging wall (Figure 3n) and a two-layer velocity structure footwall (Figure 3o) and homogeneous friction (both Land d equal 0) (Figure 3m).
3 Results
We present the simulation results of stress drop, total slip, rupture time contours, and rupture velocity on the fault plane in each Scenario. We also analyze the frequency contents of slip rate at selected on-fault stations to examine seismic radiation. We first compare the most realistic model (Scenario 2) and the reference model (Scenario 1) to examine combined effects of the two factors, namely fault friction and upper-plate rigidity. Then we examine other models and compare them with Scenario 1 and/or Scenario 2 to determine roles of the two factors in depth-dependent rupture characteristics.
3.1 Combined effects of depth-varying friction and rigidity on rupture characteristics
We examine the results of Scenarios 1 and 2 for the combined effects of heterogeneous friction and rigidity (Figure 4). The stress change distribution and slip distribution are significantly different between Scenario 1 (homogeneous rigidity and homogeneous friction) and Scenario 2 (heterogeneous rigidity and heterogeneous friction with L of 40 km and d of 30 km). The stress change distribution on the fault patch in Scenario 1, except for the boundary surrounded by velocity strengthening area (Figure 1), is negative (i.e., stress drop) from downdip toward the trench (Figure 4a). Correspondingly, slip reaches the trench in Scenario 1 with a maximum slip of ~20 m occurring at the trench. On the other hand, in Scenario 2, the updip transition from velocity-weakening behavior downdip to velocity-strengthening behavior updip (L of 40 km) diminishes slip at shallow depth, though free surface effects cause some obvious slip at the trench (Figure 4g). The maximum slip of 8m occurs at depth in this scenario. Correspondingly, stress drop (blue) mainly occurs at depth, while stress increases (red) at shallow depth (Figure 4f). Because both heterogeneous fault friction and heterogeneous wall rock properties are included in Scenario 2, we will examine contributions from each of the two factors to the above features in the slip and stress change distributions by other comparative scenarios in the sections below.
Rupture times are direct outputs from our dynamic rupture models. At the rupture time, a fault node reaches a slip rate of 0.01 m/s as the first time during the simulation. Rupture velocities are calculated from rupture times, following a method proposed by (Bizzarri and Das, 2012). In Scenario 1, the rupture time contour (Figure 4c) and the rupture velocity distribution (Figure 4d and 4e) both indicate that the rupture generally accelerates toward the trench from the nucleation patch on a subduction plane with a uniform velocity-weakening friction property embedded in a uniform medium. In Scenario 2, the rupture accelerates within the velocity-weakening patch from the nucleation patch but slows down when it propagates into the conditionally stable part (d of 30 km) and the updip transition patch (L of 40km) (Figure 4h, 4i, and 4j), in particular along the depth profile at along-strike-distance of 50km (black curve in Figure 4j), except near the trench. Supershear rupture occurs near the trench in both scenarios due to effects of free surface and shallow-dipping fault geometry. The slow rupture velocity at shallow depth in Scenario 2 may be attributed to combined effects of updip transition in friction and low-Vp and -Vs at shallow depth. The other comparative scenarios will help clarify and quantify their roles in slower rupture velocity (and thus longer duration) at shallow depth in Scenario 2 than that in Scenario 1.