Geodetic Constraints on Recent Subduction Earthquakes and Future Seismic Hazards in the Southwestern Coast of Mexico

Three major subduction earthquakes occurred on March 20, 2012 (Mw 7.4), February 16, 2018 (Mw 7.2), and June 23, 2020 (Mw 7.4) in the southwestern coast of Mexico, which caused fatalities, casualties, considerable damage, and raised safety concerns about future seismic hazards. We use satellite geodetic observations to invert for the slip distributions of the three events and then investigate their interactions. Coulomb Failure Stress (CFS) induced by their slip both on surrounding thrust and normal faults are calculated. The positive CFS changes, along with the spatial‐temporal seismicity evolution, approximate earthquake recurrence rate and interseismic coupling, collectively indicate an increased possibility of a near‐future rupture around the areas between the 2018 and 2020 events in Oaxaca. Furthermore, there is a lowered chance of shallow coastal or offshore normal earthquakes but an increased chance of inland normal ruptures.

These subduction earthquakes play an important role in shaping the regional stress field and provide important insights into future seismic hazards. First, each of these events was preceded 17-55 years ago by a similar event that struck in the same area (Figure 1), namely the 1995 and 2012 events near Copala, the 1968 and 2018 events near PN and the 1965 and 2020 events east of SPP. It is therefore concerning that the most recent large subduction earthquake in the area between PN and SPP was only observed in 1978. Second, historical records reveal a notable seismic activity gap in the western part of the area affected by the 2012 event between Acapulco and Copala where no significant (Mw > 7.0) subduction events have occurred for at least 110 years (Kostoglodov et al., 2003). Kostoglodov et al. (1996) argued that this area had the potential to undergo an earthquake greater than Mw 8.0, but GPS measurements revealed that the SSEs might have released a significant part of the accumulated strain (Radiguet et al., 2012). Another seismic gap arises between the 2018 and 1978 events around Santa Catarina Juquila (SCJ) where the density of subduction events is low compared to its surroundings according to a detailed review of historical seismic activities by Sawires et al. (2019). Finally, apart from thrust-faulting earthquakes, the subduction of the Cocos plate has produced significant intraplate normal-faulting earthquakes occurring down-dip at some distance from the coupled subduction interface (e.g., the 1980 and 1973 events) as well as near the coast (e.g., the 1999 event east of SCJ). Apart from potential seismic hazards, the high angle normal faulting events such as the 1999 event may generate substantial tsunami and/or submarine landslide hazards. Whether, and to what extent, the recent large subduction events increased the possibility of these future hazards requires detailed investigation of their fault slip and the induced stress changes.
In this study, we utilize Radarsat-2 and Sentinel-1 Interferometric Synthetic-Aperture Radar (InSAR) observations ( Figure 2) to invert for the detailed fault slip distributions of the three subduction earthquakes mentioned above. Combining the regional seismic activity from 2000 to 2020 and the Coulomb stress changes induced by the fault slip of the three events, we investigate the interaction of the closely located 2012 and 2018 events, the spatial-temporal evolution of the regional stress changes and the potential future seismic hazards.
YU ET AL.

Data and Modeling Method
Line of Sight (LOS) surface displacements of Radarsat-2 and Sentinel-1 (Table S1) were generated using the GAMMA software and European Space Agency's precise orbits with the topographic phase contributions removed using the 3-arcsec (∼90 m) Shuttle Radar Topography Mission digital elevation model (Farr et al., 2007). To reduce the long wavelength effects, we applied the tropospheric delay corrections provided by the Generic Atmospheric Correction Online Service for InSAR (Yu et al., 2018) and then subtracted a best-fit quadratic surface computed for each InSAR pair (Feng et al., 2015;Sreejith et al., 2016). A spatial multi-look filter and a Goldstein spectral filter (Goldstein & Werner, 1998) were applied to reduce phase noise, leading to an interferogram resolution of about 90 m. A minimum-cost-flow unwrapping method was used to retrieve the LOS displacements (Chen & Zebker, 2000). To reduce the spatial correlation of the deformation map due to its high resolution, we down-sampled the interferograms using a quadtree-based algorithm (Simons et al., 2002), resulting in 1,530, 9,630 and 5,210 data points respectively for the 2012, 2018, and 2020 events.
The resultant displacements were input to PSOKINV (Feng et al., 2013), assuming an elastic half-space (Okada, 1986). We first determined the fault geometry associated with each of the three events nonlinearly using a particle swarm optimization method (Feng et al., 2013), with the initial geometries set according to the slab 2.0 model (Hayes et al., 2018). We then divided the fault planes into 2 km by 2 km patches to allow the slip to vary spatially on each fault segment. An iterative grid search method (Feng et al., 2020;Fukahata & Wright, 2008) was then employed to estimate the spatially variable slip distribution. To reduce the uncertainty associated with the dip angles derived from the previous nonlinear inversion, we tested different values ranging from 5° to 35° with a 1° step and determined the optimal dip angle and slip distribution simultaneously by minimizing data misfits ( Figure S1). We also implemented a Monte-Carlo analysis to quantify the uncertainties in the slip distributions using 100 different data sets, derived from the original interferograms randomly perturbed in a normal distribution about their original values using an a-priori standard deviation of 10 mm following Wright et al. (1999) (Figure S2). The resultant slip distributions (see YU ET AL.   the supporting information) were then used to calculate the Coulomb Failure Stress (CFS) in Coulomb 3 (Lin & Stein, 2004;Toda et al., 2005) with an effective coefficient of friction of 0.4 (e.g., Graham et al., 2014;Mikumo et al., 1999). We calculated the CFS respectively for reverse receiver faults at 15 km using an average fault geometry from the slab 2.0 model, normal receiver faults at 15 km for shallow coastal events (using the geometry of the 1999 event) and normal receiver faults at 40 km for inland events (using the average geometry of the 1980, 1973, and 2017 events according to Sawires et al. (2019) for the 1980 and 1973 events and Melgar et al. (2018) for the 2017 event). We then investigated how the CFS varies with depth by calculating seven cross sections parallel or perpendicular to the trench as shown in Figure 3c2. For reverse receiver faults, we used the average strike and dip angles of the slab 2.0 model along each profile. For normal receiver faults, we used the geometry of the 1999 event for Profiles Y1 and Y2, and the average geometry of the 1980, 1973, and 2017 events for Profile Y3 (see Table S2 in the supporting information).

Slip Distributions and Interactions
The resultant slip distributions are shown in Figure Figure 3b3 (black rectangles). Patches 1 and 2 clustered respectively at the northwest and southeast edges (along the strike direction) of the 2012 event slippage area (shown in Figure 4c along Profile Y2) suggesting intense stress concentration compared to its down-dip direction (Yamashita & Knopoff, 1987). Numerous small fractures around these tips or edges could interact through the stress fields perturbated by the mainshock and coalesced with one another because of the rapid increase of the stress-intensity factors as the cracks increase in size (McGarr, 1981). The occurrence of these aftershock clusters may therefore be explained by the sharpness of the rupture termination (Das & Henry, 2003). For example, the 1992 Landers earthquake had a sharp stress concentration at its southern end of the rupture where aftershocks clustered as opposed to its northern end. Such a phenomenon can also be observed in several other major earthquakes, for example, the Mw 8.1 1998 Antarctic Plate earthquake (Henry et al., 2000), the Mw 7.3 1992 Landers earthquake (Wald & Heaton, 1994) and the Mw 8.0 1985 Mexico earthquake (Mendoza & Hartzell, 1989). Furthermore, the along strike slip termination sharpness may reflect lateral variations of interseismic coupling (which increase gradually from the southeast to the northwest [e.g., Franco et al., 2012]) and/or irregularities in structural and geological properties of the fracture zones (Kostoglodov & Ponce, 1994).
The slip of the 2012 event resulted in substantial CFS increases as shown in Figure 3b1. However, we observed an area (∼5 km wide) between Patches 2 and 3 where no co-seismic slip was observed and the increased CFS led to no apparent increase of seismicity compared to its surroundings (Figures 3b3 and 3b4). Both of its left and right sides behaved more like velocity-weakening regions and released stress via unstable slips, especially Patch 3 (right side) where the densest aftershock following the 2012 earthquake was observed despite its greater distance from the epicenter (Figure 3b1). This implies the area may have acted as a barrier that terminated the 2012 rupture. This barrier has sharp stress concentrations as evidenced by the dense aftershocks at its edges and may indicate a substantial geometrical discontinuity or a change in the strength properties of the rocks. Aki (1979) noted that this kind of sharp stress concentrator could be a barrier between twin earthquakes which describes well the relationship between the 2012 and 2018 events. The stress concentrated at the barrier could be released through nonelastic deformation or elastically in a future earthquake.
The termination of the 2012 rupture at its southeast end (Patch 2) led to large dynamic stress perturbations in the surrounding medium (Bhat et al., 2007) and a substantial increase of CFS which propagated further east to initiate numerous aftershocks, notably Patch 3. These aftershocks released the concentrated stress but the YU ET AL.  resulting slip further altered the stress field of its neighboring fault asperities (Rudnicki & Kanamori, 1981) which may be responsible for initiating the 2018 event. This event ruptured the areas down-dip of the aftershock cluster (Patch 3 in Figure 3a2) and the change in CFS (ΔCFS) predicted well the aftershock distribution northeast of the epicenter (Figure 3b4). Considering most aftershocks occurred above the subduction interface (Figure 4a2), we also investigated the CFS on the optimally oriented fault planes by considering the regional stress field (Pacheco & Singh, 2010) (Table S3) instead of using a fixed fault geometry (King et al., 1994). The resultant CFS change ( Figure S3) has a similar pattern to that in Figure 3b5 and shows a clear relationship between the aftershock distribution and positive CFS changes. The lack of seismicity north of the 2012 event and northwest of the 2018 event can be explained by frequent SSEs (Figure 3b5) such as in 2011-2012 (Colella et al., 2017), 2014 (Radiguet et al., 2016), and 2017-2018 (Maubant et al., 2020). The maximum slip of the 2011-2012 SSE around SCJ reached ∼0.1 m, which is not enough to accommodate the long-term plate convergence rate (Colella et al., 2017;Correa-Mora et al., 2008). Instead, the SSEs here are generally confined to deeper source regions which will further load the seismogenic zone upward (Colella et al., 2017).
The 2020 event ruptured the east of the SCJ region in an area with no significant CFS increase as shown in Figures 3b1 and 3b2 (cumulatively <0.1 bar), suggesting the static ΔCFS produced by the 2012 and 2018 events was unlikely the triggering mechanism. The seismic activity around the 2020 epicenter only YU ET AL.  experienced a slight increase after the 2012 earthquake ( Figure 3b1) and no substantial increase after the 2018 earthquake (Figure 3b2). However, the dynamic stress changes may have played a role in remote triggering as the passage of transient seismic waves produced by the surrounding events could have initiated a secondary mechanism (such as viscous flow, fluid flow, and afterslip) that promotes the 2020 rupture (Freed, 2005;Kilb et al., 2000). A similar delayed dynamic triggering mechanism has previously been proposed to be responsible for the triggering of the 1999 Hector Mine earthquake by the 1992 Landers earthquake (Kilb, 2003) and the 2002 Tonga deep earthquake sequence by a magnitude 7.6 earthquake occurred at a depth of 598 km (Tibi et al., 2003). No sharp slip termination and clear aftershock clusters were observed to date around the 2020 rupture probably due to the short period of observation (5 months), the pre-existence of stress shadows caused by the 1965 event ( Figure 1) and/or SSEs such as in 2007 (Graham et al., 2016).

Future Seismic Hazards
One of the intuitive ways to investigate future seismic hazards is to calculate the ΔCFS induced by recent large earthquakes. Positive ΔCFS can raise the applied tectonic loading of a fault and move the fault closer to failure (Kilb et al., 2002). Numerous studies have demonstrated the correlation between positive ΔCFS and increased seismic activity or negative ΔCFS with zones of post-earthquake quiescence (e.g., Árnadóttir et al., 2004;King et al., 1994;Lin & Stein, 2004;Pope & Mooney, 2020). The CFS changes can also affect magmatic systems by triggering volcanic unrest and eruption (Biggs et al., 2013;Jenkins et al., 2021) and increase crustal permeability by channeling fluids to the seafloor (Bonini, 2019). Figures 3b5 and 4a, the subduction earthquakes tended to increase the CFS, calculated on the reverse receiver faults, in areas immediately surrounding the major slipping patches, extending both upward and downward along the subduction slab (Profiles X1, X2, and X4) as well as to the southeast and northwest along the trench (Profile Y2). This overlaps the seismically dangerous Guerrero seismic gap (Figure 1) whose accumulated strain is comparable to at least Mw 8.0 and has not broken since 1911 (Kanamori et al., 1993). Profile X3 in Figure 4a3 also shows that most of the subduction slab near SCJ was loaded, especially the portion between 0-30 km (∼0.3 bar increase). These positive ΔCFS could be responsible for the increase of the seismicity far from the mainshock epicenters such as the offshore events and those around SCJ. Most events are shallower than 30 km as shown in Figure 4a3, suggesting the region exhibits velocity-weakening frictional behavior and is capable of holding large earthquakes (Scholz, 1998). Radiguet et al. (2016) proposed that a 0.3-0.5 bar CFS increase induced by the SSE around the same area has triggered the 2014 Mw 7.3 Papanoa earthquake. Experiencing three subsequent positive ΔCFS perturbations, notably by ∼0.24 bar at a depth of 15 km and ∼0.28 bar at 25 km after the 2018 event and cumulatively ∼0.3 bar at 25 km (Figure 3d), the fault segments around SCJ have been continuously brought closer to failure.

As shown in
Another factor related to seismic hazard is the earthquake recurrence rate. The last recorded large subduction event around SCJ was 42 years ago (Sawires et al., 2019), which is spatially closer to the 2020 event although the uncertainty of its epicenter is high due to the lack of observations. Nishenko (1991) estimated the recurrence time of large subduction earthquakes in this region to be 54 ± 8 years based on historical earthquakes in 1870 (Ms 7.9), 1928 (Ms 7.8), and 1978 (Ms 7.8) which ruptured the same fault segment according to intensity and instrumental data. This estimate agrees well with the interseismic coupling inferred from GPS measurements (Correa-Mora et al., 2008;Scholz & Campos, 2012) and a review on recent subduction earthquakes by Sawires et al. (2019). These coherent seismic and geodetic observations indicate a fully coupled subduction interface around SCJ, which could result in repeating earthquakes of very similar size (Scholz & Campos, 2012). The recurrence rate of large subduction earthquakes inferred from the most recent data as reviewed by Sawires et al. (2019) in our study region increases from northwest to southeast with SCJ in the middle (17 years around Copala, 50 years around PN and 55 years around SPP as in Figure 1), suggesting a 42-year cycle is a sensible approximation in SCJ. Considering there were fewer subduction events around SCJ (Oaxaca gap in Figure 1) compared to its surroundings (Sawires et al., 2019), the continued strain accumulation on the subduction slab portion at this stage will no doubt raise safety concerns.
Seismic refraction and gravity data revealed a strong positive seismic primary wave velocity gradient from northwest to southeast along the middle America trench (Scholz & Campos, 2012;Valdes et al., 1986), resulting in a laterally segmented subduction plate with staggered rupture cycles. This explains the existence of the abovementioned barrier between Patches 2 and 3 (Figure 3b3), the lateral sharp slip termination and the seismic gaps. Conversely, considering the short temporal and spatial intervals between the 2012 and 2018 events, their ruptured asperities may have synchronized into a longer cycle in which both asperities would be ruptured simultaneously, breaking the barrier in between and resulting in a larger magnitude due to the increased rupture size. For example, this has been observed in the Kurile subduction zone where several neighboring Mw 7.5-7.8 earthquakes repeating every 30-50 years finally synchronized resulting in the 2011 Mw 9.1 earthquake (Hashimoto et al., 2009;Scholz & Campos, 2012). However, the synchronization could be obstructed due to the narrow downdip width of the seismogenic zone in this region (<80 km), given that Herrendörfer et al. (2015) proposed that cycle synchronizations tend to occur in regions with wider widths (>100 km).
We also calculated ΔCFS on normal receiver faults for coastal and inland earthquakes at 15 and 40 km depths ( Figure 3c) and their CFS cross sections along Profiles Y1, Y2, and Y3 ( Figure 4b). Overall, the recent subduction events generated stress shadows over coastal and offshore at shallow depths (0-28 km, Profile Y1 in Figure 4b2) which reduce the possibility of shallow normal earthquakes. Otherwise, those events could be potential tsunami generators (Sugawara et al., 2019) and trigger strike-slip earthquakes in the volcanic arc (Martínez-Díaz et al., 2004). The CFS along Profile Y2 shows complicated features where there are localized positive ΔCFS immediately below the ruptured faults and at shallow depths (0-10 km, Figure 4b3). We also observe 0.2-0.5 bar of CFS increases inland (Profile Y3 in Figure 4b4) which could pose a potential risk of future normal earthquakes. For example, Segou and Parsons (2018) proposed that an average of 0.5 bar CFS increase induced by the afterslip of the 2012 subduction earthquake was enough to trigger the 2017 inland normal earthquake (Figure 1). However, it is worth noting that the 2017 rupture may have already released the increased CFS caused by the 2012 event and generated sufficient stress shadows to compensate the stress increase caused by the 2018 event.

Conclusions
We used geodetic observations to invert for the slip distributions of three subduction earthquakes in the southwestern coast of Mexico along the Middle America Trench. Numerous aftershocks clustered at both the northwest and southeast slip patch edges of the 2012 event, but there is a lack of co-seismic slip and aftershocks, despite positive ΔCFS, between the 2012 and 2018 events where there may be a sharp barrier terminating the 2012 rupture at its southeast end. This termination triggered dense aftershocks further east which favored the subsequent 2018 rupture. ΔCFS calculated from the slip distributions of the three earthquakes predicted well the aftershock distribution and provided important insight into the region's future seismic hazard. The positive ΔCFS on reverse receiver faults induced subsequently by the three subduction earthquakes, along with evidence by investigating the regional subduction event recurrence rate, the interseismic coupling, and the level of the seismic activity, collectively indicate an increased possibility of a major future rupture around SCJ between the 2018 and 2020 events. There are also slight CFS increases over the Guerrero seismic gap, which has observed no major rupture (Mw > 7.0) for at least 110 years. Furthermore, ΔCFS calculated on normal receiver faults indicate a lowered chance of shallow coastal or offshore normal earthquakes occurring, but an increased chance of inland intraplate normal ruptures such as the destructive 2017 Puebla earthquake close to Mexico City.