# The Influence of Depth-Varying Elastic Properties of the Upper Plate on Megathrust Earthquake Rupture Dynamics and Tsunamigenesis

Manel Prada, Percy Galvez, Jean Paul Ampuero, Valentí Sallarès, Carlos Sánchez-linares, Jorge Macías, Daniel Peter

► To cite this version:

Manel Prada, Percy Galvez, Jean Paul Ampuero, Valentí Sallarès, Carlos Sánchez-linares, et al. The Influence of Depth-Varying Elastic Properties of the Upper Plate on Megathrust Earthquake Rupture Dynamics and Tsunamigenesis. Journal of Geophysical Research: Solid Earth, 2021, 126 (11), pp.e2021JB022328. 10.1029/2021JB022328. hal-03462491

HAL Id: hal-03462491

https://hal.science/hal-03462491

Submitted on 6 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2021JB022328

**Key Points:**

* We use 3D simulations to evaluate the influence of realistic upper-plate rigidity on megathrust earthquake properties and tsunamigenesis

* Simulations show that realistic upper-plate rigidity variations explain the depth-dependent behavior of earthquake dynamic properties

* Overestimation of upper-plate rigidity leads to underestimation of co-seismic seafloor deformation and tsunami amplitude in our simulations

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

M. Prada, mprada@icm.csic.es

**Citation:**

Prada, M., Galvez, P., Ampuero, J.-P., Sallarès, V., Sánchez-Linares, C., Macías, J., & Peter, D. (2021). The influence of depth-varying elastic properties of the upper plate on megathrust earthquake rupture dynamics and tsunamigenesis. *Journal of Geophysical Research: Solid Earth*, 126, e2021JB022328. https://doi.org/10.1029/2021JB022328

Received 30 APR 2021
Accepted 29 OCT 2021

**Author Contributions:**

**Conceptualization:** Manel Prada, Percy Galvez, Jean-Paul Ampuero, Valenti Sallarès

**Funding acquisition:** Manel Prada

**Methodology:** Manel Prada, Percy Galvez, Jean-Paul Ampuero, Carlos Sánchez-Linares, Jorge Macías

**Project Administration:** Manel Prada

© 2021 The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution-NonCommercial License,
which permits use, distribution and
reproduction in any medium, provided
the original work is properly cited and
is not used for commercial purposes.

<figure><img src="image_2.png" /></figure>

# The Influence of Depth-Varying Elastic Properties of the Upper Plate on Megathrust Earthquake Rupture Dynamics and Tsunamigenesis

Manel Prada<fnref n="1" /> ID, Percy Galvez², Jean-Paul Ampuero<fnref n="3" /> ID, Valenti Sallarès<fnref n="1" /> ID, Carlos Sánchez-Linares⁴, Jorge Macías⁴ ID, and Daniel Peter² ID

**Abstract** Megathrust earthquakes are strongly influenced by the elastic properties of rocks surrounding the fault. In contrast to friction, these properties can be derived in situ from geophysical measurements along the seismogenic zone. However, they are often overestimated in numerical simulations, particularly in the shallow megathrust. Here we explore the influence that realistic depth-varying upper-plate elastic properties along the megathrust have on earthquake rupture dynamics and tsunamigenesis using 3D dynamic rupture and tsunami simulations. We compare results from three subduction zone scenarios with homogeneous and heterogeneous elastic media, and bimaterial fault. Elastic properties in the heterogeneous model follow a realistic depth-distribution derived from controlled-source tomography models of subduction zones. We assume the same friction properties for all scenarios. Simulations in the heterogeneous and homogeneous models show that rigidity depth variations explain the depth-varying behavior of slip, slip rate, frequency content, and rupture time. The depth-varying behavior of slip, frequency content, and rupture duration quantitatively agree with previous predictions based on worldwide data compilations, explaining the main depth-dependent traits of tsunami earthquakes and large shallow megathrust earthquakes. Large slip, slow rupture and slip rate amplification in bimaterial simulations are largely controlled by elastic rock properties of the most compliant side of the fault, which in subduction zones is the upper plate. Large shallow slip and trenchward increasing upper-plate compliance of the heterogeneous model lead to the largest co-seismic seafloor deformation and tsunami amplitude. This highlights the importance of considering realistic upper-plate rigidity variations to properly assess the tsunamigenic potential of megathrust earthquakes.

**Plain Language Summary** The largest earthquakes on Earth occur along the megathrust fault in subduction zones. When these events reach the shallow portion of the fault (<5 km depth), they are slower, generate larger slip, and release lower frequency content than in the deep region of the fault (i.e., 25 km depth). The presence of large slip at shallow depths is responsible for the large seafloor displacement that triggers the generation of devastating tsunamis. Thus, it is important to understand which factors determine the depth-dependent behavior of these rupture properties to properly assess the tsunamigenic potential of megathrust earthquakes. In this work, we use 3D numerical dynamic rupture and tsunami simulations to demonstrate that this particular depth-varying behavior is highly influenced by realistic variations of elastic rock properties overlying the fault in the upper plate. The overestimation of these properties in the shallow portion of the megathrust lead to underestimation of slip and seafloor displacement, and thus tsunami hazard.

## 1. Introduction

Megathrust earthquakes nucleate within the seismogenic zone (Byrne et al., 1988; Hyndman et al., 1997). The updip limit of this region may vary depending on material properties and thermal conditions along the megathrust (Hyndman et al., 1997), but it is commonly defined between 5 and 10 km depth (Scholz, 1998). Yet, megathrust earthquakes may occasionally rupture through the shallow and apparently aseismic region of the fault (<5 km of depth), particularly in areas with sediment-starved trenches and irregular subducting

<footnote n="1">¹Barcelona Center for Subsurface Imaging, Institut de Ciències del Mar (ICM), CSIC, Barcelona, Spain, ²Physical Science and Engineering Division, King Abdullah University of Science and Technology, Thuwal, Saudi Arabia,</footnote>

<footnote n="3">³Institut de Recherche pour le Développement, Université Côte d'Azur, Nice, France, ⁴EDANYA, Universidad de Malaga, Málaga, Spain</footnote>
**Resources:** Percy Galvez, Carlos Sánchez-Linares, Jorge Macías, Daniel Peter

**Software:** Manel Prada, Percy Galvez, Carlos Sánchez-Linares, Jorge Macías, Daniel Peter

**Supervision:** Jean-Paul Ampuero, Valenti Sallarès

**Writing – original draft:** Manel Prada

**Writing – review & editing:** Percy Galvez, Jean-Paul Ampuero, Valenti Sallarès, Carlos Sánchez-Linares, Jorge Macías, Daniel Peter

topography (e.g., Geersen, 2019; Polet & Kanamori, 2000), and heterogeneous sediment thickness and elevated pore pressure (Li et al., 2018; Tobin & Saffer, 2009). When this occurs, the rupture produces anomalously large slip near the trench, dramatic seafloor deformation and large and devastating tsunamis. This is the case of tsunami earthquakes, which are defined as a particular type of shallow events that excite disproportionally large tsunamis for their moderate seismic magnitude (Kanamori, 1972). In addition, these are relatively slow earthquakes that radiate seismic waves depleted in high frequencies, which in turn leads to a large discrepancy between their surface-wave magnitude (Ms) and moment magnitude (Mw) (e.g., Kanamori & Kikuchi, 1993; Newman et al., 2011). However, an increasing number of observations demonstrate that systematic differences of source properties are not only particular of tsunami earthquakes, but all shallow megathrust events display a depth-dependent behavior of earthquake rupture characteristics (Lay & Bilek, 2007; Lay et al., 2012). Large tsunamigenic earthquakes such as the 2004 Sumatra-Andaman (Mw9.2), 2010 Maule (Mw 8.8), and 2011 Tohoku-Oki (Mw 9.0) produced larger slip near the trench than at deeper portions of the fault (e.g., Lay et al., 2012) and radiated low-frequency energy mostly from the shallow portion of the megathrust (e.g., Kanamori & Yomogida, 2011; Koper et al., 2012; Lay et al., 2012; Meng et al., 2011; Simons et al., 2011).

These seismological observations indicate that such depth-dependent earthquake behavior shares similarities between contrasting tectonic environments (erosional and accretionary margins) and earthquake sizes, suggesting the existence of a common contributing factor. The presence of low-rigidity subducting sediments has been invoked to explain some characteristics of shallow ruptures, in particular the large shallow slip and slow rupture (Bilek & Lay, 1999; Polet & Kanamori, 2000). However, numerical simulations show that the presence of a low velocity fault zone may result in acceleration of the rupture to super-shear velocity and produce high-frequency radiation (e.g., Huang & Ampuero, 2011; Huang et al., 2016), in contrast with seismological observations of shallow events (e.g., Kanamori & Kikuchi, 1993). On the other hand, results from numerical simulations explain the large slip and slow rupture observed in tsunami earthquakes with large compliant sedimentary prisms and velocity-weakening fault friction conditions (Lotto et al., 2017). However, earlier numerical results using depth-dependent material properties and rate-and-state friction show that velocity-strengthening materials can also develop large near-trench slip (Kozdon & Dunham, 2013). In addition, some tsunami earthquakes, like the 1992 Nicaragua event, occurred in erosional convergent margins with no accretionary prisms (McIntosh et al., 2007; Sallarès et al., 2013). This implies that the cause of large shallow slip, tsunamigenesis and slow rupture is possibly controlled by factors other than well-developed sedimentary wedges or specific friction laws. More recent research proposes that, rather than fault plane properties, the depth-varying behavior of earthquake properties, including slip, rupture velocity and frequency content, may be all explained by the distribution of elastic properties of rocks overlying the fault in the upper plate (Sallarès & Ranero, 2019; Sallarès et al., 2021). Overall, all these studies share a common ground, in which the rupture properties are influenced by elastic rock properties surrounding the fault.

Over the last decade, there has been an increasing number of numerical studies attempting to understand the dynamics of megathrust earthquakes and their depth-varying behavior. Some studies invoked inelastic wedge deformation to account for the depth-dependent behavior of tsunami earthquakes (Ma, 2012; Ma & Hirakawa, 2013). Alternatively, depth-dependent initial stress conditions and/or friction laws are invoked to explain the depth-varying behavior of large earthquakes such as the 2011 Tohoku-Oki event (e.g., Huang et al., 2012; Kozdon & Dunham, 2013; Murphy et al., 2018). Fault friction is an important aspect that, in addition to elastic properties, influences some depth-dependent characteristics of megathrust earthquakes as it controls slip stability, fracture energy and stress drop. Yet, the distribution of frictional parameters along faults is highly unknown, as it depends on site-specific factors such as lithology, temperature and/or water content (Scholz, 2019), which are all difficult to quantify in situ along the entire seismogenic zone. In contrast, elastic properties can be quantified from geophysical measurements, especially from marine seismology (e.g., Sallarès et al., 2021). Yet, the distribution of elastic properties considered in previous numerical simulations are homogeneous (e.g., Huang et al., 2012), layered (e.g., Galvez et al., 2016), and/or constant for a given geological unit (e.g., Kozdon & Dunham, 2013; Moreno et al., 2012; Murphy et al., 2018; van Zelst et al., 2019) and, in most cases, it is not extracted from data. This hinders understanding of the additional role that realistic elastic properties along the fault have on the dynamics of megathrust events, and may lead to dynamic models in which depth-dependent fault friction is the only ingredient available to explain
seismological observations. On the other hand, in some models (e.g., Kozdon & Dunham, 2013; Ma, 2012), the value of rigidity considered in the shallow megathrust is one order of magnitude larger than the value estimated from drilling samples (Jepsson et al., 2018). The overestimation of elastic properties involved in the shallow rupture may lead to underestimation of slip, uplift and tsunami size (e.g., Geist & Bilek, 2001; Satake, 1994), thus underestimating related hazards.

Here, we use 3D numerical simulations to assess differences in rupture behavior between the deep and shallow megathrust attributed to depth-dependent material properties. We show that incorporating a realistic distribution of elastic properties has a first-order effect on the depth-dependent behavior of megathrust earthquakes and tsunamigenesis. We conduct 3D dynamic rupture simulations of subduction earthquakes with and without depth-dependent elastic properties, and we compare their outcomes in terms of slip, rupture duration, rupture velocity, and frequency content. To isolate differences in rupture properties attributed to differences in elastic properties, all simulated scenarios share the same friction properties and fault geometry. We find that depth-dependent elastic properties alone explain first-order observations targeted here. We then explore the influence of bimaterial fault properties, in models that include the material contrast between the overriding and downgoing plates. Finally, we assess the influence of depth-dependent rock rigidity on seafloor deformation and tsunamigenesis by carrying tsunami simulations.

## 2. Dynamic Rupture Model Setup

We consider a megathrust fault with a dip angle of 15° based on controlled-source tomographic models from the Middle American subduction zone (Sallarès et al., 2013, 2021). We model three subduction scenarios: a homogeneous elastic medium, a heterogeneous medium, and a bimaterial fault. In the homogeneous scenario, the values of P-wave velocity ($V_p$), S-wave velocity ($V_s$), and density ($\rho$) are those of rocks of the overlying the megathrust at 20–25 km depth (Figures 1d–1g). By setting the reference value at these depths, we can assess the relative behavior of rupture properties between the deep and shallow portion of the megathrust when comparing results between the homogeneous and heterogeneous model. In addition, given that this reference is the same than that used in Sallarès and Ranero (2019), this choice allows for a straight comparison between results. The heterogeneous scenario includes the generic depth-distribution of elastic properties of the upper plate inferred from worldwide controlled-source tomographic models (Sallarès & Ranero, 2019) (Figures 1d–1g). Both the upper plate and downgoing plate are characterized by the same 1D depth-dependent relation (Figure 1h). Depth variations in the heterogeneous model are not only due to changes across geological units (e.g., crystalline crust vs. sediments), but may also reflect variations in the tectonic framework within the same geological unit (e.g., density of fractures). This simplified representation of the heterogeneous distribution of elastic properties of the upper plate is consistent with first-order variations in which elastic moduli decrease toward the upper-plate toe (Figure 1h) (Calahorrano et al., 2008; Contreras-Reyes et al., 2017; Sallarès et al., 2013). This distribution of elastic properties is, therefore, ideal to explore how realistic elastic properties influence the rupture and the permanent upper-plate co-seismic deformation.

We divide the fault in three different domains based on variations in the tectonic structure of the upper plate observed from seismic data and on the rates of variations of elastic properties overlying the megathrust fault (Sallarès & Ranero, 2019). The “regular domain” extends from 30 to 10 km depth, the “transitional domain” from 10 to 5 km depth, and the “shallow domain” from 5 to 0 km depth (Figures 1d–1g). To capture the rupture properties characteristic of each domain, we performed three different sets of simulations, each confined to one of the three domains. We performed two simulations in each domain, one in the homogeneous medium and a second one in the heterogeneous medium. In addition, we ran simulations of ruptures breaking through the three domains simultaneously to assess rupture velocity and slip rate differences.

For all simulations, we assume a linear slip-weakening friction law, with constant critical slip-weakening distance (Dc) of 0.4 m, static friction coefficient ($\mu_s$) of 0.6 and dynamic friction coefficient ($\mu_d$) of 0.4. The initial normal ($\sigma_n$) and shear stresses ($\tau_0$) are 100 and 50 MPa, respectively, through the regular domain, and decrease linearly across the transitional and shallow domains to 10 and 5 MPa at the surface, respectively (Figures 2a–2c). Indeed, absolute values of rupture properties such as slip or slip rate vary depending on the choice of slip-weakening frictional parameters (i.e., $\mu_s$, $\mu_d$, Dc), but since we are comparing simulations
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Unstructured mesh of (a) the whole domain and focused view of (b) the fault domain and (c) fault surface. The rupture area is represented by the shaded area. (d–g) Depth distributions of P-wave velocity ($V_p$), S-wave velocity ($V_s$), density ($\rho$), and rigidity ($\mu$), respectively, in the homogeneous (blue) and heterogeneous (red) models. Orange, yellow and white background colors show the extent of the shallow, transitional and regular domains, respectively. (h) The spatial distribution of $V_p$ in the heterogeneous model domain.</figcaption></figure>

with the same frictional properties the relative behavior between the deep and shallow portions of the fault remains, essentially reflecting variations in material properties between models. Similarly, considering a different friction law (i.e., rate-and-state friction) would not significantly change the outcome of such comparison. We run simulations in which the rupture is confined to a particular domain of the fault by setting
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Frictional strength and initial stresses for the shallow (orange background), transitional (yellow background), and regular (white background) domains on a cross section along the center of the fault. The blue curve is the dynamic strength, black is initial shear stress, red is the static shear strength, and gray is the initial normal stress. The three areas of elevated initial shear stress (black) are the nucleation zones of each of the three simulation sets.</figcaption></figure>

$\mu_s$ to $10^4$ outside of it, which creates an artificial fault barrier that prevents the rupture from propagating further. The nucleation zone in each domain was set as a rectangular region where initial shear stress exceeds the static shear strength, and its size was estimated to be large enough to sustain rupture propagation (Day et al., 2005; Galis et al., 2015) (Figures 2a–2c). The nucleation patches are located at a depth of 24, 9 and 4.5 km for the simulations through the regular, transitional, and shallow domains, respectively (Figure 2).

To solve the dynamic rupture problem we use the open-source spectral element code SPECFEM3D (Peter et al., 2011) and its dynamic rupture capabilities (Galvez et al., 2014; Kaneko et al., 2008). The simulations use an unstructured mesh generated with the software TRELIS (Figure 1). The mesh is 240 km wide, 270 km long and 80 km deep (Figure 1a). The fault domain is a confined volume at the center of the mesh that extends over the first 30 km of depth, and is 240 km wide and 270 km long (Figure 1b). Our models rupture a square patch in the center of the fault that is 120 km wide and ~110 km long (Figure 1c). Clayton-Engquist absorbing boundary conditions were set for all simulations (Clayton & Engquist, 1977; Komatitsch & Tromp, 1999). The spectral element size in the fault domain is 500 m with 5 Gauss-Lobatto-Legendre nodes per element edge, yielding an average grid size of 125 m. This grid size is smaller than the cohesive zone, which decreases in the heterogeneous model to ~130 m in the shallowmost portion of the fault as a result of the decreasing rigidity. Outside the fault domain, the element size gradually increases to 1 km at the boundaries of the mesh.

## 3. Numerical Results

We simulated ruptures confined to each of the three fault domains, as well as whole-fault ruptures. The results in the heterogeneous medium show clear depth-varying values of slip, source duration, frequency content, rupture velocity, and slip rate, in contrast to simulations in the homogeneous medium. We report relative values between simulation results, as absolute values are of secondary importance in this study, given that our main purpose is to assess the relative rupture behavior between the deep and shallow portions of the fault.

### 3.1. Ruptures Confined to Each Fault Domain

#### 3.1.1. Slip

Figure 3 shows the spatial distribution of slip ratio between the heterogeneous and homogeneous models (Slip$_{het}$/Slip$_{hom}$) through each of the three domains (absolute values are shown in Figure S1 in Supporting Information S1). In all fault domains, the rupture through the heterogeneous model results in larger slip than in the homogeneous (Figure 3); particularly in the shallowmost 5 km of the fault (i.e., shallow domain), where the slip ratio increases from 3 to 9 toward the trench (Figure 3; max Slip$_{het}$ ~ 20 m; max Slip$_{hom}$ ~ 3 m in Figure S1 in Supporting Information S1). Despite these slip differences, similar moment magnitudes are obtained in both elastic models in the shallow ($M_0$ ~ $2 \times 10^{20}$ Nm; Mw ~ 7.5), transitional ($M_0$ ~ $2 \times 10^{20}$ Nm; Mw ~ 7.5), and regular domains ($M_0$ ~ $5 \times 10^{21}$ Nm; Mw ~ 8.4) as we are using the same rupture area (S) and stress drop ($\Delta\tau$) in all simulations, and

$$M_o = \Delta\tau S^{3/2}C^{-1}, \quad (1)$$

where C is a geometric constant of order one.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Ratio of final slip between heterogeneous and homogeneous models (Slip<sub>het</sub>/Slip<sub>hom</sub>) in the three pairs of simulations, each confined to one fault domain (absolute slip values are shown in Figure S1 in Supporting Information S1). Red stars are the nucleation points of each confined rupture. Dashed black lines depict boundaries between fault domains.</figcaption></figure>

### 3.1.2. Source Duration and Frequency Content

Simulations through the heterogeneous and homogeneous media show increasing differences in source durations and moment rate spectra toward the trench (Figures 4d–4f, S2 and S3 in Supporting Information S1). Source duration is taken at zero crossing of the moment rate functions in Figure 4, and its ratio (Duration<sub>het</sub>/Duration<sub>hom</sub>) increases from 1.3 in the transitional domain to 1.6 in the shallow domain (Figures 4a–4f and S3 in Supporting Information S1). These variations translate into a trenchward reduction of corner frequency ($f_c$), which is estimated as the inverse of the duration of each moment rate function following (Udías et al., 2014), implying an increasing depletion of high-frequency content toward the trench (Figure 4). The largest $f_c$ differences occur in the shallow domain, where $f_c$ in the heterogeneous model is twice lower than in the homogeneous model. Minor $f_c$ differences are observed in the transitional domain, while no differences are observed in the regular domain (Figures 4d–4f).

## 3.2. Whole-Fault Rupture

To assess rupture velocity and slip rate differences attributed to depth-dependent elastic properties we simulated a continuous rupture throughout the three domains (Figure 5), and compared the rupture times and slip rate values obtained in both homogeneous and heterogeneous media along a 2D cross section through the center of the fault (Figures 6d and 7). The nucleation was set at 24 km depth and friction was modified in the shallowmost 3 km of the fault to prevent supershear rupture propagation. The main reason to suppress such effect, is that, to our knowledge, shallow supershear ruptures are not reported in megathrust earthquakes. In addition, shallow supershear ruptures promote high-frequency radiation (e.g., Huang & Ampuero, 2011), in contrast to seismological observations of shallow megathrust earthquakes. This way, we increased $\mu_d$ linearly from 0.4 at 10 km depth to 0.8 at the surface (Figure 5a), resulting in a region with negative stress drop ($\tau_0 < \tau_d$) in the shallowmost 3 km and slip-strengthening ($\mu_d > \mu_s$) in the topmost kilometer (Figure 5a). Even though this results in reduction of the final slip near the trench (Figure S4 in Supporting Information S1) as shown by previous numerical studies (Duan, 2012; Galvez et al., 2016; Huang et al., 2014), the rupture through the heterogeneous model consistently produces larger slip, particularly in the shallow domain (Figures 5b and S4 in Supporting Information S1). Both ruptures have similar seismic moment (Mw ~ 8.5).

According to theory, this modified frictional setup should prevent supershear propagation (Andrews, 1976). However, supershear rupture prevails at shallow depth in the homogeneous model, in which a small daughter crack rupturing at supershear velocity appears ahead of the main rupture front in the shallowmost 6–5 km (Figure 6c).

### 3.2.1. Rupture Time and Velocity

To compare rupture time differences without the effect of the daughter crack from the homogeneous model in the comparison, we consider rupture times only for the main rupture front (Figure 6c). In addition, we spatially smoothed the rupture times of the heterogeneous model to suppress the effect of steps generated by the interaction of the rupture with the negative stress drop and slip-strengthening friction at shallow depth (Figure 6d) and capture only first-order variations.

Rupture time differences increase trenchwards (dashed colored lines in Figure 6d). The main rupture front arrives at the trench at 29 s in the homogeneous model and at ~40 s in the heterogeneous model (Figure 6d), while the total source duration is ~38 s for the homogeneous and ~43 s for the heterogeneous model (Figure 6a). Rupture velocity ($V_r$), derived from rupture time curves, ranges between 70% and 90% of $V_s$ in both models (Figures 6e and 6f), consistent with empirical observations (e.g., Bilek & Lay, 1999).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Moment rate functions of ruptures in the (a) shallow, (b) transitional, and (c) regular domains, and their corresponding moment rate spectra (d–f, respectively). Blue and red curves correspond to the results from the homogeneous and heterogeneous models, respectively. The vertical lines in the spectra are the estimated corner frequencies, $f_c$.</figcaption></figure>

Depth- variations of $V_s$ in the heterogeneous model induce similar changes in $V_r$ (Figure 6e). In both cases, rupture slows down below 70% of $V_s$ in the shallowest ~5 km of the fault, because of the increasing $\mu_d$ as- sumed in our simulations (Figure 5a).

### 3.2.2. Slip Rate

We explore here slip rate differences between the two models that rupture through the entire fault. Peak slip rate values along a 2D cross section in the center of the fault show increasing differences between elastic models toward the trench (Figure 7a), where slip velocity is up to 5 times larger in the heterogeneous model. Interestingly, when the rupture reaches the surface, it reflects back resulting in the acceleration of slip rate near the trench (Movies S1a and S1b). This process, which has been attributed to the interaction of the rupture with the free surface in previous numerical studies (e.g., Huang et al., 2012; Ma & Beroza, 2008), occurs in both elastic models (Movies S1a and S1b), but its effect is larger in the heterogeneous one (Movie S1b).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Assumed depth distributions of initial stress and strength in the simulations of rupture through the entire fault depth. The blue line is the dynamic strength, black is the initial shear stress, red is the static strength, and gray is the initial normal stress. (b) Final slip ratio between the heterogeneous and homogenous models (Slip<sub>het</sub>/Slip<sub>hom</sub>). Red star is the nucleation point. Final slip values for each elastic model are shown in Figure S4 in Supporting Information S1.</figcaption></figure>

# 4. Discussion

## 4.1. Depth-Varying Dynamic Rupture Properties

The different dynamic rupture properties presented here show a clear depth dependence that arises from variations in elastic properties. The trenchward decrease of rock rigidity enhances slip in the shallow domain, where it is 3–10 times larger than in the regular domain, where rigidity is nearly one order of magnitude larger (5 vs. 40 GPa; Figures 1g and 8a). Peak slip rate is also enhanced as rigidity decreases trenchwards (Figure 7b), an effect that is somewhat expected, as peak slip rate is proportional to Vr/μ (Ohnaka et al., 1987). Increasing slip rate toward the trench in the heterogeneous model causes stronger ground motion near the trench than in the homogeneous model (Figures 7c–7e). The corner frequency fc is up to twice lower in the shallow domain than in the deeper region along the regular domain. This is due to the longer duration of the event in the shallow region, which in turn results from the lower Vr. As shown in Figures 6e and 6f, Vr in our simulations is largely controlled by the depth distribution of Vs, which can be up to ~2 times slower in the shallow domain than in the regular domain (Figures 6e and 6f). The slip-strengthening condition imposed in the whole-fault rupture simulation causes the rupture to slow down and stop at the corners of the fault in the shallow domain (Figure 5b), implying that fault friction has additional effects slowing down Vr, although at expenses of decreasing slip near the trench (Figure S4 in Supporting Information S1). Thus, while fault friction may explain some depth-dependent properties of megathrust earthquakes, it fails to reconcile them all. This aspect is also observed in homogeneous models of the Tohoku-Oki earthquake, which managed to partially explain seismological observations of the event with fault friction variations along the fault, but found inconsistencies in the shallow region, where the amount of high frequency radiated was higher than observed (Huang et al., 2014). These observations are evidence that the depth-dependent behavior of megathrust earthquakes may be equally sensitive to
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Moment rate, and (b) moment rate spectra of ruptures through the entire fault. Blue and red lines correspond to results from the homogeneous and heterogeneous models, respectively. (c) Slip rate snapshots taken across the center of the fault at 3, 12, and 20 s of the homogeneous model simulation. (d) Rupture time as a function of depth across the center of the fault. Solid blue and red curves are the rupture time profiles for the homogeneous and heterogeneous cases, respectively. The dashed blue curve is the rupture time profile of the main rupture front through the homogeneous model, derived from slip rate profiles as shown at 20 s in (c). The red dashed curve is a polynomial fit of the first-order variations of the rupture time profile along the heterogeneous model, which yields a smoothed rupture time profile and removes from the comparison the steps in rupture time attributed to the negative stress drop and slip-strengthening region. Rupture velocity ($V_r$) profiles for the (e) heterogeneous and (f) homogeneous models. Gray bands show the range of velocity values between 70% and 90% of $V_s$, and the solid black lines show $V_s$.</figcaption></figure>

friction properties and elastic rock properties, particularly in the shallow region, where slip-strengthening (or velocity-strengthening) material may be present along the fault (e.g., Kozdon & Dunham, 2013; Moore & Saffer, 2001). To quantify the additional role that friction plays, however, requires first the incorporation of realistic depth-dependent elastic properties, which, in contrast to friction, are robustly established from tomographic models along the entire seismogenic zone (Sallarès et al., 2021) and, for rocks overlying the fault, follow a global generic downdip distribution (Sallarès & Ranero, 2019).

We compare depth-trends of slip ratio, frequency content ratio ($f_{chet}/f_{chom}$), and normalized rupture time ratio per unit length (i.e., rupture slowness ratio in Figure 8b) with those presented in Sallarès and Ranero (2019).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Depth profiles of peak slip rate along the center of the fault, for simulations rupturing the entire fault in the homogeneous (blue) and heterogeneous (red) models. Gray thick line is the depth distribution of rigidity from Sallarès and Ranero (2019). (b) Peak slip rate values of the heterogeneous model as a function of rigidity. (c) Cross section showing the location of fault points and seismic stations at which slip rate (d) and ground motion (e) values are measured. The color scale along the fault depicts rigidity of the heterogeneous model.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison between depth trends of (a) slip ratio, (b) Rupture slowness ratio (relative rupture duration ratio in Sallarès & Ranero, 2019), and (c) $f_c$ ratio from Sallarès and Ranero (2019) (white circles) and this study (red). Gray line in (a) is rigidity ratio ($\mu_{het}/\mu_{hom}$). Rupture slowness ratio (Heterogeneous/ Homogeneous) is calculated from rupture time curves in Figure 6d.</figcaption></figure>

The comparison shows that depth-trends from our study quantitatively agree with previous estimates (Figure 8). The excellent agreement of the $f_c$ ratio estimates implies that our results are also consistent with the depth-varying Ms-Mw discrepancy presented by Sallarès and Ranero (2019), explaining thereby, the depth-dependent megathrust earthquakes source characteristics (Lay et al., 2012). We show that large slip, slow rupture and high-frequency depletion are strong indicators of shallow ruptures and thus, potential indicators for high tsunami hazard, in agreement with seismological observations of tsunami earthquakes and larger shallow megathrust earthquakes (e.g., Kanamori & Kikuchi, 1993; Lay & Bilek, 2007; Newman et al., 2011). These results are also consistent with tsunami early warning studies (Lomax &MICHELINI, 2009, 2011), which reveal that tsunamigenic megathrust earthquakes are related to anomalously long apparent rupture duration (>50 s) and weak P-wave high-frequency content. However, rather than compliant fault zone sediments (Bilek & Lay, 1999; Polet & Kanamori, 2000), our study strongly suggests that this behavior should be attributed to the decreasing rigidity of rocks overlying the shallow megathrust, where rigidity decreases by one order of magnitude in the top 5 km (from ~20 to ~2 GPa; Figure 1g).

## 4.2. Bimaterial Fault

Subduction zones feature a material contrast between the overriding and downgoing plates, which may affect interplate earthquake characteristics (e.g., Ampuero & Ben-Zion, 2008; Rubin & Ampuero, 2007; Shlomai et al., 2021). To explore this, we conducted additional simulations considering a bimaterial fault interface. We assume the same 1D distribution of elastic properties from Sallarès and Ranero (2019) for the upper plate but a homogeneous media for the downgoing plate (Figure 9a). While this distribution may overestimate rigidity of the downgoing plate in its shallowmost region, it allows us to isolate the effect of each plate in rupture properties by direct comparison with previous simulations through the heterogeneous and homogeneous models (Figure 9).

To assess slip variations, we simulated three rupture scenarios, each confined to one of the three fault depth ranges. The results show similar slip values to those obtained previously in the heterogeneous model
<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) 3D view of the $V_p$ distribution of the bimaterial fault model. (b) Depth distribution of slip ratio between the heterogeneous and bimaterial models (green), the bimaterial and homogeneous models (blue), and the heterogeneous and homogeneous models (red). (c) Rupture slowness ratio between the heterogeneous and bimaterial models (green), and the heterogeneous and homogeneous models (red). The red and green curve are calculated from the corresponding dashed rupture time curves in Figure S7a in Supporting Information S1. (d) Slip rate from the bimaterial (green) and the heterogeneous (red) models at fault points 1 to 4 (deep to shallow) indicated in Figure 7c.</figcaption></figure>

(Figures 9b and S5 in Supporting Information S1). Thus, slip is mainly controlled by the medium with lower rigidity, which here is the upper plate.

Additionally, we simulate a rupture through the entire fault to compare rupture velocity and slip rate between the heterogeneous and bimaterial cases (Figures 9c, 9d and S6 in Supporting Information S1). We calculate the ratio of rupture time differences per unit of fault length between the heterogeneous and bimaterial models (rupture slowness in Figure 9c). The rupture slowness through the bimaterial model resembles that in the heterogeneous model, with some differences along the shallow domain, where the rupture accelerates and propagates at supershear velocity in the shallowmost 2 km in the bimaterial simulation. Based on results from the whole-fault rupture through the homogeneous model, this effect is possibly attributed to the
faster material of the downgoing plate assumed in our simulations. Using a more compliant and depth-dependent elastic structure for the downgoing plate could mitigate this effect. Previous theoretical and computational studies of bimaterial rupture indicate that rupture velocity $V_r$ is related to bimaterial contrast and controlled to a greater extent by the lowest $V_s$ among the two materials (Rubin & Ampuero, 2007; Shlomai et al., 2021). While supershear propagation obscures this relation in the shallow domain, it seems to be the case for the regular and transitional domains, where the rupture slowness ratio between the heterogeneous and the bimaterial scenarios is ~1 (Figures 9c, S6 and S7 in Supporting Information S1).

The bimaterial model shows larger slip rate values than the heterogeneous model in the shallowmost region, where the bimaterial contrast is larger (Figure 9d). This amplification of slip rate, and consequently of ground motion, is consistent with previous bimaterial simulations in which the overlying side of the fault (i.e., overriding plate) is more compliant than the underlying one (Ma & Beroza, 2008). This mechanism deserves further consideration in numerical simulations as it may have tsunamigenic implications. Although our setup does not have topography, the enhanced slip rate and ground motion occurs at distances of less than 10–15 km from the trench, where the continental slope is commonly found and seafloor may dip 5-15° trenchwards (e.g., Harders et al., 2011). These topographic features, in combination with strong ground shaking, may trigger slope failure processes and contribute to tsunamigenesis (e.g., Tappin et al., 2014).

Overall, results from the bimaterial simulation show a clear depth-dependent behavior of slip, rupture velocity and slip rate that is largely controlled by the elastic properties of the softer and slower material around the fault. However, as mentioned above, elastic properties of the lower plate may not be realistic. While the distribution of the lower plate properties is often poorly resolved along the entire seismogenic zone, tomographic models show a more compliant structure of the lower plate in the shallowest region of the megathrust (e.g., Calahorrano et al., 2008; Contreras-Reyes et al., 2017), which could contribute to enhance slip and decrease rupture velocity. Yet, the limited information we have indicates that variations in seismic velocities occur spatially faster in the downgoing plate along the shallowmost region of the megathrust because of sediment compaction and over-pressured fluid release (Calahorrano et al., 2008). The different seismic structure results from the different tectonic processes acting on both plates. The upper plate is controlled by contractional structures that intensify trenchwards as a result of the convergence (Kodaira et al., 2017; von Huene et al., 1994), while the downgoing plate is overlaid by subducting sediments and mostly affected by extensional bending-related faulting (Ranero et al., 2003). These differences and the fact that normal stress is higher on the downgoing plate, make the upper plate more compliant than the downgoing oceanic plate in nature. Thus, based on these observations it is likely that rocks overlying the fault exert a higher-order control on the depth-dependent behavior of megathrust earthquakes properties, as proposed by Sallarès and Ranero (2019).

### 4.3. Co-Seismic Seafloor Deformation: Implications for Tsunamigenesis

Our modeling results also show that the low rigidity of the shallow domain causes large slip near the trench, which in turn enhances large co-seismic seafloor deformation. Rupture simulations in the shallow domain show that uplift (i.e., vertical displacement of the seafloor) in the shallow domain is up to 4–6 times larger in the heterogeneous model (uplift ratios between the heterogeneous and homogeneous models, $Uplift_{net}/Uplift_{hom}$, are shown in Figure 10a; absolute values of uplift from each model are shown in Figure S8 in Supporting Information S1). These uplift differences lead to important differences in tsunami amplitudes (Figures 10c, 10d and S9 in Supporting Information S1 and Movies S2 and S3).

We used the final uplift of the rupture through the shallow domain in both models (Figure S8 in Supporting Information S1) and the Tsunami-HySea software (Macías, Castro, & Escalante, 2020; Macías et al., 2017; Macías, Ortega, et al., 2020) to calculate the tsunamigenic response. Tsunami-HySEA solves the two-dimensional non-linear hydrostatic shallow-water equations in both Cartesian and spherical coordinates systems using second order finite volume methods implemented in CUDA for multi-GPU NVIDIA architectures. Tsunami-HySEA implements a two-step numerical scheme that guarantees the mass conservation in the complete domain and prevents the generation of spurious high-frequency oscillations. The Tsunami-HySEA code has been specifically designed for tsunami simulations and has been extensively tested and validated using specific benchmarks for tsunami propagation and inundation and tsunami currents (Macías, Castro, & Escalante, 2020; Macías et al., 2017; Macías, Ortega, et al., 2020). Our dynamic rupture simulations
<figure><img src="image_12.png" /><figcaption>Figure 10.</figcaption></figure>
have no topography. Thus, to capture tsunami height variations related to a more realistic topography we assume a simplified bathymetry that contains the main topographic features of subduction zones, including the continental shelf and slope, the trench and the outer rise (Figure 10a). The extent and slopes of each of these domains are based on bathymetry data from the Nicaragua convergent margin (e.g., Sallarès et al., 2013). The spatio-temporal evolution of the tsunami wave amplitude was measured at 10 different points of interest (POI) that recorded the tsunami propagation from the trench to the shore (Figure 10c). Figure 10d shows the maximum wave ratio ($Max_{het}/Max_{hom}$) at each POI, while the absolute values of the tsunami amplitude can be found in Figure S9 in Supporting Information S1. The maximum wave amplitude of the tsunami derived from the heterogeneous model is 6–8 times larger than that from the homogeneous model, and differences increase to ~10 at the shore, because of the interaction of the shallow bathymetry with the tsunami wave. As indicated in real-case tsunami earthquake studies (Geist & Bilek, 2001; Satake, 1994; Satake & Tanioka, 1999), these results demonstrate the importance of taking into account realistic low-rigidity upper-plate values (1–10 GPa) in the shallow domain as they determine large slip and important tsunamigenic uplift.

Rupture simulations of the shallow domain also reveal significant differences in co-seismic horizontal trenchward displacement of the hanging wall (i.e., along-dip trenchward displacement) between elastic models. Similar to uplift, along-dip trenchward displacement ratios ($Along-dip_{het}/Along-dip_{hom}$) increase toward the trench, where displacement may be up to 8–10 times larger in the heterogeneous model (Figure 10b). Horizontal trenchward displacement of the slope region may promote folding of accretionary structures (Tanioka & Seno, 2001) and reactivation of pop-up structures (Hananto et al., 2020), providing additional tsunami sources as indicated in previous simulations (e.g., Murphy et al., 2016).

## 5. Conclusions

We demonstrate that depth-dependent variations of upper-plate elastic properties along the megathrust fault exert a major effect on rupture dynamics and tsunamigenesis. We performed dynamic rupture simulations through homogeneous and heterogeneous elastic models, and a bimaterial fault. The assumed friction properties are the same in all simulations. We use the generic upper-plate depth distribution of $V_p$, $V_s$, and $\rho$ from Sallarès and Ranero (2019) to build the heterogeneous scenario, while we assumed the elastic properties of rocks overlying the fault at 25 km of depth to build the homogeneous one, and set a framework for comparing with elastic properties of the deep portion of the fault.

The results show that the trenchward decrease of upper-plate rigidity is a determining factor for enhanced slip, slip rate, slow rupture, and depleted high-frequency energy radiation at shallow depth. The depth distribution of slip ratio, rupture slowness ratio and $f_c$ ratio in our heterogeneous model quantitatively agrees with empirical estimates by Sallarès and Ranero (2019), thus the depth-dependence of upper plate rigidity explains most seismological observations regarding the depth-dependent behavior of tsunami earthquakes and large shallow megathrust earthquakes (e.g., Tohoku-Oki). Our model, accounts for first-order observations of depth-dependent behavior of megathrust earthquakes, and serves as starting point for the assessment of additional observations potentially related to complex, depth-dependent frictional scenarios.

Depth-dependent elastic properties also affect the dynamics of slip rate. Peak slip rate values in the heterogeneous model anticorrelate with rigidity variations. The increment in peak slip rate correlates with enhanced ground motion in the heterogeneous model, an effect that is amplified if we consider a bimaterial contrast across the fault. In addition, bimaterial simulations show that slip, rupture velocity and slip rate are largely controlled by the softer side of the fault, which according to geophysical observations is likely to be the upper plate side.

**Figure 10.** Ratios of (a) Uplift and (b) Horizontal (or Along-dip) trenchward displacement of the hanging wall between the homogeneous and heterogeneous models. The final uplift and along-dip trenchward displacement correspond to simulations through the shallow domain. (c) Close-up of the 3D perspective of the bathymetry used as setup (small inset) to simulate the tsunamis associated with the uplift of the shallow domain in the heterogeneous and homogeneous models. Red square in the inset shows the location of the close-up. White dots depict the location of 10 points of interest (POI) used to record the temporal evolution of tsunami wave (see Figure S6 in Supporting Information S1). (d) Maximum tsunami wave ratio versus time calculated for each POI.
The anomalously large slip in the shallow domain, together with an increasing compliance of the upper plate toward the trench, result in larger co-seismic seafloor deformation than in a homogeneous medium. Uplift differences between elastic models translate into order-of-magnitude differences in tsunami amplitude at the shore. The low rigidity at the toe of the upper plate also enhances horizontal trenchward displacement, which may contribute to amplifying the tsunamigenic response.

This study shows the importance of considering realistic variations in megathrust elastic properties and upper-plate rigidity in dynamic rupture simulations and in source models for tsunami modeling. Neglecting these properties may result in significant underestimation of slip, rupture time, local ground motion, seafloor co-seismic deformation and tsunami size, leading to underestimation of the tsunami hazard of the margin.

## Data Availability Statement

The results of the simulations through each fault domain, scripts needed to build the unstructured mesh, and the files containing the frictional setup can be found here: https://figshare.com/projects/Dynamic_ruptureSimulation_results_from_Prada_et_al_preprint_https_doi_org_10_1002_essoar_10507023_1_/_/120540.

## Acknowledgments

M. Prada is funded by the Beatriu de Pinós postdoctoral programme of the Government of Catalonia's Secretariat for Universities and Research of the Ministry of Economy and Knowledge (Ref # 2017BP00170). Jean-Paul Ampuero is supported by the Investments in the Future project UCAJEDI (ANR-15-IDEX-01) managed by the French National Research Agency (ANR). First author and ICM have had funding support of the "Severo Ochoa Centre of Excellence" accreditation (CEX2019-000928-S) of the Spanish Research Agency (AEI). The authors thank Florian Le Pape for his assistance during the generation of the bimaterial unstructured mesh. The authors thank the associate editor Yehuda Ben-Zion, and two anonymous reviewers for their comments that help to improve the manuscript. Dynamic rupture simulations in this study were done at the Super Computer Shaheen II at KAUST University. Shaheen II is a Cray XC40 delivering over 7.2 Pflop/s of theoretical peak performance. Overall, the system has a total of 197,568 processor cores and 790 TB of aggregate memory. Tsunami-HySEA code development is supported by the Spanish Government-FEDER funded project MEGAFLOW (RTI2018-096064-B-C21). This is a contribution of the Barcelona Center for Subsurface Imaging that is a Grup de Recerca de la Generalitat de Catalunya (2017 SGR 1662).

## References

* Ampuero, J. P., & Ben-Zion, Y. (2008). Cracks, pulses and macroscopic asymmetry of dynamic rupture on a bimaterial interface with velocity-weakening friction. *Geophysical Journal International*, 173(2), 674–692. https://doi.org/10.1111/j.1365-246x.2008.03736.x

* Andrews, D. J. (1976). Rupture velocity of plane strain shear cracks. *Journal of Geophysical Research*, 81(32), 5679–5687. https://doi.org/10.1029/jb081i032p05679

* Bilek, S. L., & Lay, T. (1999). Rigidity variations with depth along interplate megathrust faults in subduction zones. *Nature*, 400(6743), 443–446. https://doi.org/10.1038/22739

* Byrne, D. E., Davis, D. M., & Sykes, L. R. (1988). Loci and maximum size of thrust earthquakes and the mechanics of the shallow region of subduction zones. *Tectonics*, 7(4), 833–857. https://doi.org/10.1029/tc007i004p00833

* Calahorrano, A., Sallarès, V., Collot, J. Y., Sage, F., & Ranero, C. R. (2008). Nonlinear variations of the physical properties along the southern Ecuador subduction channel: Results from depth-migrated seismic data. *Earth and Planetary Science Letters*, 267(3–4), 453–467. https://doi.org/10.1016/j.epsl.2007.11.061

* Clayton, R., & Engquist, B. (1977). Absorbing boundary conditions for acoustic and elastic wave equations. *Bulletin of the Seismological Society of America*, 67(6), 1529–1540. https://doi.org/10.1785/bssa0670061529

* Contreras-Reyes, E., Maksymowicz, A., Lange, D., Grevemeyer, I., Muñoz-Linford, P., & Moscoso, E. (2017). On the relationship between structure, morphology and large coseismic slip: A case study of the Mw 8.8 Maule, Chile 2010 earthquake. *Earth and Planetary Science Letters*, 478, 27–39. https://doi.org/10.1016/j.epsl.2017.08.028

* Day, S. M., Dalguer, L. A., Lapusta, N., & Liu, Y. (2005). Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture. *Journal of Geophysical Research*, 110(B12). https://doi.org/10.1029/2005JB003813

* Duan, B. (2012). Dynamic rupture of the 2011 Mw 9.0 Tohoku-Oki earthquake: Roles of a possible subducting seamount. *Journal of Geophysical Research*, 117(B5). https://doi.org/10.1029/2011jb009124

* Galis, M., Pelties, C., Kristek, J., Moczo, P., Ampuero, J. P., & Mai, P. M. (2015). On the initiation of sustained slip-weakening ruptures by localized stresses. *Geophysical Journal International*, 200(2), 890–909. https://doi.org/10.1093/gji/ggu436

* Galvez, P., Ampuero, J. P., Dalguer, L. A., Somala, S. N., & Nissen-Meyer, T. (2014). Dynamic earthquake rupture modelled with an unstructured 3-D spectral element method applied to the 2011 M 9 Tohoku earthquake. *Geophysical Journal International*, 198(2), 1222–1240. https://doi.org/10.1093/gji/ggu203

* Galvez, P., Dalguer, L. A., Ampuero, J. P., & Giardini, D. (2016). Rupture reactivation during the 2011 M w 9.0 Tohoku earthquake: Dynamic rupture and ground-motion simulations. *Bulletin of the Seismological Society of America*, 106(3), 819–831. https://doi.org/10.1785/0120150153

* Geersen, J. (2019). Sediment-starved trenches and rough subducting plates are conducive to tsunami earthquakes. *Tectonophysics*, 762, 28–44. https://doi.org/10.1016/j.tecto.2019.04.024

* Geist, E. L., & Bilek, S. L. (2001). Effect of depth-dependent shear modulus on tsunami generation along subduction zones. *Geophysical Research Letters*, 28(7), 1315–1318. https://doi.org/10.1029/2000gl012385

* Hananto, N. D., Leclerc, F., Li, L., Etchebes, M., Carton, H., Tapponnier, P., et al. (2020). Tsunami earthquakes: Vertical pop-up expulsion at the forefront of subduction megathrust. *Earth and Planetary Science Letters*, 538, 116197. https://doi.org/10.1016/j.epsl.2020.116197

* Harders, R., Ranero, C. R., Weinrebe, W., & Behrmann, J. H. (2011). Submarine slope failures along the convergent continental margin of the Middle America Trench. *Geochemistry, Geophysics, Geosystems*, 12(6). https://doi.org/10.1029/2010gc003401

* Huang, Y., & Ampuero, J. P. (2011). Pulse-like ruptures induced by low-velocity fault zones. *Journal of Geophysical Research*, 116(B12). https://doi.org/10.1029/2011jb008684

* Huang, Y., Ampuero, J. P., & Helmberger, D. V. (2016). The potential for supershear earthquakes in damaged fault zones—theory and observations. *Earth and Planetary Science Letters*, 433, 109–115. https://doi.org/10.1016/j.epsl.2015.10.046

* Huang, Y., Ampuero, J. P., & Kanamori, H. (2014). Slip-weakening models of the 2011 Tohoku-Oki earthquake and constraints on stress drop and fracture energy. *Pure and Applied Geophysics*, 171(10), 2555–2568. https://doi.org/10.1007/s00024-013-0718-2

* Huang, Y., Meng, L., & Ampuero, J. P. (2012). A dynamic model of the frequency-dependent rupture process of the 2011 Tohoku-Oki earthquake. *Earth, Planets and Space*, 64(12), 1061–1066. https://doi.org/10.5047/eps.2012.05.011

* Hyndman, R. D., Yamano, M., & Oleskevich, D. A. (1997). The seismogenic zone of subduction thrust faults. *Island Arc*, 6(3), 244–260. https://doi.org/10.1111/j.1440-1738.1997.tb00175.x
Jeppson, T. N., Tobin, H. J., & Hashimoto, Y. (2018). Laboratory measurements quantifying elastic properties of accretionary wedge sed- iments: Implications for slip to the trench during the 2011 Mw 9.0 Tohoku-Oki earthquake. *Geosphere, 14*(4), 1411–1424. https://doi. org/10.1130/ges01630.1

Kanamori, H. (1972). Mechanism of tsunami earthquakes. *Physics of the Earth and Planetary Interiors, 6*(5), 346–359. https://doi. org/10.1016/0031-9201(72)90058-1

Kanamori, H., & Kikuchi, M. (1993). The 1992 Nicaragua earthquake: A slow tsunami earthquake associated with subducted sediments. *Nature, 361*(6414), 714–716. https://doi.org/10.1038/361714a0

Kanamori, H., & Yomogida, K. (2011). Preface. First results of the 2011 off the Pacific coast of Tohoku earthquake. *Earth, Planets and Space, 63*(7), 511. https://doi.org/10.5047/eps.2011.07.019

Kaneko, Y., Lapusta, N., & Ampuero, J.-P. (2008). Spectral element modeling of spontaneous earthquake rupture on rate and state faults: Effect of velocity-strengthening friction at shallow depths. *Journal of Geophysical Research, 113*, B09317. https://doi. org/10.1029/2007JB005553

Kodaira, S., Nakamura, Y., Yamamoto, Y., Obana, K., Fujie, G., No, T., et al. (2017). Depth-varying structural characters in the rupture zone of the 2011 Tohoku-oki earthquake. *Geosphere, 13*(5), 1408–1424. https://doi.org/10.1130/GES01489.1

Komatitsch, D., & Tromp, J. (1999). Introduction to the spectral element method for three-dimensional seismic wave propagation. *Geo- physical Journal International, 139*, 806–822. https://doi.org/10.1046/j.1365-246x.1999.00967.x

Koper, K. D., Hutko, A. R., Lay, T., & Sufri, O. (2012). Imaging short-period seismic radiation from the 27 February 2010 Chile (Mw 8.8) earthquake by backprojection of P, PP, and PKIKP waves. *Journal of Geophysical Research, 117*, B02308. https://doi.org/10.1029/2011JB008576

Kozdon, J. E., & Dunham, E. M. (2013). Rupture to the trench: Dynamic rupture simulations of the 11 March 2011 Tohoku earthquake. *Bulletin of the Seismological Society of America, 103*(2B), 1275–1289. https://doi.org/10.1785/0120120136

Lay, T., & Bilek, S. L. (2007). Anomalous earthquake ruptures at shallow depths on subduction zone megathrusts. In T. Dixon, & C. Moore (Eds.), *The seismogenic zone of subduction thrust faults* (pp. 476–511). Columbia University Press. https://doi.org/10.7312/dixo13866-015

Lay, T., Kanamori, H., Ammon, C. J., Koper, K. D., Hutko, A. R., Ye, L., et al. (2012). Depth-varying rupture properties of subduction zone megathrust faults. *Journal of Geophysical Research, 117*, B04311. https://doi.org/10.1029/2011JB009133

Li, J., Shillington, D. J., Saffer, D. M., Bécel, A., Nedimović, M. R., Kuehn, H., et al. (2018). Connections between subducted sediment, pore-fluid pressure, and earthquake behavior along the Alaska megathrust. *Geology, 46*(4), 299–302. https://doi.org/10.1130/G39557.1

Lomax, A., & Michelini, A. (2009). Tsunami early warning using earthquake rupture duration. *Geophysical Research Letters, 36*, L09306. https://doi.org/10.1029/2009GL037223

Lomax, A., & Michelini, A. (2011). Tsunami early warning using earthquake rupture duration and P-wave dominant period: The importance of length and depth of faulting. *Geophysical Journal International, 185*(1), 283–291. https://doi.org/10.1111/j.1365-246x.2010.04916.x

Lotto, G. C., Dunham, E. M., Jeppson, T. N., & Tobin, H. J. (2017). The effect of compliant prisms on subduction zone earthquakes and tsunamis. *Earth and Planetary Science Letters, 458*, 213–222. https://doi.org/10.1016/j.epsl.2016.10.050

Ma, S. (2012). A self-consistent mechanism for slow dynamic deformation and tsunami generation for earthquakes in the shallow subduc- tion zone. *Geophysical Research Letters, 39*, L11310. https://doi.org/10.1029/2012GL051854

Ma, S., & Beroza, G. C. (2008). Rupture dynamics on a bimaterial interface for dipping faults. *Bulletin of the Seismological Society of Amer- ica, 98*(4), 1642–1658. https://doi.org/10.1785/0120070201

Ma, S., & Hirakawa, E. T. (2013). Dynamic wedge failure reveals anomalous energy radiation of shallow subduction earthquakes. *Earth and Planetary Science Letters, 375*, 113–122. https://doi.org/10.1016/j.epsl.2013.05.016

Macías, J., Castro, M. J., & Escalante, C. (2020). Performance assessment of Tsunami-HySEA model for NTHMP tsunami currents bench- marking. Laboratory data. *Coastal Engineering, 158*, 103667. https://doi.org/10.1016/j.coastaleng.2020.103667

Macías, J., Castro, M. J., Ortega, S., Escalante, C., & González-Vida, J. M. (2017). Performance benchmarking of tsunami-HySEA model for NTHMP's inundation mapping activities. *Pure and Applied Geophysics, 174*(8), 3147–3183.

Macías, J., Ortega, S., Castro, M. J., & González-Vida, J. M. (2020). Performance assessment of Tsunami-HySEA model for NTHMP tsunami currents benchmarking. Field cases. *Ocean Modelling, 152*, 101645. https://doi.org/10.1016/j.oocemod.2020.101645

McIntosh, K. D., Silver, E. A., Ahmed, I., Berhorst, A., Ranero, C. R., Kelly, R. K., & Flueh, E. R. (2007). The Nicaragua convergent margin. In T. H. Dixon, & J. C. Moore (Eds.), *The seismogenic zone of subduction thrust faults, part III* (Chapter 9, pp. 257–287). Columbia Uni- versity Press. https://doi.org/10.7312/dixo13866-009

Meng, L., Inbal, A., & Ampuero, J.-P. (2011). A window into the complexity of the dynamic rupture of the 2011 Mw 9 Tohoku-Oki earth- quake. *Geophysical Research Letters, 38*, L00G07. https://doi.org/10.1029/2011GL048118

Moore, J. C., & Saffer, D. (2001). Updip limit of the seismogenic zone beneath the accretionary prism of southwest Japan: An effect of diagenetic to low-grade metamorphic processes and increasing effective stress. *Geology, 29*(2), 183–186. https://doi.org/10.1130/0091-7 613(2001)029<0183:ulotsz>2.0.co;2

Moreno, M., Melnick, D., Rosenau, M., Baez, J., Klotz, J., Oncken, O., et al. (2012). Toward understanding tectonic control on the Mw 8.8 2010 Maule Chile earthquake. *Earth and Planetary Science Letters, 321*, 152–165. https://doi.org/10.1016/j.epsl.2012.01.006

Murphy, S., Di Toro, G., Romano, F., Scala, A., Lorito, S., Spagnuolo, E., et al. (2018). Tsunamigenic earthquake simulations using experi- mentally derived friction laws. *Earth and Planetary Science Letters, 486*, 155–165. https://doi.org/10.1016/j.epsl.2018.01.011

Murphy, S., Scala, A., Herrero, A., Lorito, S., Festa, G., Trasatti, E., et al. (2016). Shallow slip amplification and enhanced tsunami hazard unravelled by dynamic simulations of mega-thrust earthquakes. *Scientific Reports, 6*(1), 1–12. https://doi.org/10.1038/srep35007

Newman, A. V., Hayes, G., Wei, Y., & Convers, J. (2011). The 25 October 2010 Mentawai tsunami earthquake, from real-time discriminants, finite-fault rupture, and tsunami excitation. *Geophysical Research Letters, 38*, L05302. https://doi.org/10.1029/2010GL046498

Ohnaka, M., Kuwahara, Y., & Yamamoto, K. (1987). Constitutive relations between dynamic physical parameters near a tip of the propa- gating slip zone during stick-slip shear failure. *Tectonophysics, 144*(1–3), 109–125. https://doi.org/10.1016/0040-1951(87)90011-4

Peter, D., Komatitsch, D., Luo, Y., Martin, R., Le Goff, N., Casarotti, E., et al. (2011). Forward and adjoint simulations of seis- mic wave propagation on fully unstructured hexahedral meshes. *Geophysical Journal International, 186*, 721–739. https://doi. org/10.1111/j.1365-246x.2011.05044.x

Polet, J., & Kanamori, H. (2000). Shallow subduction zone earthquakes and their tsunamigenic potential. *Geophysical Journal Internation- al, 142*, 684–702. https://doi.org/10.1046/j.1365-246x.2000.00205.x

Ranero, C. R., Morgan, J. P., McIntosh, K., & Reichert, C. (2003). Bending-related faulting and mantle serpentinization at the Middle America trench. *Nature, 425*(6956), 367–373. https://doi.org/10.1038/nature01961

Rubin, A. M., & Ampuero, J. P. (2007). Aftershock asymmetry on a bimaterial interface. *Journal of Geophysical Research, 112*(B5). https://doi.org/10.1029/2006JB004337
Sallarès, V., Melendez, A., Prada, M., Ranero, C. R., McIntosh, K., & Grevemeyer, I. (2013). Overriding plate structure of the Nicaragua convergent margin: Relationship to the seismogenic zone of the 1992 tsunami earthquake. *Geochemistry, Geophysics, Geosystems*, **14**(9), 3436–3461. https://doi.org/10.1002/ggge.20214

Sallarès, V., Prada, M., Riquelme, S., Melendez, A., Calahorrano, A., Grevemeyer, I., & Ranero, C. R. (2021). Large slip, long duration, and moderate shaking of the Nicaragua 1992 tsunami earthquake caused by low near-trench rock rigidity. *Science Advances*, **7**, https://doi.org/10.1126/sciadv.abg8659

Sallarès, V., & Ranero, C. R. (2019). Upper-plate rigidity determines depth-varying rupture behaviour of megathrust earthquakes. *Nature*, **576**(7785), 96–101. https://doi.org/10.1038/s41586-019-1784-0

Satake, K. (1994). Mechanism of the 1992 Nicaragua tsunami earthquake. *Geophysical Research Letters*, **21**, 2519–2522. https://doi.org/10.1029/94GL02338

Satake, K., & Tanioka, Y. (1999). Sources of tsunami and tsunamigenic earthquakes in subduction zones. *Pure and Applied Geophysics*, **154**, 467–483. https://doi.org/10.1007/978-3-0348-8679-6_5

Scholz, C. H. (1998). Earthquakes and friction laws. *Nature*, **391**, 37–42. https://doi.org/10.1038/34097

Scholz, C. H. (2019). *The mechanics of earthquakes and faulting*. Cambridge University Press.

Shlomai, H., Kammer, D. S., Adda-Bedia, M., Arias, R. E., & Fineberg, J. (2021). Unstable cracks trigger asymptotic rupture modes in bimaterial friction. *Journal of the Mechanics and Physics of Solids*, **149**, 104330. https://doi.org/10.1016/j.jmps.2021.104330

Simons, M., Minson, S. E., Sladen, A., Ortega, F., Jiang, J., Owen, S. E., et al. (2011). The 2011 magnitude 9.0 Tohoku-Oki earthquake: Mosaicking the megathrust from seconds to centuries. *Science*, **332**(6036), 1421–1425. https://doi.org/10.1126/science.1206731

Tanioka, Y., & Seno, T. (2001). The sediment effect on tsunami generation of the 1896 Sanriku tsunami earthquake. *Geophysical Research Letters*, **28**, 3389–3392. https://doi.org/10.1029/2001gl013149

Tappin, D. R., Grilli, S. T., Harris, J. C., Geller, R. J., Masterlark, T., Kirby, J. T., et al. (2014). Did a submarine landslide contribute to the 2011 Tohoku tsunami? *Marine Geology*, **357**, 344–361. https://doi.org/10.1016/j.margeo.2014.09.043

Tobin, H. J., & Saffer, D. M. (2009). Elevated fluid pressure and extreme mechanical weakness of a plate boundary thrust, Nankai Trough subduction zone. *Geology*, **37**(8), 679–682. https://doi.org/10.1130/g25752a.1

Udías, A., Vallina, A. U., Madariaga, R., & Buforn, E. (2014). *Source mechanisms of earthquakes: Theory and practice*. Cambridge University Press.

van Zelst, I., Wollherr, S., Gabriel, A.-A., Madden, E. H., & van Dinther, Y. (2019). Modeling megathrust earthquakes across scales: One-way coupling from geodynamics and seismic cycles to dynamic rupture. *Journal of Geophysical Research: Solid Earth*, **124**, 11414–11446. https://doi.org/10.1029/2019jb017539

von Huene, R., Klaeschen, D., Cropp, B., & Miller, J. (1994). Tectonic structure across the accretionary and erosional parts of the Japan Trench margin. *Journal of Geophysical Research*, **99**, 22349–22361. https://doi.org/10.1029/94jb01198