# Interseismic coupling along the Mexican subduction zone seen by InSAR and GNSS

Louise Maubant, Mathilde Radiguet, Erwan Pathier, Marie-Pierre Doin, Nathalie Cotte, Ekaterina Kazachkina, Vladimir Kostoglodov

► To cite this version:

Louise Maubant, Mathilde Radiguet, Erwan Pathier, Marie-Pierre Doin, Nathalie Cotte, et al.. Interseismic coupling along the Mexican subduction zone seen by InSAR and GNSS. Earth and Planetary Science Letters, 2022, 586, pp.117534. 10.1016/j.epsl.2022.117534. hal-04986732

HAL Id: hal-04986732

https://hal.science/hal-04986732v1

Submitted on 24 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Interseismic coupling along the Mexican subduction zone seen by InSAR and GNSS

Louise Maubant¹,², Mathilde Radiguet¹, Erwan Pathier¹, Marie-Pierre Doin¹,
Nathalie Cotte¹, Ekaterina Kazachkina³, Vladimir Kostoglodov³

¹Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, IRD, Université Gustave Eiffel,
ISTerre, Grenoble, France

²Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology,
Cambridge, MA, USA

³Institute of Geophysics, National Autonomous University of Mexico, Mexico City, Mexico

## Keywords:

* InSAR time series

* Interseismic coupling

* Mexican subduction

Corresponding author: Louise Maubant, maubantl@mit.edu
**Abstract**

In the Mexican subduction zone, the heterogeneous spatial density of the GNSS network is limiting the capability to measure and study the variability of the surface deformation field all along the subduction zone. To improve the spatial coverage of interseismic deformation estimates, we combine InSAR and GNSS data. We use four Sentinel-1 descending tracks covering a large part of the Mexican Pacific coast from Jalisco to Oaxaca, over the period from 2016 to 2019.5. The InSAR and GNSS time series, over the same 3.5-year period, are corrected for the coseismic deformations of two $M_w > 7$ earthquakes but not for slow slip events (SSE) and postseismic surface displacements. We present estimates of mean surface velocities over the study period, which combine the linear inter-SSE trend and transient aseismic slip. By inverting these velocities, we estimate the short term coupling at the plate interface, which shows large lateral variations. InSAR data significantly improves the resolution of the inversion, compared to a model based only on GNSS data. Although not representative of long-term coupling in areas with large transient events, the short-term coupling estimates highlight the role of transient signals, such as slow slip events, on the coupling variability measured by geodesy. In agreement with GNSS time series over longer time period, the data also suggest that the Jalisco area (longitude from -106°E to -104°E), where no large transient events are observed, is highly coupled with a coupling ratio higher than 0.8.

**1 Introduction**

Interseismic deformation in continental domains and in subduction zones can be highly variable in time and space (e.g. Rousset et al., 2016; Radiguet et al., 2012; Nishimura et al., 2004). In some subduction zones such as Mexico, a succession of loading and releasing periods have been observed during the interseismic phase at different timescales (daily to several months), corresponding to the occurrence of slow slip events (SSEs). Taking into account these phases of aseismic slip in the estimation of the coupling during the interseismic period changes the average loading rate (e.g. Radiguet et al., 2012; Wallace et al., 2012).
Because the SSEs surface displacements are in a direction opposite to those caused by elastic
loading produced by the plate convergence, the occurrences of SSEs tend to decrease the
interseismic coupling estimation. In subduction regions with important transient slip, it is
crucial to consider the observational period over which the coupling is estimated.

The possibility to estimate interseismic coupling along long sections of subduction zones
is then important to better understand its variability despite the relatively limited obser-
vation periods allowed by geodetic observations. In this study, we target a 1000-km-long
segment of the Mexican subduction zone (Figure 1) prone to large earthquakes and SSEs,
using a combination of InSAR and GNSS data. Along the Mexican subduction zone, the
Rivera and Cocos plates converge towards the North American plate with a northeast dir-
ection (Figure 1). The convergence velocity increases eastward from 4 cm/year in Jalisco to
7 cm/yr in Chiapas (DeMets et al., 2010). The subduction interface geometry varies along
strike with a flat slab at 40 km depth between -101°E and -96°E, which extends about 200
km from the coast. Many historical earthquakes have been observed since the beginning
of the 20th century (e.g. Singh et al., 1981; Anderson et al., 1989). In the eastern part of
the subduction zone (in Oaxaca) the recurrence time of Mw ~ 7.0 events is about 40 years
(Nishenko & Singh, 1987). In the western part, much larger earthquakes with magnitudes
above 8 have been observed in Michoacan (1985) and Jalisco (1932, 1995) (Courboulex
et al., 1997; Anderson et al., 1986), with larger but not well constrained recurrence time
(Nishenko & Singh, 1987). Along this subduction, SSEs are known and monitored by geo-
detic measurements in the regions of Guerrero and Oaxaca (see Figure 1). Several large
periodic SSEs have been detected by continuous GNSS stations in Guerrero in 1998 (Lowry
et al., 2001), 2002 (Kostoglodov et al., 2003), 2006 (Radiguet et al., 2011; Kostoglodov et
al., 2010), 2009-2010 (Walpersdorf et al., 2011), 2014 (Radiguet et al., 2016), and 2017-2018
(Cruz-Atienza et al., 2021), and/or using InSAR data (e.g. Cavalié et al., 2013; Maubant
et al., 2020). These SSEs have equivalent magnitudes of Mw 7-7.5 and a recurrence time
around 4 years. Several periodic large SSEs have also been detected in the Oaxaca region.
These SSEs have equivalent magnitudes of $M_w$ 6.6-6.9 and a recurrence time of 1-2 years (Brudzinski et al., 2007; Graham et al., 2016).

Based on previous geodetic studies providing a slip distribution of the SSEs in this region (Cavalié et al., 2013; Radiguet et al., 2011), we know that these aseismic events are located between 15 to 40 km down-dip from the seismogenic zone in the vicinity of the Guerrero seismic gap (Kostoglodov et al., 2003; Radiguet et al., 2012). This seismic gap is a 200km-long segment, where no large earthquake has been observed during at least one century (Kostoglodov et al., 2003). It is suspected that SSEs penetrate into the Guerrero Gap and account for 75% of the strain release in this region, and thus that these SSEs play a role in the observed lateral variations of the coupling in this area (Radiguet et al., 2012). Tremors and Low Frequency Earthquakes (LFEs) have been detected downdip of the large SSE areas (Figure 1). Their occurrences are strongly correlated in time with the SSEs activity in the Guerrero and Oaxaca regions (Husker et al., 2019; Frank et al., 2015). Episodic burst of tremors are also observed in Jalisco and Michoacan, but they are sparse and occur in isolated clusters (Maury et al., 2018), and do not seem to be accompanied by SSEs (Brudzinski et al., 2016).

In 2017-2018 a sequence of earthquakes and SSEs occurred in Mexico (Figure 1). Two of these earthquakes were intra-slab: the first one was the $M_w$ 8.2 Chiapas earthquake, which occurred on September 7th 2017, the second was the $M_w$ 7.1 Puebla-Morelos earthquake on September 19th (e.g. Melgar et al., 2018; Segou & Parsons, 2018). The third event was a $M_w$ 7.2 interface earthquake which occurred on February 16th 2018, in Pinotepa. Afterwards, two other large interface earthquakes followed (Oaxaca, 2020 $M_w$ 7.4; Acapulco, 2021 $M_w$ 7.1). In the Guerrero region two SSEs have been observed in 2017-2018 and in 2019 (Maubant et al., 2020; Cruz-Atienza et al., 2021). In the Oaxaca region a SSE has been observed before the Pinotepa earthquake (Cruz-Atienza et al., 2021). All these seismic and aseismic events show the richness of the seismic cycle along the Mexican subduction zone on the time scale of geodetic observations. However, the spatial and temporal variations of
these phenomena also challenge interseismic coupling estimation, as shown by Villafuerte et al. (2021).

Several previous studies have analyzed the interseismic deformation using the continuous GNSS network in the regions of Oaxaca (Franco et al., 2005; Correa-Mora et al., 2008; Villafuerte et al., 2021) and Guerrero (e.g. Kostoglodov et al., 2003; Radiguet et al., 2012; Kazachkina et al., 2020) and a few studies have been carried out in the Michoacan and Jalisco regions (Azúa et al., 2002; Marquez-Azua & DeMets, 2009). Even though lateral variations of the plate coupling were observed with GNSS data in those studies, the network coverage was very sparse and the estimates are biased by a poor spatial constraint in the Jalisco and Michoacan states (Azúa et al., 2002). Up to now, no studies have been conducted spanning the whole Mexican subduction zone. Some studies combining InSAR and GNSS data to get interseismic velocity field over wide area have already been conducted in different tectonic contexts (e.g Weiss et al., 2020; Xu et al., 2021), but none on the Mexican subduction zone.

Since the start of the Sentinel-1 satellite SAR mission in 2014, which has a global and systematic regular acquisition strategy (down to 6-day revisit time), satellite SAR Interferometry (InSAR) can help to greatly improve the spatial resolution of geodetic measurements and provide a more homogeneous coverage at the scale of the subduction zone. In this study, we will focus on the period between 2016 and 2019.5, using the freely available Sentinel-1 data from the European Space Agency (ESA). We combine InSAR and GNSS data to obtain for the first time a homogeneous coverage of the interseismic coupling along a 1000 km long portion of the Mexican subduction zone (from Jalisco to Oaxaca) (Figure 1). Over such a large area along this subduction zone, the plate convergence velocity and obliquity, the subduction interface geometry, and the history and magnitude of seismic, postseismic and transient slow slip events (SSE) (Figure 1) significantly vary in space and time.

In the next section we will present the data used and the processing implemented to extract the linear term from the geodetic data. Then we present the estimation of the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Tectonic setting of the Mexican subduction zone and InSAR and GNSS data coverage.</figcaption></figure>

The black arrows indicate the convergence velocities of the Rivera (Riv.) and Cocos plates with respect to the North America Plate from the MORVEL 2010 model (DeMets et al., 2010). The black dashed-line with triangles indicates the Middle American Trench (MAT). Black contour lines represent the depth of the subduction slab between 0 and 70km (every 10km) (from Hayes et al., 2018). Dark blue lines are the 10 cm slip contours of SSE: to the west, the 2002, 2006, 2010 and 2017-2018 Guerrero SSEs (Radiguet et al., 2016, 2011; Maubant et al., 2020) and to the east the Oaxaca SSEs (Graham et al., 2016). The magenta lines are the tectonic tremor contours modified from Maury et al. (2018). The stars represent the epicenters of the large earthquakes sequence: Chiapas (07/09/2017), Puebla-Morelos (17/09/2017), Pinotepa (16/02/2018), Oaxaca (23/06/2020) and Acapulco (08/09/2021). Subduction earthquakes ruptures since the 20th century are represented in cyan. The dashed black rectangles represent the coverage of the four descending Sentinel-1 tracks processed in this study. The black dots represent the GNSS stations from different networks (UNAM, IGF, SSN, TLALOCNet).
linear trend extracted from the data, which allows retrieving by inversion the coupling
along the subduction. Then, we discuss the contribution of InSAR in a large scale study of
a subduction zone deformation. Finally, we compare the linear velocities calculated over a
short period of time (3.5 years) and longer term (5 to 16 years) velocities from GNSS data
in order to assess the representativeness of short term coupling over longer period of time.
This comparison allows us to discuss the spatio-temporal complexity of the coupling in the
Mexican subduction.

## 2 Data and Methods

### 2.1 InSAR data and processing

We use radar images from the Sentinel-1 European constellation operated by ESA.
Images are acquired in TOPS IW mode that is a large swath acquisition mode (~250 km
wide) with a return time of 6 to 12 days. We process four 250-400 km long descending
tracks, which together cover a ~1000 km long band of land from Jalisco to Oaxaca state
(from west to east, track relative orbit numbers 012, 114, 041 and 143, Figure 1) between
January 2016 and June 2019. Dates before 2016 were excluded because the time interval
between successive acquisitions was at least 24 days, which decreases the coherence and
makes the phase unwrapping difficult. We did not use ascending tracks due to a large gap
between May 2017 and September 2017, which makes it difficult to construct a reliable time
series and to estimate a robust trend. The InSAR data used in this study are presented in
table 1 and Figure S1. InSAR time-series are calculated using a SBAS strategy (Berardino
et al., 2002). We use the NSBAS processing chain, which partly relies on a set of modified
scripts and libraries originating from the legacy InSAR processing software Repeat Orbit
Interferometry PACkage (ROI\_PAC) (Doin et al., 2011; Rosen et al., 2004), expanded to
allow for the processing of a stack of interferograms and time series, and adapted to the
specificity of Sentinel-1 TOPS IW mode. We use a processing strategy similar to Maubant et
al. (2020), which is summarized below. To increase the signal-to-noise ratio, interferograms
<table><caption>**Table 1.** Table of the four Sentinel-1 descending tracks used in this study and the corresponding number of images and interferograms selected for each track (from West to East).</caption><thead><tr><th>Track name</th><th># images</th><th># interferograms</th></tr></thead><tbody><tr><td>D012</td><td>88</td><td>318</td></tr><tr><td>D114</td><td>78</td><td>319</td></tr><tr><td>D041</td><td>126</td><td>816</td></tr><tr><td>D143</td><td>114</td><td>681</td></tr></tbody></table>

are filtered using complex multi-looking on a window of 64 pixels in range × 16 pixels in azimuth, resulting in a spatial resolution of about 160 m × 240 m. We use empirical corrections for phase ramps and tropospheric signals to flatten the interferograms before unwrapping in order to minimize unwrapping errors. The final tropospheric correction will be done in the decomposition step described hereafter. After the unwrapping step, we add back the empirical corrections, as they might contain signal of interest. The phase of unwrapped interferograms is referenced with respect to a large area located 280km away from the Trench (Figure S4) far enough not to be significantly impacted by the subduction zone deformation. Then we invert the interferogram network with an SBAS approach to get phase time-series (Doin et al., 2015; López-Quiroz et al., 2009).

The InSAR time series are then decomposed with a parametric approach, with the goal to isolate the average velocity in line of sight (LOS) over the period of interest. For each pixel, the temporal evolution of the phase Φ(t) is decomposed into a linear term, a term corresponding to the evolution of the tropospheric signal (Atmospheric Phase Screen APS), and a co-seismic offsets when necessary, using the following equation:

$$ \phi_{\text{model}}(t) = V.t + b.APS(t) + \sum_{i=1}^{N} c_i.H(t - t_i) + d \quad (1) $$
where *t* is the acquisition time of each image of the time series for one track, *V* is the velocity, *b* is the amplitude of the tropospheric signal, *c*ᵢ the amplitude of the co-seismic offsets (defined with a Heaviside function *H*) of the event *i* occurring at the date *t*ᵢ, and *d* is a constant. The number *N* of considered co-seismic events can vary depending on the tracks, and here only the track D143 is corrected for the coseismic offsets of the September 17th 2017 Puebla and February 16th 2018 Pinotepa earthquakes. For each track, we impose the temporal evolution of the APS using tropospheric delays estimated by the GNSS stations present in the footprint of that track as in Maubant et al. (2020) where it has been shown that this approach is more efficient to correct for tropospheric perturbations in our study area than using atmospheric models like ERA-I. For each pixel, the parameters *V*, *b*, *c*ᵢ and *d* of the phase evolution model are estimated using a least-square adjustment.

Transient signals, such as SSEs or postseismic slip are present in our time series. It would be also interesting to separate the transient signals from the linear trend before the SSE (inter-SSE velocity), (e.g. Rousset et al., 2016; Radiguet et al., 2012). However, the duration of our InSAR time series and the number of dates available are not sufficient to separate them with confidence. Therefore, we do not correct our time series from SSEs or postseismic signals, and only estimate an average velocity over the study period.

For each InSAR track, we obtain a LOS velocity map, that is the ground velocity along the radar Line of Sight (LOS), positive toward satellite, corresponding to the linear term (*V* in eq. 1, Figure 2). We observe a north-south velocity gradient in the LOS displacement along the tracks (Figure 2). Several cities are known to be subsiding in central Mexico (Chaussard et al., 2014) and have been masked automatically in the InSAR data to show exclusively the tectonic signal caused by the subduction. We asked pixel in urban areas when they fall outside the range of two standard deviations of the entire velocity map (values below 8 cm/yr). The uncertainty of the LOS velocity for each pixel (Figure S4) is computed as the root mean square error between the model and the data using equation 2:
$$RMSE = \sqrt{\frac{\sum_{t=1}^{N} (\phi_{data}(t) - \phi_{model}(t))^2}{N}} \quad (2)$$

where N is the number of acquisition dates, and $\phi_{model}$ is the result of the parametrized decomposition.

## 2.2 GNSS data

In addition to InSAR data, we use GNSS displacement time series. In the region of interest, the data of 62 continuous GNSS stations are available during the 2016 to 2019.5 time period. We processed the GNSS data with the GipsyX software (Bertiger et al., 2020) and referenced them with respect to the North America plate (NA). To estimate the interseismic velocity of the GNSS stations, we follow a similar approach as for the InSAR data. For each station, the displacement in each component (East, North, Up) is modelled as:

$$disp_{GNSS}(t) = V.t + \sum_{i=1}^{N} c_i.H(t - t_i) + d \quad (3)$$

where t is the time, V is the velocity, $c_i$ the coseismic offset (defined by an Heaviside function) of the i th event occuring at the date $t_i$, with N the number of co-seismic events considered for each station when it is necessary, and d is a constant. The number N vary between 0 and 2. We consider the Puebla-Morelos and Oaxaca earthquakes, and not the Chiapas earthquake that we consider too far away. We corrected the GNSS stations present in the footprint of the same InSAR track where the post-seismic has been corrected (D143). To ensure consistency between the GNSS and InSAR velocities, we do not correct the data for contributions from transient aseismic displacements. Note also that no tropospheric term is considered for the GNSS model as corrections of tropospheric delays are already accounted for in the GNSS time series processing. The uncertainties associated with the velocity estimates are computed as for the InSAR data, using the residuals between the model from equation (3) and the GNSS time series.
Figure 2 (top) shows the GNSS velocity map with the horizontal and vertical velocities estimated from the time series. Most of the GNSS velocity vectors are pointing landward in the northeast direction, as expected in the interseismic period for a coupled subduction thrust; however in the area of large episodes of transient aseismic displacement (SSE or post-seismic), the velocity vectors are deviating from the general interseismic trend, and some of them are pointing to the southwest, toward the trench. This deviation corresponds to the dominant unloading process during the transient events.

## 2.3 GNSS - INSAR comparison

To check that our two datasets are in agreement we compare the velocity obtained from GNSS data projected in the Line-Of-Sight (LOS) to the velocity extracted from InSAR data. InSAR and GNSS results are in good agreement (Figure S3).

## 2.4 Inversion model

To model the displacement rates observed at the surface, we set up a static inversion to retrieve the slip deficit rates on the subduction interface, using a back-slip approach (Savage, 1983). Since the velocity field we modeled corresponds to a elastic strain accumulation and/or release of the upper plate motion during the three years considered, we do not constrain the sign of the slip inverted, positive slip deficit corresponds to back slip (i.e. plate coupling) and negative slip deficit to forward slip (during transients event). For simplicity, the results of the inversion are shown in the figures as a map of the coupling ratio, with negative coupling corresponding to forward slip. For the inversion, we use a linear least-squares algorithm (Tarantola, 2005) with the regularization scheme of Radiguet et al. (2011). We chose the Slab2.0 model for our 3D model geometry (Hayes et al., 2018). This model features significant lateral variations of the slab geometry along the Mexican subduction zone. In our forward model, the Green's functions are computed for a homogeneous elastic half-space using the analytical formulation of Okada (1992). To alleviate the inversion's computational cost, we reduce the number of InSAR velocity measurements in
each track by averaging pixel values with a 10 × 10 km² window. As our interseismic signal is characterized by a long spatial wavelength, we perform a uniform downsampling. The associated InSAR uncertainties are given from the map of RMSE (eq. 2 and Figure S4 using the same downsampling method). We neglect the covariance between pixels, and covariance terms between our GNSS and InSAR datasets in this inversion. The slip direction is fixed in the inversion, and is computed for each fault patch as the projection of the plate velocity vectors from the block model from DeMets et al. (2010) (see Text S1 and Figure S5).

We tested two different plate velocity models: PVEL and MORVEL 2010 (DeMets et al., 2010). Overall, the differences between the two models are not significant in term of slip direction or plate velocities, except for the Rivera plate velocity (between -107°E and -104°E), for which the plate velocities are smaller for the PVEL model than for the MORVEL 2010 model (Figure S5) (DeMets et al., 2010). Bandy et al. (1997) proposed that the velocity of the Rivera subduction could be around 4cm/yr, that is adequacy with the MORVEL 2010 model. Even if it is not clear yet which Rivera plate motion model is preferable, in this study we chose the MORVEL 2010 model as our plate velocity reference, to be consistent with Bandy et al. (1997).

The relative weight between the two datasets is introduced by a weighting factor in the covariance matrix of the data $C_d$ such that:

$$C_d = \begin{bmatrix} \alpha^2 C_{d\text{InSAR}} & 0 \\ 0 & (1 - \alpha)^2 C_{d\text{GNSS}} \end{bmatrix} \quad (4)$$

where $C_{d\text{GNSS}}$ is the covariance of the GNSS data, $C_{d\text{InSAR}}$ is the covariance of the InSAR data and $\alpha$ is the weighting coefficient. If $\alpha = 0$ the weight of the InSAR data is null, if $\alpha=1$ the weight of the GNSS data is null.

To check the influence of each data type, we inverted each dataset separately. We set first $\alpha=0.001$ and then $\alpha=0.999$, which effectively corresponds to the inversion of GNSS or InSAR data alone, respectively (Figure S6). The comparison of the inverted coupling maps
(Figure S6) shows that using only InSAR data, the coupling is lower in the Jalisco region with some negative values close to the trench, whereas the GNSS model shows negative values at greater depth (35-45 km). In this area where the GNSS network is sparse the resolution of the model using only GNSS is poor.

We explore different $\alpha$ values to evaluate the best compromise to get the lowest misfit $\chi_2$ for each dataset (Figure S7). We choose $\alpha = 0.4$, as it gives the lowest value of the cumulative misfit (sum of $\chi_2$ values for InSAR and GNSS see Text S2).

Once the optimal value of $\alpha$ is found ($\alpha=0.4$), we test the variability of our model to the regularization parameters through L curves. We perform an exploration of the damping value ($\sigma_{m0}$) and of the correlation length value $\lambda$. A large value of $\lambda$ allows only a large scale spatial variation of the coupling. We choose the best compromise between the misfit ($\chi_2$) and model roughness (evaluated with the $L_2$ norm of the model parameters) for these two parameters (Figure S8). We first select the optimal value of $\log_{10}(\sigma_{m0}) = -1.4$, and then search for the optimal $\lambda$ value (Figure S8) which we fix to $\lambda=120\text{km}$. The optimal model has a $\chi_2 = 0.44$.

## 3 Results and Discussion

To obtain a coupling map of the subduction zone of Mexico, we divide the back slip velocity map (Figure S9) by the velocities of the Rivera and Cocos plates (Figure S5). The coupling ranges between -0.8 and 0.8 (Figure 4), a null coupling value corresponds to slip at the plate convergence rate. Most of the shallow part of the subduction interface is coupled (0-40km depth), except in areas where large SSEs or post-seismic slip occur. From Jalisco to the western part of Guerrero, the coupling is high (0.8), which is consistent with the fact that no large SSE displacements have been observed so far in this area. This high coupling region coincides with previous earthquake slip contours (light blue contours in Figure 4), for example the one of the 1985 Michoacan earthquake. In Jalisco, despite uncertainties about the plate velocity, with significant differences between MORVEL and PVEL models,
our result shows a high coupling ratio for both MORVEL and PVEL models (0.5 and 1 respectively). This high coupling ratio is consistent with recurrence time of 60-80 yr or longer of large earthquakes in this region (Nishenko & Singh, 1987).

In the Guerrero and Oaxaca regions, there are negative values of coupling ratio due to the presence of SSEs and post-seismic motions during the analyzed period. However between the longitudes -99°E and -100°E, a zone of positive coupling can be observed, indicating that the two SSE regions are separated or at least do not propagate into this seismogenic interface. This is consistent with the fact that the Mw 7.1 September 8th 2021, Acapulco earthquake occured at the boundary of this coupled area. Similarly, in the eastern limit of the study area, the *M*<sub>w</sub> 7.4 June 23rd 2020 Oaxaca earthquake occured. There, the resolution of our coupling map is low due to a lack of GNSS and InSAR data, nevertheless, we can observe that before this Oaxaca event the coupling had been quite high at the epicenter of the earthquake (around 0.7) (Figure 4).

The present study allows to estimate the contribution of the Sentinel-1 InSAR data to improve the capability of assessing the interseismic coupling over long sections of subduction zones where the GNSS station coverage is not dense enough. We compared the modeled surface displacements using different weighting of InSAR and GNSS datasets with the observations. Comparing the surface displacements predicted by the coupling model obtained by inverting only GNSS data (α=0.001, Figure S6) and our optimally weighted model (α=0.4, Figure 4), we can see that in regions where the spatial coverage of GNSS is low, the GNSS-only model (α=0.001) does not reproduce InSAR observations (Figure S10). In particular, the GNSS-only model predicts an area of subsidence between -100°E and -98°E which is not observed in the InSAR data (blue pattern in the predicted LOS displacement map, Figure S10). This observation is due to the lack of GNSS data in the region which does not allow to fully constrain the model. The optimal joint model prediction (α=0.4) shows a much better agreement with InSAR data, and a reasonable fit to the GNSS data (Figure 5). This highlights the key role of InSAR data to stabilize the inversion
results in areas where GNSS coverage is low. Nonetheless, some significant residuals remain locally in our optimal model like on the western part of the D143 track (profile DD' Figure 5) where transient deformation signal from Guerrero SSE and postseismic of the Pinetopa Earthquake occured. There, the spatial smoothing chosen to our model is probably strong to reproduce such deformation signals which are at relatively short wavelength in space.

One key question related to this work is the reliability of our coupling map estimated in a short 3.5-year period for longer time scales, especially in a region where deformation transients are common. To find an answer to this question, we investigate the difference between long term and short term velocities from GNSS time series. We compare the "long-term" horizontal velocities estimated from 5 to 16 years time series from Kazachkina et al. (2020) with our "short-term" 3.5-year horizontal velocities. Figure 6 shows the ratio between short and long-term velocities $V_{st}/V_{lt}$ for GNSS stations sorted by longitude. It reveals that in the Jalisco and Michoacan regions (from -106°E to -102°E longitude) $V_{st}/V_{lt} \sim 1.0$ which means that short and long-term velocities are similar. This is due to the absence of large transients during both the short and long-term observation periods. Previous studies, based on seven years of GNSS observation from 1993 to 2001, estimated that in the Jalisco region the subduction was fully locked before the 1995 Colima-Jalisco earthquake and re-coupled after a post-seismic period of two years (Marquez-Azua & DeMets, 2009; Cosenza-Muralles et al., 2021a; Azúa et al., 2002). The short term coupling retrieved in our study is also in agreement with the long term coupling retrieved by Cosenza-Muralles et al. (2021b). Although this study is performed with a denser GNSS network, the average coupling recovered is between 0.8 and 1. We can conclude that in Jalisco and Michoacan, our 3.5-year coupling model is representative of the longer-term coupling.

In Guerrero and Oaxaca however, where aseismic transients occur during our 3.5-year observation period, the difference between short and long-term observations is substantial, and the negative GNSS velocity ratio (Figure 6) is a marker of forward slip during the
SSE episodes. In these two regions, previous studies revealed that the long-term coupling, averaged over several SSE cycles is low, due to the impact of SSEs (Radiguet et al., 2012; Rousset et al., 2016). However, our 3.5-year time period, which includes the SSEs and postseismic transient from the 2018 Pinotepa earthquake, is not long enough to provide a reliable estimate of the long-term coupling. One can note however that from -100°E to -99°E the short to long term velocity ratio is again close to one, which corresponds to the location of the highly coupled patch, corresponding to the rupture area of the $M_w$ 7.6, 1957, "Angel" earthquake (Singh et al., 1982), identified in Figure 4.

We conclude that outside the two areas with SSEs or recent post-seismic deformation, our interseismic coupling map over 3.5-years is representative of the coupling over longer period with the benefit of being homogeneously computed over a long section of the subduction zone. The map is based on both GNSS and InSAR data and thus is better constrained in the areas with a poor GNSS station coverage.

# 4 Conclusions and perspectives

Thanks to the spatial and temporal densification of the InSAR data with the Sentinel-1 constellation (6-12 days of intervals), it is now possible to create time series of high enough quality to extract reliable interseismic velocity estimates along the Mexican subduction zone within a relatively short period of time (3.5-years). We show that the contribution of InSAR improves the spatial resolution of the coupling model, especially in regions where the GNSS network is sparse.

Our model confirms a large lateral variation of coupling along the Mexican subduction zone. In both the Guerrero and Oaxaca areas, aseismic slip events occurred between 2016 and 2019, and post-seismic deformation occurred in the Oaxaca region. These events are responsible for the apparent negative coupling (i.e. reverse slip on the plate interface). Our results show that the rupture areas of the $M_w$7.4, 23 June 2020 (Oaxaca) and $M_w$7.1, 8 September 2021 (Acapulco) earthquakes were coupled before these events (coupling is ~ 60
366

367

368

369

370

371

372

373

374

375

- 70%, Figure 4). A study of stress variations in both areas would allow us to better understand what kind of interactions may have occurred between the weak and strong coupling regions. In the Michoacan area the results show a high degree of coupling, which means that this region is continuously accumulating elastic strain and could host a future large earthquake. The Jalisco region is also a strongly coupled area. However, large uncertainty exists on the plate velocity model in Jalisco and hence on the exact coupling values. A comparison between the GNSS interseismic velocities estimated within our three and a half year period and velocities estimated over 5-16 year period shows a good consistency in areas where no transient events occur. This gives us confidence in the reliability of our interseismic coupling model.

376

377

378

379

380

381

382

383

384

385

386

387

The present study confirms the great potential of Sentinel-1 InSAR time series to constrain the interseismic deformation in subduction zones. In a few years, longer InSAR time series will be available, and similar approaches could be used to investigate coupling variations not only spatially but also temporally. In this study, we could only exploit data from the descending orbit without being able to take advantage of the combination of ascending and descending tracks to better constrain the vertical component of displacements. A more complete Sentinel-1 dataset should make it possible to process tracks located further south-east and thus extend the zone of estimated coupling as far as the southernmost state of Chiapas in Mexico. In the future we will be able to produce a long-term coupling map including InSAR and GNSS data. Longer time series will also allow the observation of several SSE cycles which will make possible the investigation of both inter-SSE and interseismic coupling.

388

Acknowledgement

389

390

391

Sentinel-1 data are available online https://scihub.copernicus.eu, GNSS data were provided by the National Seismological Service of Mexico (SSN) and Department of Seismology, UNAM (PAPIIT G100921, and CONACYT 284212 grants), whose staff are acknow-
ledged for station maintenance and data acquisition. Plate velocity models are available on https://www.unavco.org/. Data are available at: https://doi.org/10.5281/zenodo.3666702. USGS models are available at: https://earthquake.usgs.gov/earthquakes/eventpage/us6000ah9t/executive and https://earthquake.usgs.gov/earthquakes/eventpage/us7000f93v/executive. All computation results presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is partly supported by the EquipMeso project (ANR-10-EQPX-29-01) of the PIA supervised by the French Agence Nationale pour la Recherche. This research was supported by the CNES (L. Maubant Ph.D. funding and Grant APR-CNES for projet SSEMEX), Labex OSUG@2020. We thank Simon Daout for discussion about InSAR processing. The Editors, Roland Bürgmann and an anonymous reviewers are gratefully acknowledged for their valuable comments.

## References

Anderson, J., Bodin, P., Brune, J., Prince, J., Singh, S., Quaas, R., & Onate, M. (1986). Strong ground motion from the michoacan, mexico, earthquake. *Science*, 233(4768), 1043–1049.

Anderson, J., Singh, S. K., Espindola, J., & Yamamoto, J. (1989). Seismic strain release in the mexican subduction thrust. *Physics of the earth and planetary interiors*, 58(4), 307–322.

Azúa, B. M., DeMets, C., & Masterlark, T. (2002). Strong interseismic coupling, fault afterslip, and viscoelastic flow before and after the oct. 9, 1995 colima-jalisco earthquake: Continuous gps measurements from colima, mexico. *Geophysical research letters*, 29(8), 122–1.

Bandy, W., Kostoglodov, V., Singh, S., Pardo, M., Pacheco, J., & Urrutia-Fucugauchi, J. (1997). Implications of the october 1995 colima-jalisco mexico earthquakes on the rivera—north america euler vector. *Geophysical research letters*, 24(4), 485–488.

Berardino, P., Fornaro, G., Lanari, R., & Sansosti, E. (2002). A new algorithm for surface
deformation monitoring based on small baseline differential SAR interferograms. *IEEE Transactions on Geoscience and Remote Sensing*. doi: 10.1109/TGRS.2002.803792

419

420

Bertiger, W., Bar-Sever, Y., Dorsey, A., Haines, B., Harvey, N., Hemberger, D., ... others (2020). Gipsyx/rtgx, a new tool set for space geodetic operations and research. *arXiv preprint arXiv:2004.13124*.

421

422

423

Brudzinski, M., Cabral-Cano, E., Correa-Mora, F., Demets, C., & Márquez-Azúa, B. (2007). Slow slip transients along the Oaxaca subduction segment from 1993 to 2007. *Geophysical Journal International*, 171(2), 523-538. doi: 10.1111/j.1365-246X.2007.03542.x

424

425

426

Brudzinski, M., Schlanser, K. M., Kelly, N. J., DeMets, C., Grand, S. P., Márquez-Azúa, B., & Cabral-Cano, E. (2016). Tectonic tremor and slow slip along the northwestern section of the Mexico subduction zone. *Earth and Planetary Science Letters*, 454, 259-271. Retrieved from http://dx.doi.org/10.1016/j.epsl.2016.08.004 doi: 10.1016/j.epsl.2016.08.004

427

428

429

430

431

Cavalié, O., Pathier, E., Radiguet, M., Vergnolle, M., Cotte, N., Walpersdorf, A., ... Cotton, F. (2013). Slow slip event in the Mexican subduction zone: Evidence of shallower slip in the Guerrero seismic gap for the 2006 event revealed by the joint inversion of InSAR and GPS data. *Earth and Planetary Science Letters*, 367, 52-60. doi: 10.1016/j.epsl.2013.02.020

432

433

434

435

436

Chaussard, E., Wdowinski, S., Cabral-Cano, E., & Amelung, F. (2014). Land subsidence in central Mexico detected by ALOS InSAR time-series. *Remote Sensing of Environment*, 140, 94-106. Retrieved from http://dx.doi.org/10.1016/j.rse.2013.08.038 doi: 10.1016/j.rse.2013.08.038

437

438

439

440

Correa-Mora, F., Demets, C., Cabral-Cano, E., Marquez-Azua, B., & Diaz-Molina, O. (2008). Interplate coupling and transient slip along the subduction interface beneath Oaxaca, Mexico. *Geophysical Journal International*, 175(1), 269-290. doi: 10.1111/j.1365-246X.2008.03910.x

441

442

443

444

Cosenza-Muralles, B., DeMets, C., Márquez-Azúa, B., Sánchez, O., Stock, J., Cabral-Cano, E., & McCaffrey, R. (2021a). Coseismic and postseismic deformation for the 1995 colima-

445

446
447

jalisco and 2003 tecomán thrust earthquakes, mexico subduction zone, from modelling of
gps data. *Geophysical Journal International*.

448

449

Cosenza-Muralles, B., DeMets, C., Márquez-Azúa, B., Sánchez, O., Stock, J., Cabral-Cano, E., & McCaffrey, R. (2021b). Gps-derived interseismic fault locking along the jalisco-colima segment of the mexico subduction zone. *Geophysical Journal International*.

450

451

452

Courboulex, F., Singh, S. K., Pacheco, J. F., & Ammon, C. J. (1997). The 1995 colima-jalisco, mexico, earthquake (mw 8): A study of the rupture process. *Geophysical Research Letters, 24*(9), 1019-1022.

453

454

455

Cruz-Atienza, V. M., Tago, J., Villafuerte, C., Wei, M., Garza-Girón, R., Dominguez, L. A., ... others (2021). Short-term interaction between silent and devastating earthquakes in mexico. *Nature communications*, 12(1), 1-14.

456

457

458

DeMets, C., Gordon, R. G., & Argus, D. F. (2010). Geologically current plate motions. *Geophysical Journal International*, 181(1), 1-80. doi: 10.1111/j.1365-246X.2009.04491.x

459

460

Doin, M.-P., Lasserre, C., & Grandin, R. (2015). InSAR Processing Sentinel 1 data Case study of subsidence in Mexico city.
doi: 10.5270/Fringe2015.pp116

461

462

463

Doin, M.-P., Lodge, F., Guillaso, S., Jolivet, R., Lasserre, C., Ducret, G., ... Pinel, V. (2011). Presentation of the small baseline NSBAS processing chain on a case example: the Etna deformation monitoring from 2003 to 2010 using Envisat data. *Proceedings of the ESA 'Fringe 2011 Workshop', Frascati, Italy, (19-23 September 2011), 2011*(January), 19-23.

464

465

466

467

468

Franco, S. I., Kostoglakov, V., Larson, K. M., Manea, V. C., Manea, M., & Santiago, J. A. (2005). Propagation of the 2001-2002 silent earthquake and interplate coupling in the Oaxaca subduction zone, Mexico. *Earth, Planets and Space*, 57(10), 973-985. doi: 10.1186/BF03351876

469

470

471

472

Frank, W. B., Radiguet, M., Rousset, B., Shapiro, N. M., Husker, A. L., Kostoglakov, V., ... Campillo, M. (2015). Uncovering the geodetic signature of silent slip through repeating earthquakes. *Geophysical Research Letters*, 42(8), 2774-2779. doi: 10.1002/

473

474
2015GL063685

Graham, S., DeMets, C., Cabral-Cano, E., Kostoglodov, V., Rousset, B., Walpersdorf, A., ... Salazar-Tlaczani, L. (2016). Slow Slip History for the MEXICO Subduction Zone: 2005 Through 2011. *Pure and Applied Geophysics*, **173**(10-11), 3445-3465. doi: 10.1007/s00024-015-1211-x

Hayes, G. P., Moore, G. L., Portner, D. E., Hearne, M., Flamme, H., Furtney, M., & Smoczyk, G. M. (2018). Slab2, a comprehensive subduction zone geometry model. *Science*, **362**(6410), 58-61.

Husker, A., Frank, W. B., Gonzalez, G., Avila, L., Kostoglodov, V., & Kazachkina, E. (2019, jan). Characteristic Tectonic Tremor Activity Observed Over Multiple Slow Slip Cycles in the Mexican Subduction Zone. *Journal of Geophysical Research: Solid Earth*, **124**(1), 599-608. doi: 10.1029/2018JB016517

Kazachkina, E., Kostoglodov, V., Cotte, N., Walpersdorf, A., RAMIREZ-HERRERA, M.-T., Gaidzik, K., ... Santiago, J. A. (2020). Active 650-km long fault system and xolapa sliver in southern mexico. *Frontiers in Earth Science*, **8**, 155.

Kostoglodov, V., Husker, A., Shapiro, N. M., Payero, J. S., Campillo, M., Cotte, N., & Clayton, R. (2010, dec). The 2006 slow slip event and nonvolcanic tremor in the Mexican subduction zone. *Geophysical Research Letters*, **37**(24), n/a-n/a. Retrieved from http://doi.wiley.com/10.1029/2010GL045424 doi: 10.1029/2010GL045424

Kostoglodov, V., Singh, S. K., Santiago, J. A., Franco, S. I., Larson, K. M., Lowry, A. R., & Bilham, R. (2003). A large silent earthquake in the Guerrero seismic gap, Mexico. *Geophysical Research Letters*, **30**(15), 513-530. doi: 10.1029/2003GL017219

López-Quiroz, P., Doin, M.-p., Tupin, F., Briole, P., & Nicolas, J. M. (2009). Time series analysis of Mexico City subsidence constrained by radar interferometry. *Journal of Applied Geophysics*, **69**(1), 1-15. Retrieved from http://dx.doi.org/10.1016/j.jappgeo.2009.02.006 doi: 10.1016/j.jappgeo.2009.02.006

Lowry, A. R., Larson, K. M., Kostoglodov, V., & Bilham, R. (2001). Transient fault slip in guerrero, southern mexico. *Geophysical Research Letters*, **28**(19), 3753-3756.
503

Marquez-Azua, B., & DeMets, C. (2009). Deformation of mexico from continuous gps from 1993 to 2008. *Geochemistry, Geophysics, Geosystems, 10*(2).

504

505

Maubant, L., Pathier, E., Daout, S., Radiguet, M., Doin, M.-P., Kazachkina, E., ... Walpersdorf, A. (2020). Independent component analysis and parametric approach for source separation in insar time series at regional scale: Application to the 2017–2018 slow slip event in guerrero (mexico). *Journal of Geophysical Research: Solid Earth*, 125(3), e2019JB018187.

506

507

508

509

510

Maury, J., Ide, S., Cruz-Atienza, V. M., & Kostoglodov, V. (2018). Spatiotemporal Variations in Slow Earthquakes Along the Mexican Subduction Zone. *Journal of Geophysical Research: Solid Earth*, 123(2), 1559–1575. doi: 10.1002/2017JB014690

511

512

513

Melgar, D., Pérez-Campos, X., Ramirez-Guzman, L., Spica, Z., Espíndola, V. H., Hammond, W. C., & Cabral-Cano, E. (2018). Bend Faulting at the Edge of a Flat Slab: The 2017 Mw7.1 Puebla-Morelos, Mexico Earthquake. *Geophysical Research Letters*, 45(6), 2633–2641. doi: 10.1002/2017GL076895

514

515

516

517

Nishenko, S., & Singh, S. (1987). The acapulco-ometepec, mexico, earthquakes of 1907-1982: Evidence for a variable recurrence history. *Bulletin of the Seismological Society of America*, 77(4), 1359–1367.

518

519

520

Nishimura, T., Hirasawa, T., Miyazaki, S., Sagiya, T., Tada, T., Miura, S., & Tanaka, K. (2004). Temporal change of interplate coupling in northeastern japan during 1995–2002 estimated from continuous gps observations. *Geophysical Journal International*, 157(2), 901–916.

521

522

523

524

Okada. (1992). Internal deformation due to shear and tensile faults in a half space. *Bulletin of the Seismological Society of America*, 82(2), 1018–1040. Retrieved from http://bssa.geoscienceworld.org/content/82/2/1018.short

525

526

527

Radiguet, M., Cotton, F., Vergnolle, M., Campillo, M., Valette, B., Kostoglodov, V., & Cotte, N. (2011). Spatial and temporal evolution of a long term slow slip event: The 2006 Guerrero Slow Slip Event. *Geophysical Journal International*, 184(2), 816–828. doi: 10.1111/j.1365-246X.2010.04866.x

528

529

530
531

Radiguet, M., Cotton, F., Vergnolle, M., Campillo, M., Walpersdorf, A., & Cotte, N. (2012).

532

Slow slip events and strain accumulation in the Guerrero gap, Mexico., *117*, 1-13. doi:

533

10.1029/2011JB008801

534

Radiguet, M., Perfettini, H., Cotte, N., Gualandi, A., Valette, B., Kostoglakov, V., ...

535

Campillo, M. (2016). Triggering of the 2014 Mw7.3 Papanoa earthquake by a slow slip

536

event in Guerrero, Mexico. *Nature Geoscience*, *9*(11), 829-833. doi: 10.1038/ngeo2817

537

Rosen, P. A., Hensley, S., Peltzer, G., & Simons, M. (2004). Updated Repeat Orbit

538

Interferometry Package Released. *Eos*, *85*, 47.

539

Rousset, B., Lasserre, C., Cubas, N., Graham, S., Radiguet, M., DeMets, C., ... others

540

(2016). Lateral variations of interplate coupling along the mexican subduction interface:

541

Relationships with long-term morphology and fault zone mechanical properties. *Pure and*

542

*Applied Geophysics*, *173*(10-11), 3467-3486.

543

Savage, J. C. (1983). A dislocation model of strain accumulation and release at a sub-

544

duction zone. *Journal of Geophysical Research*, *88*(B6), 4984-4996. doi: 10.1029/

545

JB088iB06p04984

546

Segou, M., & Parsons, T. (2018). Testing Earthquake Links in Mexico From 1978 to the

547

2017 M = 8.1 Chiapas and M = 7.1 Puebla Shocks. *Geophysical Research Letters*, *45*(2),

548

708-714. doi: 10.1002/2017GL076237

549

Singh, S., Astiz, L., & Havskov, J. (1981). Seismic gaps and recurrence periods of large

550

earthquakes along the mexican subduction zone: A reexamination. *Bulletin of the Seis-

551

*mological Society of America*, *71*(3), 827-843.

552

Singh, S., Espindola, J., Yamamoto, J., & Havskov, J. (1982). Seismic potential of acapulco-

553

san marcos region along the mexican subduction zone. *Geophysical Research Letters*, *9*(6),

554

633-636.

555

Tarantola, A. (2005). *Inverse problem theory and methods for model parameter estimation*.

556

Society for Industrial and Applied Mathematics.

557

Villafuerte, C., Cruz-Atienza, V., Tago, J., Solano, D., Garza-Girón, R., Franco, S., ...

558

Kostoglakov, V. (2021). Slow slip events and megathrust coupling changes reveal the
559

earthquake potential before the 2020 mw 7.4 huatulco, mexico, event. *Earth Space Sci.*
*Open Archive*, 56 (10.1002).

560

561

Wallace, L. M., Beavan, J., Bannister, S., & Williams, C. (2012). Simultaneous long-term and short-term slow slip events at the hikurangi subduction margin, new zealand: Implications for processes that control slow slip event occurrence, duration, and migration.

562

563

564

*Journal of Geophysical Research: Solid Earth*, 117 (B11).

565

Walpersdorf, A., Cotte, N., Kostoglodov, V., Vergnolle, M., Radiguet, M., Santiago, J. A., & Campillo, M. (2011). Two successive slow slip events evidenced in 2009-2010 by a dense GPS network in Guerrero, Mexico. *Geophysical Research Letters*, 38(15), 1-6. doi: 10.1029/2011GL048124

566

567

568

569

Weiss, J. R., Walters, R. J., Morishita, Y., Wright, T. J., Lazecky, M., Wang, H., ... others (2020). High-resolution surface velocities and strain for anatolia from sentinel-1 insar and gnss data. *Geophysical Research Letters*, 47(17), e2020GL087376.

570

571

572

Xu, X., Sandwell, D. T., Klein, E., & Bock, Y. (2021). Integrated sentinel-1 insar and gnss time-series along the san andreas fault system. *Journal of Geophysical Research: Solid Earth*, 126(11), e2021JB022579.

573

574
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Velocity maps for the period 2016-2019 along the Mexican subduction zone. Top: Horizontal (black arrows) and vertical (colored dots) velocities estimated from GNSS time series. Bottom: Velocities estimates in LOS for the four InSAR tracks. Black dots: GNSS stations, Magenta dots: GNSS stations used in figures 3 and 6. Black contour: 1 cm slip contour of the Guerrero (left) and Oaxaca (right) SSE between 2017 and 2019 (Cruz-Atienza et al., 2021). Degree indicate the incidence angle near range to far range.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Comparison between InSAR and GNSS in LOS (cm) for two GNSS stations, TNCM and TCOL, labeled in Figure 2. The black dots represent the GNSS data in LOS (cm), red dots are the InSAR data corrected from APS, for the same location as the GNSS station (averaged over 5x5 pixels). As the InSAR data are more scattered than the GNSS data due to residual atmospheric signal to help the comparison with GNSS, we smooth the InSAR time series using a sliding window of 3 months (red line). The blue and the green lines represent the linear trend estimated from InSAR and GNSS data respectively. Top: station TNCM in the Jalisco region where no transient deformation has been observed between 2016-2019. The linear term of the model (eq. 1) represents well the station motion during the study period. Bottom: TCOL station located in Guerrero region, the SSE is clearly visible in the GNSS data in 2018. The modeled linear term represents the average velocity over the period of this study on the raw InSAR time series.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Coupling map of the Mexican subduction zone using MORVEL 2010 model. Blue contours: Guerrero's SSEs (left) and Oaxaca's SSEs (right). Light blue contours: historical ruptures. Green and Orange stars are Puebla-Morelos and Pinotepa earthquake respectively. Yellow and red contours: slip at 10 cm of the 23/06/2020 Oaxaca and September 9th 2021 Acapulco earthquakes (from USGS).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison between observation and predicted surface displacement from our optimal model. Black arrows are the horizontal GNSS velocities, colored dots the vertical velocities from GNSS stations. Right column shows four profiles with a width of 0.2°. Grey dots represents the GNSS data, red dots are synthetic GNSS. Black lines represent the InSAR data, grey contour is the corresponding error, red lines are the synthetic InSAR.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Ratio between the short term and the long term velocities (horizontal components) from GNSS time series as a function of the longitude. $V_{st}$ are the short term velocities corresponding to the 3.5-year time periods considered in this study, and $V_{lt}$ are the long-term velocities from Kazachkina et al. (2020), estimated on periods from 5 to 16 years. Red and blue rectangles correspond to the regions of the Guerrero and Oaxaca SSEs respectively.</figcaption></figure>