Slow slip events (SSEs) have been observed in spatial and temporal proximity to megathrust earthquakes in various subduction zones, including the 2014
Transient slow deformation of faults, slow-slip events, or silent earthquakes have been observed at convergent plate boundaries (Douglas et al., 2005;Dragert et al., 2001;Peng & Gomberg, 2010;Schwartz & Rokosky, 2007;Shelly et al., 2006) and at large continental faults, for example, the San Andreas fault (Linde et al., 1996;Rousset et al., 2019). Slow slip events (SSEs) may be accompanied by low-frequency seismic radiation, including tectonic tremors, low-frequency earthquakes, and very-low-frequency earthquakes (Khoshmanesh et al., 2020;Shelly et al., 2007). SSEs usually slip 10-100 times faster than the tectonic loading and last from days to years at depths close to the brittle-ductile transition (Dragert et al., 2001;Peng & Gomberg, 2010;Schwartz & Rokosky, 2007). The physical mechanisms underlying SSEs and their interaction with earthquakes are debated (Bürgmann, 2018):
The spatial viability of both fast and slow earthquakes on plate-boundary faults has been attributed to several factors, including structural and material heterogeneity (Lay et al., 2012;D. Li & Liu, 2016;Tobin & Saffer, 2009;Ulrich et al., 2022;Wang, 2010), rheological variability with depth (Gao & Wang, 2017;Saffer & Wallace, 2015) and fluid migration within oceanic sedimentary layers (W. B. Frank et al., 2015;Zhu et al., 2020).
Whether transient slow slip can serve as a universal precursor of eminent megathrust earthquake initiation is essential for seismic and tsunami hazard assessments in metropolitan margins (Bürgmann, 2018;Obara & Kato, 2016;Pritchard et al., 2020;Ruiz et al., 2014). However, the spatial and temporal interactions between slow and fast earthquakes, specifically the potential of slow-slip triggering megathrust earthquakes, remain enigmatic. Due to the observational challenges associated with the large variability of space and time scales, physics-based models are indispensable to illuminate the physics and in-situ fault properties, rendering SSE triggering of large earthquakes plausible.
On 18 April 2014, a M w 7.3 megathrust earthquake struck the coast of Mexico at the western edge of the Guerrero Gap, which had experienced no significant seismic events since 1911 (Kostoglodov et al., 1996;Radiguet et al., 2012). Geodetic inversions suggest that long-term slow-slip cycles have accommodated most of the plate convergence on the sub-horizontal oceanic slab between 20 and 45 km depth in Guerrero (Kostoglodov et al., 1996;Radiguet et al., 2012Radiguet et al., , 2016) ) (Figure 1a). In addition to long-term SSEs, transient bursts of short-term low-frequency earthquakes and tectonic tremors have been detected at different depths along the slab (Husker et al., 2012;W. B. Frank et al., 2015;W. Frank et al., 2015;Pérez-Campos et al., 2008). Slow-slip and slow earthquakes have been attributed to the elevated pore fluid pressure associated with an ultra-low velocity layer atop the subducting plate derived from dense-array seismic imaging (Song et al., 2009). Recent off-shore seismic observations have revealed a combination of co-seismic earthquake, aseismic and creeping deformation, suggesting the existence of multiple asperities across the slab interface (Plata-Martinez et al., 2021). Considering the unique slip characteristics of the Guerrero Gap, the initiation of the 2014 M w 7.3 earthquake has been related to the accumulated static Coulomb stress changes cast by an ongoing slow-slip event below 20 km depth that eventually accumulated an equivalent moment magnitude of M w 7.6 on the megathrust interface (Radiguet et al., 2016;Gualandi et al., 2017).
Integrated modeling of long-term tectonic loading and coseismic rupture advances the understanding of the dynamics of interseismic and coseismic slip, as well as their interplay (Cattania, 2019;Kaneko et al., 2011;Liu et al., 2020). While a few implementations have been developed to integrate long-term slow interseismic loading and fast coseismic rupture (Cattania & Segall, 2021;Segall et al., 2010;Yang & Dunham, 2023), they typically omit inertia effects during the interseismic period. Liu et al. (2020) couple two 3D finite element methods, one for long-term seismic cycle modeling and another for short-term dynamic earthquake rupture, linking stress and frictional parameters in geometrically simple setups. Cattania and Segall (2021) use 1D fractally rough faults and heterogeneous effective normal stress to model the spatiotemporal relationships between precursory slow slip and clusters of foreshocks. Due to algorithmic complexity and computational cost (e.g., Jiang et al., 2022;Lapusta & Liu, 2009;Thomas et al., 2014;Uphoff et al., 2023), it remains challenging to model the complete dynamics of 3D seismic cycles using a single code for a heterogeneous, geometrically complex subduction zone (see Text S1 in Supporting Information S1). Such modeling should also allow for observational data validation, as we undertake here.
In this study, we present 3D numerical models of the dynamic rupture of the 2014 M w 7.3 Guerrero earthquake, linked to 3D episodic slow-slip cycles under long-term tectonic loading, ensuring consistent stress states across the fault interface. Physics-based models of earthquake initiation, propagation, and arrest require choices regarding the pre-existing state of stress and fault strength governing frictional sliding (Harris et al., 2021;Oglesby & Mai, 2012;Ramos et al., 2021;van Zelst et al., 2019). Our SSE cycle and dynamic rupture models account for the same geophysical and geological observational inferences, such as the regional slab geometry, elevated pore fluid pressure, and depth-dependent frictional strength constrained from laboratory experiments and thermal modeling (Section 2). We bridge time scales from decades governing four episodes of long-term SSEs to fractions of seconds during earthquake rupture within the Guerrero Gap using the SSE cycle results to inform the dynamic earthquake rupture scenario models. The modeled, observationally constrained, transient stress evolution of the 2014 SSE event can lead to spontaneous co-seismic failure in the hypocentral region of the Guerrero earthquake. However, the episodic non-linear variability in shear stress caused by the three preceding SSEs, which correspond to the 2002, 2006, and 2009-2010 SSEs, remains too small compared to the high static fault strength required to match observations in the dynamic rupture model (Section 3). We also find that, in addition to SSE-induced stress heterogeneity, the complex propagation and arrest of the Guerrero earthquake require preexisting variable friction properties. Our study provides a mechanically self-consistent model for slow-slip triggered megathrust earthquakes and has important implications for the interaction between earthquakes and slow-slip in subduction zones and at large continental faults worldwide (Section 4).
We model episodic slow-slip cycles spontaneously emerging under long-term geological loading along the curved slab interface of the Guerrero Gap (Section 2.1). The long-term tectonic loading, which accumulates fault shear stresses, is balanced by the fault strength which is defined from a laboratory-derived rate-and-state friction , 2010). The so-called Guerrero Seismic Gap is a 100-km long segment between 100.2°W and 101.2°W (yellow bar) that lacks recent large earthquakes (Lowry et al., 2001). Purple shades indicate large (M w ≤ 6.8) earthquakes after 1940 (Lowry et al., 2001). The focal mechanism of the 2014 M w 7.3 Guerrero earthquake is shown in red (strike:304°, dip:21°, rake:99°, Global Centroid Moment Tensor catalog (GCMT) (Dziewonski et al., 1981;Ekström et al., 2012). A finite coseismic source model using teleseismic inversion is shown as yellow-to-red-to-black rectangles (Ye et al., 2016). The orange contours indicate the 10 and 20 cm aseismic levels of fault slip during the 2014 M w 7.3 slow-slip events (Radiguet et al., 2016). The blue triangles mark the permanent Global Positioning System stations used in a geodetic inversion of both the coseismic and slow slip (Gualandi et al., 2017). Depth contours from 5 km depth (trench) to 80 km depth are shown as dashed lines with 5 km depth spacing. Mexico City is shown in black. (b) Slab surface geometry extending from the trench to a depth of 60 km in both slow-slip cycle and dynamic rupture simulations. This slab geometry is inferred from the Middle America Seismic Experiment (MASE) (Pérez-Campos et al., 2008). We use the standard global projection WGS84/UTM, zone 11N to Cartesian coordinates. The detailed description of mesh generation and convergence analysis can be found in Text S2 in Supporting Information S1. Tetrahedral elements are color-coded by a 1D layered velocity model from seismic imaging (Dougherty & Clayton, 2014) that is used in the dynamic rupture model. law (Section 2.1.2). We constrain the fault frictional parameters by combining laboratory experiments on wet gabbro gouges (He et al., 2007) with a 2D steady-state thermal model constrained by P-wave seismic tomography (Manea & Manea, 2011). We extend a previous model that focused on the deeper part (10-60 km depth) of the slab covering episodic SSEs only (Perez-Silva et al., 2021). Here, we consider the geometrically complex slab up to the trench and thus include the entire seismogenic zone (5-60 km depth). We account for elevated pore fluid pressure atop the oceanic plate which locally reduces fault strength and eventually leads to episodic slow-slip emerging between depths of 20 and 45 km (Section 2.1.1, Figure 2). This elevation of pore fluid pressure has been suggested based on the seismically inferred high Vp/Vs ratios in central Mexico (Song et al., 2009) as well as in other subduction zones (Audet et al., 2009;Shelly et al., 2006).
Direct observations of slow-slip cycles are limited, motivating numerical simulations to elucidate the underlying mechanics of SSE and earthquake interactions. We simulate long-term slow-slip sequences on a convergent plate boundary and analyze the time-dependent evolution of slip rates and shear stresses on the fault interface in 3D (Figure 1b). We use a quasi-dynamic formulation and the Boundary Element Method. Our forward model adopts a laboratory-derived rate-and-state friction law and a 3D realistic subducting slab geometry beneath central Mexico. The governing equations relate the temporal shear stress evolution of an individual element in response to fault slip and long-term plate convergence following Rice (1993) as
where δ i (t) is the fault slip and K i,j is the shear stress in element j due to a unit dislocation in dip direction of element i. The static Green's function K i,j is calculated using triangular dislocations in a uniform half-space (Stuart et al., 1997) assuming a homogeneous shear modulus of μ = 30 GPa and density ρ = 2,670 kg/m 3 . The plate convergent rate V pl is set to be uniformly 61 mm/year based on a global plate motion model, the Relatvie Motion Model (DeMets et al., 2010).
We use the open-source code TriBIE (https://github.com/daisy20170101/TriBIE) (D. Li & Liu, 2016;Perez-Silva et al., 2021), which is parallelized with OpenMPI and has been verified in 2D and 3D community benchmark exercises (Erickson et al., 2023;Jiang et al., 2022). We here use the quasi-dynamic approach approximating inertia effects with radiation damping for our SSE cycle simulations. To this end, the radiation damping factor η = μ/(2c s ) (with c s being the shear wave speed) has been introduced (Rice, 1993). Compared to fully dynamic
10.1029/2023AV000979
simulations, the quasi-dynamic approach can lead to similar overall seismic cycle behavior but differing rupture dynamics (Jiang et al., 2022;Lapusta & Liu, 2009;Thomas et al., 2014). We detail all slow-slip cycle modeling parameters in the following.
Figure 2b shows the along-depth profiles of our assumed effective normal stress σn , pore fluid pressure (p f ), hydrostatic (0.37*σ z ) and lithostatic pressures (σ z ). We assume that lithostatic pressure is depth-dependent with a constant overburden gradient (i.e., σ z = ρg (-z)). The effective normal stress, defined as the difference between lithostatic pressure and pore fluid pressure, increases with depth at a constant gradient σn = 28 MPa/km until a depth of 2.7 km. At lower depths, effective normal stress remains constant as σn = 50 MPa except at the SSE source depth between 20 and 45 km. An effective normal stress of 50 MPa at seismogenic depth is a common assumption used in community benchmark studies (Jiang et al., 2022).
To reproduce the relatively low stress drops inferred for SSEs, we assume a low effective normal stress of σSSE n = 2.5 MPa at depths between 20 and 45 km based on our previous work for a narrower slab geometry (Perez-Silva et al., 2021) and linked to elevated pore fluid pressure. Such high, near-lithostatic pore fluid pressure is supported by the observed elevated ratio between V p and V s from seismic imaging along the coast of southwest Japan, Cascadia, and central Mexico (Audet & Burgmann, 2014;Song et al., 2009).
Fault shear strength in the quasi-dynamic SSE simulation is governed by a laboratory-derived rate and statedependent friction law, the aging law (Dieterich, 1979;Ruina, 1983). The effective friction coefficient f depends on the fault slip rate v and a single state variable θ as
Here, a and b are non-dimensional friction parameters for the direct effect and evolution effect, respectively, D RS is the characteristic slip distance over which θ evolves in response to velocity steps, f 0 is the friction coefficient at a reference velocity v 0 at steady state, and σn = σ np f is the effective normal stress, defined as lithostatic loading stress minus the pore fluid pressure.
At steady state θ = D RS /v, the friction coefficient is
. Slip remains stable, and any slip perturbation evolves toward a steady state when the friction stability parameter (ab) is positive (velocitystrengthening, VS). Slip can be either unstable or conditionally stable when (ab) is negative (velocityweakening, VW). We use uniform distributions for the initial slip rate V ini and the initial state variable θ ini on the entire fault.
We adopt the definition of the critical nucleation length h * RA based on the fracture energy balance for a quasistatically expanding crack (Rubin & Ampuero, 2005),
Here, we assume a shear modulus of μ = 30 GPa and Poisson's ratio of ν = 0.25. The ratio between the maximum width of the velocity-weakening portion of the slab and the critical nucleation length (h * RA ) significantly affects the slip behavior of modeled SSEs (Lapusta & Liu, 2009;Y. Liu & Rice, 2009) (D. Li & Liu, 2017;Perez-Silva et al., 2021).
For faults governed by rate-and-state friction, the quasi-static process zone at a non-zero rupture speed can be estimated as Λ 0 = C μ * D RS bσ n , where C is a constant of order 1 (Day et al., 2005;Jiang et al., 2022;Lapusta & Liu, 2009), μ* = μ for antiplane strain and μ* = μ/(1ν) for plane strain, where ν is Poisson's ratio. We note that our mesh size is considerably smaller than Λ 0 which ensures numerical stability and accuracy.
We adopt the empirical "aging" law that can be interpreted to account for time-dependent healing of microscopic stationary frictional contacts (Beeler et al., 1996, e.g.), for describing the temporal evolution of state variable (θ):
To regularize the solution at low slip rates we use the modification proposed by Rice and Ben-Zion (1996):
which is Equation 2 when V ≫ 0.
A distribution of (ab) at different temperatures has been obtained from laboratory experiments for wet gabbro gouges (He et al., 2007). We project this temperature-dependent (ab) distribution onto the slab interface using the thermal profile from a 2D steady-state thermal model constrained by P-wave seismic tomography in central Mexico (Manea & Manea, 2011). We assume a downdip transition temperature, (ab) = 0, of 415°C, which coincides with the maximum down-dip extent of long-term SSEs inferred from Global Positioning System (GPS) inversions (Radiguet et al., 2012). Velocity-strengthening conditions (ab) > 0 are imposed at the two lateral sides of the model domain to stabilize slip toward the plate convergence rate. The distribution of (ab) across the entire slab is shown in Figure 2a. The physical parameters including friction, initial stress, and elastic material properties aforementioned are listed in Table 1.
We use the open-source software SeisSol (https://github.com/SeisSol), which is based on the Arbitrary Highorder Derivative (ADER) Discontinuous Galerkin (DG) finite element method, to perform simulations of earthquake rupture dynamics and seismic wave propagation (Dumbser & Käser, 2006;Käser & Dumbser, 2006;Pelties et al., 2012). SeisSol has been optimized for modern high-performance computing architectures including
10.1029/2023AV000979
an efficient local time-stepping algorithm (Breuer et al., 2014;Heinecke et al., 2014;Krenz et al., 2021;Uphoff et al., 2017) and has been validated against several community benchmarks following the SCEC/USGS Dynamic Rupture Code Verification exercises (Harris et al., 2018;Pelties et al., 2014). Stress and particle velocities are approximated with 3rd-degree polynomials, yielding 4th-order accuracy in space and time during wave propagation simulation. We detail all dynamic rupture modeling parameters in the following.
We constrain the initial stresses in the dynamic rupture model from a snapshot of the shear and effective normal stresses across the fault interface in the 2014 SSE model. We track the traction ratio as the slow-slip fronts migrate along-strike and find that the local peak in the hypocentral region appears on day 317 (Figures 3f and 4a). This local peak of traction ratio is associated with the acceleration of the migrating front from 0.5 km/day to 3 km/day (Figures 4b and 4c). The shear traction and effective normal stress on day 317 of the 2014 SSE quasi-dynamic model are saved and spatially interpolated onto the higher-resolution dynamic rupture mesh of the subduction fault surface using the package ASAGI (Rettenberger et al., 2016). The resulting ratio between the initial shear and effective normal stress is shown in Figure 3f. The time-dependent evolution of the traction ratio parameter on the fault during the modeled SSE is shown in Movie S2.
We use a 1D depth-dependent model of the density and seismic velocities to set the elastic properties (μ and λ) in the dynamic rupture model, as shown in Figure S9 in Supporting Information S1 and Figure 1b. This 1D velocity model is based on seismic imaging of the central Mexico subduction zone (Dougherty & Clayton, 2014) using the Mapping the Rivera Subduction Zone (MARS) seismic array, which consists of 50 broadband seismic instruments with a station spacing of ∼40 km deployed from January 2006 to June 2007. This 1D layered velocity structure captures the major features of the subsurface (Kim et al., 2010;Song et al., 2009).
In the 3D dynamic rupture simulations, we adopt a linear slip-weakening (LSW) friction law to constrain the fault frictional strength which has been shown to largely depend on the fault slip distance in laboratory experiments (Ida, 1972;Palmer & Rice, 1973) (Table 2). LSW friction laws have been widely used in dynamic rupture simulations including models of large megathrust earthquakes such as the 2004 M w 9.1-9.3 Sumatra-Andaman earthquake (Ulrich et al., 2022;Uphoff et al., 2017), 2011 M w 9.0 Tohoku-Oki earthquake (Galvez et al., 2014), and rupture scenarios for the Cascadia subduction zone (Ramos et al., 2021). While SeisSol offers using various rate-and-state-friction laws, we opt for LSW friction due to its computational efficiency and fewer parameters. Although using rate-and-state friction as in the SSE cycle simulation may seem more consistent, differences in time stepping and time integration methods across numerical techniques can introduce inconsistencies as well (Liu et al., 2020).
Fault friction initial conditions are difficult to constrain on the scale of megathrust slip but play an important role in dynamic rupture nucleation and propagation (Ulrich et al., 2022;van Zelst et al., 2019). Based on several trial dynamic rupture scenarios we set the static friction coefficient to μ s = 0.626 and the dynamic friction coefficient to μ d = 0.546 within the assigned rupture asperities which yield realistic co-seismic rupture dynamics and arrest as well as spontaneous nucleation at a depth of 22 km due to the 2014 SSE stressing.
Our choice of static friction allows for a smooth nucleation process at the hypocenter without introducing additional overstress and is within the range of effective static friction typically used in dynamic rupture megathrust scenarios (Galvez et al., 2014;Madden et al., 2022;Ramos & Huang, 2019). We assume depth-dependent frictional cohesion c 0 and constant critical slip distance d c (Text S1 in Supporting Information S1).
We assume a statically strong fault (static friction coefficient μ s = 0.626) in agreement with the high static frictional strength of rocks (Byerlee, 1978) but effectively weakened by high pore fluid pressure. This specific choice of μ s allows us to model realistic co-seismic rupture dynamics and arrest, including realistic levels of slip, rupture speed, and stress drop, as well as spontaneous nucleation at 22 km due to the modeled 2014 SSE event.
The selection of dynamic friction is constrained by matching both the seismic source time function and the geodetic static surface displacements while ensuring a smooth rupture arrest. Figure S8 in Supporting Information S1 shows that the steady state rate-state friction at coseismic slip rates in the seismogenic zone is corresponding to the dynamic friction value in the LSW law. In our preferred model (referred to as Model A1), we include two asperities, constrained by the two peaks in moment rate function revealed in kinematic source inversion (Ye et al., 2016). We use a constant μ d within each asperity. An increase in μ d outside the asperities is required for smooth and spontaneous rupture arrest (Text S2 in Supporting Information S1). We find that by increasing μ d to values 30% (μ d = 0.826) higher than μ s , dynamic rupture gradually stops at the edges of the asperities. This setup results in a comparable duration and peak of moment release to teleseismic inversion (Ye et al., 2016) (Figure 6a). The on-fault distribution of μ d following 0.826-0.28 × G 1 (r 1 , r 2 ) is shown in Figure 6f.
We model cycles of long-term SSEs (Text S2 in Supporting Information S1) and select four sequential events that occur repeatedly every 4 years. During the 200-year simulation, the recurring times range between 1 and 5 years (Figure S9 in Supporting Information S1). Figure 3 shows snapshots of the fault slip rate in the modeled scenario of the 2014 SSE. Each SSE episode lasts for up to 12 months (Radiguet et al., 2012) and reaches a peak slip rate of up to 10 -6 m/s (Figures 3a-3c and 3e). Our numerical results match the region-specific source characteristics of long-term SSEs inferred from geodetic inversion using the regional GPS network (Radiguet et al., 2016) (Table S1 in Supporting Information S1). We attribute the good match of the first-order SSE characteristics to the realistic flat slab geometry and assumed near-lithostatic pore fluid pressure (D. Li & Liu, 2016;Perez-Silva et al., 2021). We select four sequential SSE episodes of our model, closely corresponding to the four geodetically recorded events in 2001/2002, 2006, 2009/2010, and 2014. We calculate the horizontal and vertical components of synthetic surface displacements at regional GPS stations and compare them with geodetic inversions (Gualandi et al., 2017;Radiguet et al., 2012). The comparison between the synthetic and observed GPS vectors during the 2014 SSE is shown in Figures 3g and 3h and for the three earlier SSE episodes in Figure S7 in Supporting Information S1. All modeled SSE events yield good agreement with geodetic observations, although only dip-slip is considered in our simulations (D. Li & Liu, 2016).
The 2014 SSE initiates simultaneously at the eastern and western edges of the Guerrero Gap at a depth of 40 km. Both slip fronts migrate toward the center at a rate of 0.5 km/day (Figures 3a and 4b). The megathrust slips at a higher rate after the coalescence of the migrating fronts in the center, and the SSE then bilaterally propagates across the entire fault between 25 and 40 km depth. However, we observe no immediate coseismic slip nucleating upon coalescence of the SSE fronts (between a depth of 20-45 km). This is different from the results of earlier 2D planar fault simulations (Kaneko et al., 2017) but in agreement with recent on-and off-shore observations that find no evidence of coseismic rupture due to collapsed slow-slip migrating fronts in the Guerrero Gap (Plata-Martinez et al., 2021).
Figure 4 shows the time-dependent evolution of the on-fault shear-to-effective-normal traction ratio and alongstrike migration speed during the cycle of all four SSEs. During the quasi-periodic emergence of the SSEs, we find that fault shear tractions overall increase down-dip of the seismogenic zone (below a depth of 20 km). However, this increase is not steady and varies considerably with the acceleration of the migrating slip fronts. The space-time evolution of the traction ratio, defined as the shear over effective normal stress during the modeled transient slip, is shown in Figures 3b-3d and 3f. Here, the traction ratio increases gradually from down-dip (30 km depth) to up-dip (20 km depth) and eventually reaches 0.64 in the hypocentral area of the 2014 M w 7.3 earthquake at a depth of 22 km, which is slightly shallower than that inferred by the USGS (Figures 3f and 4a).
The migrating 2014 SSE front moves slowly until day 267 and accelerates to 3.0 km/day at day 317 (Figure 4b). This acceleration, associated with rapid strain energy release, eventually increases shear stress at the down-dip end of the seismogenic zone in our model (see Figure 4c and Movie S2). The migration speed can vary depending on the temporal evolution of stress and stressing rate during the modeled SSE, which results in various values of traction ratio below the locked zone between different slow-slip cycles (Figure S5 in Supporting Information S1). Accelerating SSE fronts, as in our 2014 SSE model, have been observed before the 2014 Chile earthquake (Socquet et al., 2017a) and before larger earthquakes in Japan (Uchida et al., 2016), which was suggested as a potential precursory signal initiating megathrust earthquake nucleation.
In contrast, traction ratios increase considerably less during the earlier three modeled SSEs (blue lines in Figure 4a and blue-to-purple lines in Figure S4 in Supporting Information S1). Shear stresses temporally increase during the 2001/2002 and 2006 SSEs but decrease during the 2009/2010 event. For example, the peak traction ratio in the 2014 episode is about 3.23% higher than in the preceding 2009-2010 event, corresponding to a 0.1 MPa increase in shear stress. We highlight that the long-term increase of the peak traction ratio at the hypocentral depth during the 20-year-long simulation is small compared to the transient traction changes during the 2014 SSE (Figure 4a). None of the three earlier events leads to traction ratios large enough to overcome the (prescribed) frictional fault strength in the seismogenic part of the slab in our preferred dynamic rupture model.
In our 200-year long-term SSE cycle simulation there appear no earlier SSEs with comparable magnitude and recurrence intervals to our selected sequence and earlier transient stresses are insufficient to initiate a megathrust rupture in our model configuration (Figure S18 in Supporting Information S1). The long-term stress loading is accommodated by very long-term, low-amplitude slow slip episodes within the seismogenic zone. This modulates the stressing at seismogenic depth with a recurrence time of 100 years but causes no coseismic rupture (Figure S1 in Supporting Information S1). These long-lasting events accommodate a considerable fraction of the total accumulated strain within the shallow seismogenic zone, consequently limiting the shallow peak slip rates in the dynamic rupture simulation. The lack of shallow coseismic slip in our slow slip cycle simulation aligns with recent evidence for shallow fault creep off-shore Guerrero (Plata-Martinez et al., 2021). However, due to uniform plate loading rate and a lack of earlier geodetic constraints, we cannot rule out alternative models in which potential earlier SSEs may meet the megathrust's frictional yielding criteria.
We present the first 3D dynamic rupture model of the 2014 M w 7.3 Guerrero earthquake. Our rupture scenarios are informed by the transient stress of preceding SSEs when the peak of the traction ratio reaches the hypocenter (Figure 4) and additional predefined frictional heterogeneity on the fault. We focus on a preferred model (Section 2.2; Figure 5) which uses an LSW friction law (Andrews, 1985) to describe the co-seismic fault strength and yielding. The specific choice of a critical slip-weakening distance of d c = 0.05 m and a statically strong fault (static friction coefficient μ s = 0.626) ensures that the model that reproduces the key features of geophysical observations and provides physically self-consistent descriptions of earthquake initiation, dominantly governed by SSE-induced shear stress changes, and its dynamics and arrest, which are predominantly governed by predefined frictional asperities. We discuss alternative rupture scenarios, including one less realistic model with smaller μ s as shown in Figure S15 in Supporting Information S1, probing sensitivity to initial conditions in Section 4.2.
Although earthquake nucleation is linked to the transient stress of the SSE cycle, we show that capturing realistic rupture propagation and arrest requires additional heterogeneity of the megathrust slab. We show that including two circular frictional asperities (Section 2.2.3) can reproduce the observed co-seismic characteristics to firstorder. We vary the maximum possible frictional strength drop smoothly within each asperity: the dynamic friction coefficient μ d gradually increases at the edge of the asperities. High variability of dynamic friction has been reported in relation to fault materials and sliding rates in laboratory experiments (Collettini et al., 2019;Di Toro et al., 2004) and has been shown to largely affect coseismic rupture dynamics on crustal faults in numerical models (Aochi & Twardzik, 2020;Ramos & Huang, 2019).
10.1029/2023AV000979
In our earthquake model, self-sustained dynamic rupture nucleates spontaneously at a depth of 22 km, where the modeled 2014 SSE front acceleration leads to a local increase in shear traction (Figures 4a-4c). This location agrees with the observationally inferred hypocenters within their uncertainties (Figures 5a and 5b). Unlike typical dynamic rupture models, where nucleation is prescribed ad hoc (e.g., Galis et al., 2014), spontaneous runaway rupture is initiated merely by the locally increased shear stress of the preceding SSE transient. Our rupture model dynamically breaks the central asperity and subsequently migrates to the second patch under slightly increasing slip rates (Figure 5 and Movie S3). The rupture arrests smoothly at the boundaries of the prescribed frictional asperities. The final rupture area is located up-dip from the hypocenter and has no clear overlap with the area that hosts aseismic rupture during slow slip (Figure S10 in Supporting Information S1).
Our preferred earthquake simulation resembles the key observed seismic and geodetic characteristics within observational uncertainties (Figures 6a-6e). Two broad peaks in the moment release rate emerge in our dynamic rupture model, as inferred from teleseismic inversion using more than 70 stations across 35°-80°epicentral distance (Ye et al., 2016) (Figure 6a). This suggests a multi-asperity rupture process, including dynamic triggering and delays between different asperities (Figure 6f). In our rupture dynamics model, the first and second peaks appear closer in time than inferred in the inversions which may reflect additional complexities on natural faults and observational uncertainties. For example, the shape of the second asperity area may be varied in our
10.1029/2023AV000979 dynamic rupture model to better match the observed moment rate release timing. However, teleseismic inversion lacks the adequate resolution to better inform on the spatial extent of slip (Ye et al., 2016). Our modeled total cumulative moment release is 9.41 × 10 19 Nm, which corresponds to a moment magnitude of M w 7.28 and agrees well with the observations (Figure 6a). An alternative dynamic rupture model with only a single asperity (Section 4.2; Figure 7) fails to reproduce a realistic moment magnitude and the pronounced two-peak character of the moment rate release. Because both dynamic rupture models spontaneously initiate due to the same transient SSE stresses but strongly differ in co-seismic dynamics, we conclude that additional frictional heterogeneity is required to model the propagation dynamics and arrest of the Guerrero earthquake.
Geodetic inversion using permanent on-shore GPS stations yields smaller slip amplitudes (Gualandi et al., 2017) but a larger rupture area extending up to the trench, compared to teleseismic inversion (Ye et al., 2016) (Figures 6c and 6d). Similarly, our modeled dynamic rupture features shallow fault slip up-dip of the hypocenter, while our maximum slip amplitude is 2.5 m (Figure 6e), which is consistent with teleseismic inversion assuming V r = 2.5 km/s (Ye et al., 2016). We note that the differences in geodetic and teleseismic fault slip inversions are likely affected by limitations in data resolution and differences in the assumed source time functions, velocity et al., 2016) and SCARDEC (http://scardec.projects.sismo.ipgp.fr) (Vallee et al., 2011). (b) Mapview with horizontal surface displacements observed at continuous Global Positioning System stations (black Gualandi et al., 2017) and in our simulation (red). The red star marks the USGS catalog hypocenter. Accumulated fault slip from (c) regional geodetic inversion (Gualandi et al., 2017), (d) teleseismic inversion (Ye et al., 2016), and (e) preferred dynamic rupture scenario. The maximum slip is 0.25 m, 2.5 and 2.5 m, respectively. (f) Distribution of the prescribed heterogeneous dynamic friction coefficient μ d which gradually increases from 0.546 within to 0.826 at the edge of the asperities following an exponential function (see Methods: "Linear slip-weakening friction"). models, and/or fault geometries. Figure 6b shows modeled static surface deformation at 80 s after the rupture initiation and its comparison with geodetic observations (Gualandi et al., 2017). There are only two GPS stations (ZIHP and PAPA) with clear recorded signals close to the rupture area and one station (TCPN) with a smalleramplitude signal distant from the epicenter. Our synthetic surface displacements at ZIHP and PAPA are consistent with the reverse plate movement direction but slightly higher in amplitude than those observed.
Our preferred two-asperity dynamic rupture model reproduces both seismic and geodetic characteristics and is consistent with the localized slip heterogeneity inferred from seismic imaging using regional networks (Plata-Martinez et al., 2021;Song et al., 2009). Given the sparsity of co-seismic seismic and geodetic observations, we judge our forward model as data-justified first-order illumination of rupture dynamics and arrest. We note that future incorporation of a high-resolution regional velocity model, affecting the non-linear, coupled dynamics of rupture dynamics process and seismic wave propagation, may improve the achieved observational match.
We analyze the stress drop and energy budget of our preferred dynamic rupture model accounting for the preceding slow-slip cycle with respect to event-specific and global observations (Text S2 in Supporting Information S1). We calculate the average co-seismic stress drop in two different ways: (a) by spatially averaging the onfault stress drop, and (b) by averaging the modeled stress drop based on energy considerations (Noda et al., 2013;Perry et al., 2020). The two approaches result in average model stress drops of 1.74 and 2.1 MPa, respectively. These values are within the expected uncertainties (Abercrombie, 2021) of the seismological inference of 2.94 MPa (Ye et al., 2016) and are consistent with the global average of the inferred megathrust earthquake stress drops (Abercrombie & Rice, 2005).
Next, we analyze the earthquake initiation energy budgets accounting for the transient stress shadowed by the preceding SSE. We calculate the average fracture energy across the effective nucleation area directly induced by our modeled 2014 SSE in the hypocentral area as 0.17 MJ/m 2 (Text S2 in Supporting Information S1).
This inference is comparable to the range of nucleation energies (0.1-1 MJ/m 2 ) estimated for most M > 8 Nankai earthquakes in southwestern Japan (N. Kato, 2012), implying that the transient stresses of aseismic slip may play a ubiquitous role in the nucleation of megathrust earthquakes. In comparison, the dynamic rupture fracture energy averaged across the entire co-seismically slipping fault is only 0.11 MJ/m 2 . This is about 35% lower than the SSE fracture energy at the hypocenter governing the nucleation stage and similar to a seismologically inferred global average of 0.1-10 MJ/m 2 (Abercrombie & Rice, 2005), but 45% lower than the range of 0.2-2.0 MJ/m 2 measured on natural crustal faults (Tinti et al., 2005). This relatively low overall fracture energy is consistent with the low average stress drop, which results from the assumed elevated pore fluid pressure constrained by regional seismic imaging (Song et al., 2009). The elevated pore fluid pressure at depth is crucial for recovering faulting dynamics during both the long-term SSE and short-term initiation of our dynamic rupture model.
In addition to shear stress amplitudes, also the shear stressing rate increases significantly with increasing slip rate during the 4th SSE, and we observe a pronounced peak five days before the linking date (day 317, Figure S6 in Supporting Information S1). Shear stressing rates also change at the onset of the first and second SSE, but remain smaller or negative, and the peak amplitude of shear stress is lower during the 3rd event. Although temporal changes in shear stressing rate are not included in the dynamic rupture nucleation process, our linked model may suggest that the increasing stressing rate associated with the migrating fronts might be a proxy for an accelerating aseismic signal (Uenishi & Rice, 2003). have shown that the spontaneous nucleation governed by LSW friction is independent of the distribution of loading stresses or stressing rates as long as stress reaches the peak fault strength over a sufficiently wide region. However, the critical nucleation size of real events may depend on loading rate according to laboratory and numerical experiments using rate-and-state friction laws (Guérin-Marthe et al., 2019;Kaneko et al., 2008).
Our dynamic rupture models of the M w 7.3 Guerrero earthquake initiated by quasi-dynamic models of the preceding long-term SSE cycles illustrate the interaction between aseismic and co-seismic fault slip. It has been suggested that slow slip at the down-dip end of the seismogenic zone transfers shear stresses up-dip (Y. Liu & Rice, 2007) or temporally aid up-dip pore fluid migration (W. Frank et al., 2015), both of which potentially destabilize the locked portion of the megathrust, eventually triggering co-seismic rupture (e.g., Cattania & Segall, 2021) and increasing regional seismicity (e.g. Y. Liu & Rice, 2009). The kinematic migration patterns of off-shore aseismic slip are often challenging to constrain due to the lack of dense geodetic observations. Sequences of foreshocks and migrating seismicity before large events such as the 2011 Tohoku-Oki earthquake have been interpreted as proxies for aseismic fault slip and as potential long-term precursory signals of megathrust earthquake nucleation processes (A. Kato et al., 2012). Other observations of possible precursory signals include the acceleration of a M w 6.5 SSE that was recorded by the land-based GPS stations 8 months before the 2014 M w 8.1 North Chile earthquake (Socquet et al., 2017a).
We find that the transient increase in the shear-to-effective-normal-stress ratio resulting from the accelerating migration of the preceding slow-slip events can lead to the spontaneous initiation of realistic earthquake rupture and that this process is sensitive to the dynamics of the long-term transient SSE cycle. In our model, the increasing transient shear stress is sufficiently high for spontaneous dynamic rupture without additional weakening mechanisms, such as the effects of thermal pressurization (Noda et al., 2009). The total SSE-induced shear stress increase is ≈0.021 MPa, the difference between shear stress and yielding strength, in the hypocentral area. Figure S17 in Supporting Information S1 shows an alternative 3D dynamic rupture scenario in which, instead of using the transient stresses induced by slow slip, we prescribe an ad hoc time-dependent rupture initiation (following, e.g., Harris et al., 2018) as a weaker, spherical patch, centered at the hypocenter. The SSE transient stresses are not only large enough to nucleate earthquake dynamic rupture but also affect 3D rupture dynamics. Figures S17c and S17d in Supporting Information S1 show the resulting in shorter rupture duration, lower moment magnitude, and less complex moment rate release function due to reduced rupture complexity.
However, accounting for additional co-seismic weakening may further aid the slow-slip transient initiation of dynamic rupture (Hirono et al., 2016) inherently capturing our here prescribed variability of co-seismic frictional strength drop (Perry et al., 2020). Similarly, a recent conceptual model combining shallow SSEs and two asperities finds that the time-dependent balance between stress and strength is complex and not all SSEs directly lead to the nucleation of an earthquake (Meng & Duan, 2022), even when no geometrical complexity or pore fluid variation is considered.
For simplicity, we assume constant pore-fluid pressure during our geodetically constrained slow slip cycle modeling. Future work may explore the additional effects of dilatancy that may stabilize co-seismic slip (Segall et al., 2010) and may affect the overall slip budget at the downdip limit of the seismogenic zone (Y. Liu &
10.1029/2023AV000979 Rubin, 2010;Y. J. Liu, 2013). The effects of dilatancy and permeability enhancement in highly permeable fault zones may alter aseismic slip (Yang & Dunham, 2023). Dal Zillo et al., 2020 consider dilatancy to model SSEs in a planar Cascadia model and find slightly slower down-dip rupture speed and longer event durations, which may affect megathrust earthquake nucleation.
Accounting for megathrust asperities in our co-seismic dynamic rupture model is important for reproducing observationally inferred first-order source characteristics. Our preferred dynamic rupture scenario includes two frictional asperities (Figure 6f), which vary in their local dynamic friction coefficient from the surrounding slab interface, as proxies of megathrust heterogeneity governing the co-seismic rupture complexity. Simpler numerical model setups lend themselves to parameter space exploration (Y. Liu & Rubin, 2010;Ampuero & Rubin, 2008) While we here do not aim to cover the range of all possible initial condition variations in our complex model setup, we show two selected alternative dynamic rupture scenarios that illustrate the sensitivity of our SSEinitiated co-seismic rupture dynamics to prescribed frictional asperities. Our SSE cycle model is the preferred model out of five different long-term SSE cycle simulations (SI, Perez et al., 2019).
In dynamic rupture simulations, asperities due to locally reduced dynamic frictional strength lead to similar rupture behavior as asperities of elevated initial shear stress due to the equivalent fracture energy. Here, we present an alternative dynamic rupture model, Model A2, with a constant dynamic friction coefficient but heterogeneous initial shear stress. The initial shear stress is smoothly reduced outside both rupture asperities, which leads to spontaneous rupture arrest. We use the same spatial exponential function G 1 (r 1 , r 2 ) defined in Text S2 in Supporting Information S1 to decrease shear stresses smoothly outside the two geometrically equivalent preassigned rupture asperities. We set the initial shear stress as τ A2 0 = τ sse × G 1 (r 1 ,r 2 ) where τ sse refers to the onfault shear stress linked from the SSE cycle model (Figure S12a in Supporting Information S1). This setup leads to a localized distribution of the shear-to-effective-normal-stress ratio near the USGS catalog hypocenter (Figure S12b in Supporting Information S1).
The modeled source characteristics of the earthquake, including moment release, magnitude, slip distribution, and surface deformation, are all similar to our preferred model (Figure S14 in Supporting Information S1), except for a slightly sharper peak in moment release, corresponding to rupture arrest, than that of our preferred model (Model A1). We conclude that, in principle, local shear-stress asperities can lead to equivalent SSE-initiated rupture dynamics compared to frictionally-weak asperities.
Next, we demonstrate the sensitivity of rupture dynamics and synthetic observables (e.g., moment rate release) to megathrust heterogeneity using a single circular asperity wherein the dynamic frictional strength locally decreases (Model B1).
We examine a model with a single asperity with varying μ d on the fault. We manually introduce an exponential taper function, called G2 (r 1 ), similar to G 1 defined in Text S2 in Supporting Information S1 on the fault. The distribution of dynamic friction shaped according to function G2 is shown in Figure 7a. The resulting moment magnitude is only M w 7.15, and the moment rate release features a single sharp peak instead of reproducing the observed characteristic two-peak shape (Figure 7c). The modeled spatial extents of the fault slip and surface displacement amplitudes are significantly smaller (Figures 7b-7d).
10.1029/2023AV000979
Megathrust asperities have been related to depth-varying seismic and aseismic faulting behaviors (Lay et al., 2012;Walton et al., 2021). While we here parameterize both asperities as dynamically weak (low μ d ), heterogeneity in the initial stresses, structure, effective static fault strength, or pore fluid pressure (Bilek & Lay, 1999;Bürgmann, 2018) may serve as dynamically viable asperities (Harris et al., 2021;Madden et al., 2022;Ramos & Huang, 2019) and additional observations are required to distinguish between them. We show that local shear-stress asperities can lead to equivalent rupture dynamics in Section 4.2 and Figure S14 in Supporting Information S1. Our parameterization of frictional asperities is relatively simple but effective in reproducing firstorder characteristics within the uncertainties of sparsely observed earthquake kinematics. With improved observational coverage, better-constrained seismic and geodetic fault slip inversion may provide better information on frictional asperities. Such smaller-scale stress or frictional strength heterogeneity may lead to a more complex rupture process: Laboratory experiments, geodetic measurements, and seismological observations imply that additional small-scale heterogeneity and physical processes, such as variations in rheology (Gao & Wang, 2017), frictional properties (Lay et al., 2012), as well as pore fluid effects (Zhu et al., 2020) may impact the coseismic behavior. Denser regional seismic and geodetic instrumentation along the central Mexican coast and off-shore, allowing for better imaging of coseismic fault slip, would be crucial to inform and validate dataintegrated and physics-based modeling.
Our choice of frictional parameters in the dynamic rupture model allows for balancing the depth-dependent fault strength, heterogeneous initial shear stresses, and heterogeneous frictional strength drop to achieve realistic levels of coseismic slip and moment release across a relatively small rupture area in dynamic rupture simulations.
Varying the, in LSW friction well-defined, static friction coefficient impacts our match of the observed smooth acceleration in the moment rate function. Given the heterogeneous shear stress perturbation of the preceding SSEs, a well-defined yielding strength is helpful to understand spontaneous dynamic rupture nucleation to first order. This sensitivity is exemplified in Figure S15a in Supporting Information S1 where a slightly lower μ s results in delayed rupture arrest, a larger rupture area, and over-prediction of the amplitude and arrival of the first peak in the modeled moment release. Although simpler than the rate-and-state friction law used in the long-term SSE cycle simulations, we yield a similar range in reference friction coefficients (Figure S15b in Supporting Information S1) and comparable behavior in coseismic slip.
Our models help interpret geodetic and seismological observations of slow slip and coseismic megathrust rupture and help to unravel their interaction using available observations in Guerrero. We identify the acceleration of slow slip migration fronts as a driving mechanism preceding the initiation of coseismic rupture in our models. This may have important implications for enhancing the understanding of precursory slow slip, seismicity, and megathrust earthquakes in other subduction zones, such as in Japan (A. Kato et al., 2012). While our models do not enable the prediction of the relationship between long-term slow slip and future earthquakes, we anticipate our findings will also enhance the understanding of observed signals associated with the spectrum of megathrust faulting.
Our modeled SSE and coseismic fault slip are located largely off-shore in central Mexico, where a dense array of ocean bottom seismometers has discovered episodic shallow tremors, suggesting small-scale slow-slip events or low-frequency earthquakes (Plata-Martinez et al., 2021) potentially linked to small asperities up-dip of the slowslip region. Accounting for additional small-scale heterogeneity on the fault may help explain high-resolution observations, such as complexity in moment release rate and strong ground motions (Galvez et al., 2016) Here, we focus on the one-way interaction between the SSE cycle and dynamic rupture and omit the respective influence of coseismic rupture on slow-slip transients. Modeling 3D fully dynamic earthquake cycles on geometrically complex faults (Erickson et al., 2023;Jiang et al., 2022) that incorporate spontaneous (aseismic) nucleation, dynamic rupture, and post-seismic deformation are computationally challenging but are becoming achievable at realistic scales and levels of complexity to allow for direct observational verification. Extending our approach to a unified and fully coupled slow-slip and dynamic rupture framework is a promising future step.
We discuss the choice of linear-slip weakening friction in our dynamic rupture simulation by comparing key controlling factors of earthquake nucleation, the equivalent static friction coefficient (μ RS s ) and slip-weakening rate (W, as defined by Uenishi and Rice (2003)) between our 3D slow slip cycle and dynamic rupture models.
Coseismically, the slip-dependent fault weakening behavior governed by aging law rate-and-state friction is similar to that governed by LSW friction as has been shown in theoretical and numerical analysis (e.g., Bizzarri & Cocco, 2003;Garagash, 2021;Kaneko et al., 2008). We estimate an equivalent peak rate-and-state static friction coefficient μ RS s using the relation μ RS s ≈ f p = f 0 + a ln(V sr / V 0 ) (Garagash, 2021) and assuming slip rates ranging between 10 -9 -10 -7 m/s during the slow slip cycle simulation and a = 0.01. The such estimated peak value is μ RS s = 0.62, comparing well with μ s = 0.626 used in our linear-slip weakening dynamic rupture model.
Following Garagash (2021), we can estimate the equivalent linear-slip weakening Dc from aging law rate-andstate frictional weakening near the rupture front as D c ≈ 5.8 m, with constant b = 0.0135, σ = 50 MPa. We can also compare the equivalent critical slip distances assuming slip-law rate-and-state friction, following Uenishi and Rice (2003) by equaling the slip weakening rates for our frictional parametrizations of both models, defined as Δτ/ D c = W LSW and as b σ/ D RS = W RS with Δτ = (μ sμ d ) σ, which results in D c /D RS 5.93, implying an equivalent D c ≈ 1.5 m in linear slip weakening friction.
However, we find that our linked dynamic rupture model requires a small D c = 0.05 m (cf. Figure S16 in Supporting Information S1), resulting in a slip weakening rate of 77.9 MPa/m. This discrepancy may express different megathrust frictional behavior governing regions hosting SSE and dynamic rupture and could be further explored in future work, including additional physics or heterogeneity, for example, scale-dependent fracture energy (Gabriel et al., 2023;Ide & Aochi, 2005), alternative long-term friction evolution models (T. Li & Rubin, 2017), or analytically accounting for the rupture speed dependence of the aging law equivalent linear-slip weakening estimates. We note that matching dynamic friction may be less crucial since additional weakening mechanisms can be active at coseismic slip rates (e.g., Di Toro et al., 2011) and we caution that we here do not fully explore the effects of self-consistent parameterization on the interaction between slow slip and dynamic rupture simulations.
We simplify the complex physics and initial conditions in our models of SSEs and dynamic rupture in several ways. The long-term slow slip model initial conditions are not observationally constrained. Our model results in a series of quasi-periodic SSEs that vary considerably over time. For example, the recurrence intervals range between one and 5 years (Text S3 in Supporting Information S1). Our approach neglects the (small) volumetric stress changes induced by slow slip outside the megathrust interface, which may lead to inconsistencies when extending the linked dynamic rupture models to include off-fault plasticity (Ma & Nie, 2019;Ulrich et al., 2022) in future work. Although inertia effects of slow slip are expected to be mostly minor, the complex long-term stress evolution and short-lived changes in stressing rate that we find here motivate future work, for example, using an integrated dynamic switch between inter-and co-seismic stages (e.g., Liu et al., 2020).
By coupling porosity and permeability evolution to elastic fault deformation, Yang and Dunham (2023) demonstrate the potentially critical role of pore fluid transportation and permeability evolution on slow slip and seismic cycles in a 2D antiplane fault model. Using a two-phase flow model that couples solid rock deformation and pervasive fluid flow, dal Zilio et al. ( 2020) investigate the effect of poroelastic coupling on long-term fault evolution in a solid-fluid constitutive framework but restricted to 2D. Focusing on the geodetically-constrained SSE source characteristics and for computational efficiency, we here omit potential SSE-underlying poroelastic effects (e.g., Heimisson et al., 2019). These can be caused, for example, by the dynamics of fluid migration and pressure variations within porous materials and will be important to study, specifically in 3D, in future work.
We construct a 3D dynamic rupture model of the 2014 Guerrero earthquake initiated solely by a geodetically constrained long-term model of the 2014 SSE and not by three preceding events. Our chosen frictional parameters balance slow slip transient stressing with depth-dependent fault strength and frictional strength drop, resulting in realistic co-seismic dynamics, especially when compared to alternative models with differing friction coefficients. Our mechanically self-consistent and data-driven 3D models of long-term SSE cycles, megathrust earthquake initiation, and rupture dynamics in the Guerrero Seismic Gap contribute to a better understanding of the earthquake generation process. They can potentially lead to improved time-dependent operational earthquake forecasting (Uchida & Bürgmann, 2021). By incorporating the transient stress evolution of slow-slip before coseismic rupture and asperities in co-seismic friction drop, our models reproduce the kinematic and dynamic characteristics of both aseismic slip and co-seismic rupture and reveal their physical link. Although long-term
10.1029/2023AV000979 stress does not continuously accumulate, the accelerating migrating SSE fronts transiently increase shear stress at the down-dip end of the seismogenic portion of the megathrust. The SSE-induced transient stresses are not only large enough to nucleate megathrust earthquakes but also increase the complexity of 3D rupture dynamics. Improvements in the detection of transient aseismic slip deformation will aid in assessing seismic hazards in coastal regions (A. Kato et al., 2012;Socquet et al., 2017b). Furthermore, identifying distinct acceleration signals might be routinely possible in future regionally dense networks, specifically off-shore (Hilley et al., 2022).
2576604x, 2024, 2, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023AV000979, Wiley Online Library on [25/08/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
2576604x, 2024, 2, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023AV000979, Online Library on [25/08/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License