3D Modeling of Long‐Term Slow Slip Events Along the Flat‐Slab Segment in the Guerrero Seismic Gap, Mexico

During the last two decades, quasi‐periodic long‐term slow slip events (SSEs) of magnitude up to Mw7.5 have been observed about every 4 years in the Guerrero Seismic Gap, Mexico. We present numerical simulations of the long‐term SSE cycles along the 3D slab geometry of central Mexico. Our model accounts for the hydrated oceanic crust in the framework of rate‐and‐state friction and captures the major source characteristics of the long‐term SSEs occurring between 2001 and 2014, as inferred from geodetic observations. Synthetic surface deformation calculated from simulated fault slip is in good agreement with the cumulative GPS displacements. Our results suggest that the flat‐slab segment of the Cocos plate aids the large magnitudes and long recurrence interval of the long‐term SSEs. We conclude that 3D slab geometry is an important factor in improving our understanding of the physics of slow slip events.

The effect of slab geometry has been studied in previous non-planar rate-and-state friction modeling of SSEs in the Nankai, Cascadia, eastern Alaska and Hikurangi subduction zones (Li & Liu, 2016Li et al., 2018;Matsuzawa et al., 2013;Shibazaki et al., 2012Shibazaki et al., , 2019. Here we investigate the importance of 3D variations in slab geometry for the dynamics of SSEs in Guerrero. We perform numerical simulations of long-term SSEs using realistic 3D slab geometry of the Cocos plate in central Mexico. We compare modeled long-term SSE source characteristics and surface deformation with geodetic inversion results from twodecade-long continuous GPS records. Our model has important implications for our understanding of the physics of long-term SSEs in relation to slab geometry and fault strength. PEREZ-SILVA ET AL. Orange patches indicate the rupture area of major thrust earthquakes (adapted from Figure 1 in Radiguet et al. (2012). Black arrows show direction and rate (in cm/yr) of plate convergence (DeMets et al., 2010). The green contour represents the 20 cm cumulative slip contour of SSEs in 2001/2002from Radiguet et al. (2012. The blue contour denotes the location of SSEs in the Oaxaca segment from Fasola et al. (2016). Cyan patches show the regions of NVTs (after Figure 1 Husker et al. (2019)). The dashed red line shows the 4 mm slip contour of a short-term SSE (Villafuerte & Cruz-Atienza, 2017). Black dashed lines indicate the 10-km spacing depth contours of the Cocos plate from Pérez-Campos et al. (2008), with tags at every 20 km. Yellow triangles denote regional permanent GPS stations. The thick black line indicates the location of the trench; its magenta part highlights the along-strike extension of the model geometry used in this study and detailed in Figure 2. Mexico City is shown as a black square.

Methods and Model Setup
To model SSE cycles, we use the code developed by Li and Liu (2016). There are three main ingredients to this approach: (a) a quasi-dynamic formulation of traction and slip as defined by Rice (1993), (b) the fault constitutive response is given by a laboratory-derived rate-and-state friction law, and (c) it enables to incorporate non-planar, 3D fault geometry. We describe each ingredient in detail in the Supplementary Text S1.
We incorporate 3D plate geometry of the Cocos plate inferred from a dense broadband seismic array experiment (Pérez-Campos et al., 2008). The model domain extends 430 km along the Cocos-North America plate boundary and covers a depth range from 10 to 60 km ( Figure 2a). The slab is assumed locked from the trench to 15 km depth by imposing zero slip boundary conditions. While near-trench observations such as seafloor geodesy data is lacking, we base this assumption on the shallow locations of historical major thrust earthquakes (Figure 1) and the updip limit of inverted SSE distributions (Radiguet et al., 2016). We define a uniform plate convergence velocity that is directed N63°E at a rate of V pl = 6.1 cm/yr (DeMets et al., 2010). The slab geometry is discretized into triangular elements with edge lengths no longer than 1,500 m using the commercial software Cubit/Trelis (https://www.coreform.com/).
The friction parameter (a − b) is adapted from frictional experiments in wet gabbro gouges (He et al., 2007) and mapped onto the 3D slab interface (Figure 2b) based on the thermal structure beneath Central Mexico (Manea & Manea, 2011). Velocity-strengthening (VS) conditions (a − b > 0) are also imposed at the edges of the model domain to stabilize slip at the plate convergence rate ( Figure S1). The assumption of a VS band at the eastern end of the model is based on geodetic coupling estimates that indicate a low long-term coupling in the region eastward from the 2014 SSE (Radiguet et al., 2016). We demonstrate that the slip rate evolution does not change significantly in an alternative model setup without velocity-strengthening bands at the edges of the model domain (Text S3, Figure S8).
Following previous studies (e.g., Liu & Rice, 2007;Li & Liu, 2016;Shibazaki et al., 2019), we account for the inferred high pore fluid pressure condition atop the Cocos plate where SSEs occur (Jödicke et al., 2006;Kim et al., 2010;Song et al., 2009) by assuming low effective normal stress  within the SSE source region.  is set to be 2.5 MPa from 20 to 45 km depth and it is 50 MPa elsewhere. We refer to the depth range of low  as the "SSE zone" (Figure 2b). The model parameter W is defined as the downdip width of the SSE zone under velocity-weakening conditions (Figure 2b).
The ratio W/h*, where h* is the critical nucleation size given by the definition of Rubin and Ampuero (2005) (Text S1), has been shown to be a key parameter that controls the periodic behavior of the fault and the PEREZ-SILVA ET AL.  emergence of SSEs (Barbot, 2019;Liu & Rice, 2007). Previous studies have also shown that W/h* ∼ 1 reproduces slow slip characteristics in Cascadia (Liu & Rice, 2007;Li & Liu, 2016). Here we vary W and h* independently and report the source properties (recurrence, magnitude and slip velocity) for each parameter configuration ( Figures S3 and S4). Specifically, we explore various W by assuming different updip limiting depths of the SSE zone (Text S2). The simulation assuming W/h* = 1.18 (d c = 10.1 mm) is selected as the preferred model since it best reproduces the characteristics and geodetic signature of long-term SSEs in Guerrero. In the following we describe the results of this simulation. Other parametrizations and the respective sensitivities are discussed in Supplementary Text S2.
Of the preferred SSE cycle model we select four events as representative of the observed long-term SSEs in Guerrero and calculate their synthetic surface deformation. The required Green's functions are calculated assuming a homogeneous elastic half-space and implemented for triangular dislocation elements using the MATLAB code developed by Meade (2007).

Spatio-Temporal Evolution of Slip Rate
The model produces spontaneous slow slip events under constant plate loading until 145 years, at which time we terminate the simulation at the appearance of a seismic event (V > 5 mm/s). This earthquake's peak velocity occurs at 35 km depth (i.e., within the SSE zone) and is omitted in the following analysis since there is no observational evidence for such an event in the Guerrero Seismic Gap, yet. However, from previous planar fault SSE models (Liu & Rice, 2007) we expect that the major source characteristics of later SSEs would not be affected by the emergence of earthquakes in case these are separated by several years.
The slip rate of the modeled SSEs varies by several orders of magnitude on the fault (Figure 3a). Rupture migrating fronts where V > 3V pl (red contours in Figure 3a) indicate the occurrence of slow slip events. In the time interval between these events the fault is locked with slip velocities ranging from 0.03V pl to 0.3V pl (dark blue areas in Figure 3a).
In this model, we identify three types of SSEs with different along-strike extent. Type I events rupture most of the slab along-strike, extending approximately 300 km (green arrows in Figure 3a). The evolution of these events starts with a slow nucleation phase that takes place close to VS-VW transition at two distant locations (y = 150 km and y = −150 km), from which two slip fronts migrate toward each other (convergent yellow arrows in Figure 3a). These slip fronts merge into a velocity peak (V max > 10 3 V pl ) that occurs at ∼35 km depth; afterward they propagate bilaterally along-strike (divergent yellow arrows in Figure 3a) and alongdip (both updip and downdip). The downdip-migrating slip front propagates across the flat-slab segment. During this propagation phase, the peak velocities of the slip fronts progressively decrease. The migration fronts of these events are asymmetric, in that the slip front migrates slightly slower and over a longer distance eastwards than westwards.
Type II events (black arrows in Figure 3a) show a similar evolution pattern to Type I SSEs, except for a shorter along-strike extent (150-200 km) and a more symmetric migration path. Type III events have the shortest along-strike extent (<100 km) (red arrows in Figure 3a) and the lowest peak slip rates, V max < 10 2 V pl . Type II and Type III events may occur close to each other in space and time (e.g., the SSE during year ∼20 in Figure 3a) and in this case two distinct velocity peaks arise (Movie S3).
The peak slip rates during SSEs appear at different along-strike locations; peak velocities arise mostly in the eastern part of the fault (y < 0 km) for Type I and II SSEs and in the western part (y > 0 km) for Type III SSEs. This difference is important as most of the final slip accumulates during the period of peak slip velocities. Thus, the final slip of Type I and II SSEs concentrates in the eastern part of the fault, while for Type III it concentrates on the western part of the fault.
The magnitude and recurrence interval of peak slip rates also vary along the strike. The history of slip rates at two points (P1 and P2, colored circles in Figures 3a and 3d) is shown in Figure 3c, respectively. Peak slip velocities at P2 are 1-2 orders of magnitude lower than those at P1. The time interval between two sequential peaks is shorter at P2 (∼2-3 years) than that at P1 (∼4 years). The along-strike change in the slip rate evolution intensifies in the next 50 simulation years ( Figure S5a), where Type III SSEs persistently emerge Red stars indicate the data from four long-term SSEs in Guerrero (2001Guerrero ( /2002Guerrero ( , 2006Guerrero ( , 2009Guerrero ( /2010Guerrero ( , 2014 taken from Radiguet et al. (2012Radiguet et al. ( , 2016. Best fit scaling of modeled SSEs shown in black (M ∼ T 1.76 ). M ∼ T and M ∼ T 3 scaling are shown as reference.
in the western part of the fault (between 50 to 150 km along-strike); whereas Type I and II events concentrate further to the east.

Comparison With Geodetic Observations
We compare simulated SSEs in the preferred model with those estimated from geodetic inversion in terms of duration, magnitude and recurrence interval. We select the modeled SSEs in the time period between 10 to 50 years that occur within the GSG and calculate their source properties assuming a slip rate threshold of 10V pl (i.e., 61 cm/yr or 10 −7.7 m/s). The SSE duration is defined as the time period over which this slip velocity threshold is exceeded. We then calculate the total cumulative slip and moment magnitude within the estimated duration. The recurrence interval is given by the time between peak slip rates of successive SSEs. We assume a minimum slip of 1 cm to calculate the SSE magnitude, consistent with the threshold used in geodetic inversion (Radiguet et al., 2012).
Our modeled SSEs capture the major characteristics of the four long-term SSEs that occurred in 2001/2002, 2006, 2009/2010 and 2014, as inferred from geodetic inversion. They have an average duration, magnitude and recurrence interval of 8.7 ± 3 months, M w 7.44 ± 0.08, and 4.2 ± 0.2 years, respectively; all within the range of observations (Graham et al., 2016;Radiguet et al., 2012Radiguet et al., , 2016, as shown in Table S1. Figure S2 shows daily time series at GPS station CAYA, located along the Guerrero coast, and the cumulative slip at a fault node projected vertically from the station. The modeled recurrence interval agrees well with that indicated by the permanent geodetic records. To further validate the modeled SSEs, we show in Figure 4 the slip distribution of four modeled SSEs that best capture the characteristics of the four long-term SSEs in Guerrero. Movies S1-S4 show the spatial and temporal evolution of the slip rate on the fault during these events. The synthetic surface deformation of modeled SSEs (red arrows in Figure 4) is calculated as described in Section 2. To quantify the comparison between observed and modeled surface displacement, we define a misfit function as where  mod S and  obs S are the modeled and observed GPS displacement vectors at the jth station and N is the number of stations that detect the SSEs within our model domain. We report the respective misfits in Figure 4. The synthetic vectors match magnitude and direction of observations reasonably well, although the direction of the horizontal components along the coast shows a slight anti-clockwise rotational offset. The latter may indicate a strike-slip component, which contributes to the observed displacements. This component is not captured by our model, as we assume pure trench-normal slip.
The modeled SSEs exhibit along-strike migration rates of 0.5 ± 0.3 km/day, which is comparable to the slow migration speeds (0.8 km/day) reported for the 2006 SSE , but lower than the 6-9 km/ day estimated during the 2002 SSE (Kostoglodov et al., 2003). Low migration speeds, close to our model results, have also been reported for both observed and modeled long-term SSEs in Upper Cook Inlet in Alaska (Fu et al., 2015;Li et al., 2018), and in southwest Japan (Liu et al., 2010).

Long-Term Slip Budget
To estimate the long-term slip budget within the GSG from our model, we sum up the total cumulative slip released by long-term SSEs over 20 years and divide it by the total amount of slip accumulated due to plate convergence over the same period. The total slip released (Figure 3d) is calculated as follows: where N = 5 and δ i are the number of SSE episodes and the cumulative slip of each episode, respectively, and T = 20 years. The slip deficit equals to 1 − ξ.
We find that within the GSG, the fault releases >80% of the plate convergence loading via slow slip. This result is comparable to the geodetic estimate in Radiguet et al. (2012), which indicate that SSEs release 75% of the accumulated strain within the GSG over three SSE cycles. The slip budget varies minimally if we consider a time interval of 40 years to calculate the slip budget ( Figure S6).
PEREZ-SILVA ET AL.

Moment-Duration Scaling Relation
We calculate the moment-duration scaling relation of the modeled SSEs (triangles in Figure 3e) assuming a velocity threshold of 10V pl (i.e., 61 cm/yr) and a threshold slip of 1 cm to calculate the moment. The best-fit scaling of the modeled SSEs follows M ∼ T 1.76 . The moment and duration of the four long-term SSE episodes reported by Radiguet et al. (2012Radiguet et al. ( , 2016 fall within the upper bound of our model (red stars in Figure 3e). The wider range in magnitude and duration of modeled SSEs may result from the different spatio-temporal behaviors of all three types of SSEs as described in Section 3.1. This scaling relation changes only slightly when including more events with different h* values (M ∼ T 1.56 in Figure S7a). Assuming a noise-free slip threshold of 3 V pl , we find a slightly reduced exponent of 1.34 ( Figure S7b).

Geometric Effects on the Source Properties of SSEs
The emergence of long-term SSEs of large magnitudes, M w ≥ 7.0, observed along the flat-slab shallowly dipping segment beneath Guerrero suggests that variations in fault geometry may play a key role in understanding the variability of slow slip dynamics (e.g Maury et al., 2018). Our numerical findings support the importance of 3D slab geometry in reproducing the spontaneous emergence of realistic SSE cycles. In our model, the velocity-weakening portion of the fault under near-lithostatic pore fluid pressure conditions (defined as W) is inversely correlated to the average dipping angle at specific depths (20-45 km). As a result, W varies significantly (between 60 and 160 km) along strike, as shown Figure 3b.
Previously modeled SSE source characteristics (e.g., recurrence, slip rate, cumulative slip, etc.,) roughly scale with W/h* (Liu & Rice, 2009). In our preferred model, h* is kept constant along the entire slab and the relatively large W is a dominant factor that leads to large magnitudes and long recurrence interval characterizing long-term SSE dynamics in Guerrero. We perform additional simulations, which confirm the effect of W assuming constant h*. Increasing W by 6 km, which represents only a 4% increase with respect to its preferred value (W = 144.4 km), leads to a notable increase in the median periodicity, magnitude, and peak slip rate of the emerging SSEs ( Figure S4). Thus, even small changes in W/h* have an effect on the characteristics of modeled SSEs.
The lateral curvature of the slab influences the shear stress evolution on the fault. In previous models, larger SSEs tend to appear where the fault is flatter; and steepening of the slab promotes SSE arrest (Li & Liu, 2016). Here, we additionally find that the lateral patterns of the modeled SSEs, especially smaller events appearing on the western part of the fault, vary moderately over time; this can be seen by comparing the along-strike migration patterns of SSEs shown in Figure 3a and Figure S5. We note that this is a direct result of the complex evolution of shear stresses, since our simulation is independent of initial conditions after a short spin-up period (Rubin, 2008;Liu, 2014). This along-strike variation reflects an additional effect of the non-planar fault (Matsuzawa et al., 2013), as the western part of the fault has narrower W, which promotes more frequent and smaller SSEs (e.g., Liu & Rice, 2009). However, this influence on the lateral pattern is only moderate, as the termination of an individual modeled SSE is largely affected by the stress field left behind from previous events and is thus transient over time ( Figure S5b).

Implications for Diverse Slow Slip Along Central Mexico
In our preferred SSE cycle model, smaller SSEs nucleate northwest of the GSG (Figure 3a). In contrast, no SSEs have been observed in this region (∼101.5-103.5 W) between 2001 and 2014 (Maury et al., 2018). Recent time-dependent GPS modeling of the 2019 M w 7.0 SSE resolved aseismic slip starting northwest of the GSG (Cruz-Atienza et al., 2021), implying that this region may host slow slip.
To assess whether the modeled Type III SSEs would be detected given current GPS noise thresholds (<5 mm) we calculate the synthetic surface displacement of several events ( Figure S9). We find notable displacements of 10-25 mm at several close GPS stations. This inconsistency between observation and modeling may suggest that the region northwest of the GSG requires more complex frictional properties or pore pressure distribution than what our model assumes.
The emergence of SSEs also depends on the frictional properties and pore fluid pressure distribution on the plate interface, which in our model are assumed to be only depth-dependent. However, the heterogeneous distribution of the ultraslow velocity layer (USL) imaged atop the Cocos plate (Dougherty & Clayton, 2014;Song et al., 2009), implies that along-strike variations in pore fluid pressure exist beneath both the Guerrero and Oaxaca regions. Additionally, the degree of coupling of the plate interface, which also gives insight into the distribution of the frictional properties, also exhibits along-strike changes (Radiguet et al., 2012(Radiguet et al., , 2016. Our model does not include short-term SSEs associated with low-frequency earthquakes (LFEs) at the socalled sweet spot further down-dip (Frank et al., 2015;Husker et al., 2012), due to limited geodetic resolution by only two close continuous GPS stations (Figure 1). This along-dip variation of SSE recurrence may reflect the pore fluid increasing with depth modulated by temperature-dependent silica deposition as suggested by seismic imaging in northern Cascadia (Audet & Burgmann, 2014). Future application of our modeling approach may help to understand the along-dip variation in SSE source characteristics in future work.
Smaller and more frequent SSEs have been observed in Oaxaca (Graham et al., 2016), southeast of Guerrero (blue contours in Figure 1). The occurrence of these SSEs may be related to the narrower width of the USL in this region (Song et al., 2009), which results in a smaller W and thus in SSEs with lower magnitudes and shorter recurrence intervals ( Figure S4). The convergence rate of the Cocos plate under the North American plate, which increases southeastwards (DeMets et al., 2010), may also contribute to the shorter recurrence period of SSEs in Oaxaca, as this factor has been shown to inversely correlate with the recurrence interval of simulated SSEs (Li et al., 2018;Shibazaki et al., 2012;Watkins et al., 2015).

Implications for Source Scaling Relation
Our scaling falls in between the M ∼ T 3 scaling found for a wide range of regular earthquakes (Kanamori & Anderson, 1975) and the M ∼ T scaling inferred from a global compilation of SSEs (Ide et al., 2007). The differences in scaling relations between slow slip and regular earthquakes has been documented for many subduction zones (Ide et al., 2007;Gao et al., 2012;Peng & Gomberg, 2010) and is typically attributed to fundamental differences in the underlying physical mechanisms. For the four long-term SSEs inferred from geodetic inversion (Radiguet et al., 2012(Radiguet et al., , 2016, the observational scaling remains difficult to constrain due to the insufficient range of magnitudes and durations (Figure 3e).
It has been shown that simultaneous SSEs tend to have a different scaling relation than temporally non-overlapping, distinct SSEs, regardless of fault geometry and friction properties (Liu, 2014). Future modeling is required to include smaller SSEs further downdip to understand the factors controlling the observed scaling of SSEs in Guerrero. For instance, Rousset et al. (2017) constrained the moment-duration scaling of smaller SSEs and long-term SSEs and find that the moment scales with the duration with a power of 1-2, with the exponent being closer to 2 when accounting for the loading and release periods of long-term SSEs.
Recently, a cubic moment-duration scaling has been reported for deep SSEs in the Nankai (Takagi et al., 2019), Cascadia (Michel et al., 2019) and Mexico (Frank & Brodsky, 2019) subduction zones from geodetic and seismic observations. We note that an apparent shift of the scaling from M ∼ T to M ∼ T 3 may result from breaking a large slow slip event (as the 2006 SSE) into a cluster of daily slow transients calibrated by seismic LFE records (Frank & Brodsky, 2019). The identification based on cut-off slip rate may also considerably influence the geodetically resolved moment-duration scaling ( Figure S7b) (Li & Liu, 2017). Our results suggest that the separation between the two scaling relations may be not distinct. Rather, dynamic variability of natural fault slip (Peng & Gomberg, 2010) may also reflect in, potentially regional specific, continuous variability in SSE scaling relations.

Conclusions
We present the first 3D sequence simulations of long-term slow slip events in the Guerrero Seismic Gap, Mexico. Our model accounts for a realistic 3D fault geometry and laboratory-derived rate-and-state friction, and assumes the presence of high-fluid pressure regions atop the subducting slab at SSE source depths, supported by the existence of ultralow velocity layer revealed by high-resolution seismic imaging. The simulation produces spontaneously emerging long-term SSEs under constant geological plate convergence.
Our preferred model successfully reproduces the main source characteristics of long-term SSEs along the flat-slab segment beneath Guerrero as well as surface deformation obtained from continuous GPS measurements. In particular, we find that the source characteristics of the simulated SSEs agree well with those of the long-term SSEs. Our model results suggest that the unusually large magnitudes (M w ≥ 7.0) and long recurrence intervals (∼4 years) of SSEs in Guerrero are favored by the shallow-dipping flat-slab segment of the Cocos plate.
In addition, three distinct types of SSEs emerge in the model, which suggests that along-strike changes in the slab dip angle may affect the lateral distribution of SSEs. Modeled SSEs follow a moment-duration scaling of M ∼ T 1.76 , which is between the originally proposed linear scaling and the recently reported cubic relation for SSEs in Nankai, Cascadia and Mexico. Future work may be directed toward understanding the origin of the scaling trend of both long-term and short-term SSEs in Guerrero and Oaxaca, Mexico.

Data Availability Statement
The model parameter setup and simulation data is available at https://doi.org/10.5281/zenodo.4561753.

Acknowledgments
The work presented in this paper was supported by the