Reassessing the 2006 Guerrero slow‐slip event, Mexico: Implications for large earthquakes in the Guerrero Gap

In Guerrero, Mexico, slow‐slip events have been observed in a seismic gap, where no earthquakes have occurred since 1911. A rupture of the entire gap today could result in a Mw 8.2–8.4 earthquake. However, it remains unclear how slow‐slip events change the stress field in the Guerrero seismic region and what their implications are for devastating earthquakes. Most earlier studies have relied on a sparse network of Global Navigation Satellite Systems measurements. Here we show that interferometric synthetic aperture radar can be used to improve the spatial resolution. We find that slip due to the 2006 slow‐slip event enters the seismogenic zone and the Guerrero Gap, with ∼5 cm slip reaching depths as shallow as 12 km. We show that slow slip is correlated with a highly coupled region and estimate that slow‐slip events have decreased the total accumulated moment since the end of the 2001/2002 slow‐slip event (4.7 years) by ∼50%. Nevertheless, even accounting for slow slip, the moment deficit in the Guerrero Gap increases each year by Mw∼6.8. The Guerrero Gap therefore still has the potential for a large earthquake, with a slip deficit equivalent to Mw∼8.15 accumulated over the last century. Correlation between the slow‐slip region and nonvolcanic tremor, and between slow slip and an ultraslow velocity layer, supports the hypothesis of a common source potentially related to high pore pressures.


Introduction
In southern Mexico, the Cocos Plate subducts beneath the North America plate at a rate of ∼6.1 cm/yr [DeMets et al., 2010], posing a huge hazard for coastal cities like Acapulco, and also for nearby Mexico City. A M w 8.0 earthquake in 1985 caused more than 10,000 deaths in Mexico City and 3-4 billion U.S. dollars of damage (U.S. Geological Survey). Since then, the city population has swelled to 21 million. Of particular concern is a region of the subduction interface where no significant earthquake has occurred since 1911 (101.2 • W-99.2 • W), referred to as the "Guerrero Gap" (Figure 1). The largest slip deficit is assumed to be located on the northwestern section of the gap (101.2 • W-100.4 • W) [Kostoglodov et al., 1996], as some slip is thought to have been released on the southeastern section by the M w 7.7 earthquake in 1957 and the M w 7.1 and 7.0 earthquakes in 1962 [Ortiz et al., 2000]. Prior to the discovery of slow-slip events by Lowry et al. [2001], Singh and Mortera [1991] estimated that a rupture of the Guerrero Gap today could result in an earthquake of M w 8. 2-8.4.
Slow earthquakes, also denoted as slow-slip events (SSEs), release strain over days to years. While SSEs are not in themselves hazardous, they alter the stress field and can therefore have an impact on the timing of subsequent, potentially damaging earthquakes. Since 1995, one large event has been observed in Guerrero approximately every 3-5 years: in 1995 Vergnolle et al., 2010], 1998 [Lowry et al., 2001;Larson et al., 2004;Vergnolle et al., 2010/2002[Kostoglodov et al., 2003Iglesias et al., 2004;Vergnolle et al., 2010], 2006 [Larson et al., 2007;Correa-Mora et al., 2009;Vergnolle et al., 2010;Radiguet et al., 2011;Bekaert, 2011;Hooper et al., 2012;Cavalie et al., 2013], and 2009/2010 [Walpersdorf et al., 2011]. Of the more recent events, the 2006 SSE started around April 2006 and lasted till December 2006 [Larson et al., 2007;Correa-Mora et al., 2009] or perhaps February 2007 [Vergnolle et al., 2010]. Continuous Global Navigation Satellite Systems (GNSS) observations by Vergnolle et al. [2010] (Figure 1) show that the surface displacements associated with this SSE are small (maximum ∼5.6 cm) and extend over a large area. Sites near the coast (CAYA and COYU) uplifted by ∼5.6 cm, over a period of 9-10 months. About 110 km inland, sites subsided up to ∼2.6 cm (MESC), with displacements due to the SSE decreasing to zero beyond 250-300 km. The cumulative horizontal displacements for most of the sites are in a southwesterly direction, with the largest displacement reaching approximately 5.5 cm at CAYA. In general, the surface displacements associated with the slow-slip events extend from the coast to Mexico City, 300 km inland. We focused the data using the ROI_PAC software [Rosen et al., 2004], followed by the interferometric processing using the DORIS software [Kampes et al., 2003]. In our processing we used the 90 m Shuttle Radar Topography Mission digital elevation model (SRTM DEM) [Farr et al., 2007] to reduce topographic contributions. The change in surface properties between two SAR acquisitions as well as the specific satellite acquisition geometry adds decorrelation noise [Zebker, 1992]. This makes unwrapping more difficult, particularly in the more vegetated regions near the coast. We reduce the noise contribution by applying the StaMPS time series InSAR technique [Hooper et al., 2007], in order to select only those pixels that demonstrate stable phase characteristics in time. A summary of the applied DORIS and StaMPS processing parameters are given in Table 2. Spatial-temporal variations in the troposphere leads to an additional signal in the interferograms, with a turbulent component (mostly affecting length scales up to a few kilometers) and a long-wavelength component (>10 km). The turbulent component results from troposphere dynamics, while the long-wavelength component is related to variations in pressure, temperature, and relative humidity. Tropospheric delays in interferograms record the integrated delay along the entire traveled path of the radar signal; hence, the long-wavelength component exhibits a correlation with the topography. In our study, tropospheric delay masks the smaller slow-slip displacement signal [Hooper et al., 2012]; we observe large tropospheric delays with an average correlation between the unwrapped phase and the topography of K Δ = 2.8 rad/km (∼1.2 cm/km) [Bekaert et al., 2015]. K Δ varies spatially and using a single value over a large region fails to improve the result. Instead, we account for the spatially varying tropospheric signal by applying a power law correction technique, assuming an exponent of 1.4 and a reference height of 7 km, both estimated from sounding data for the observed period of time. An in-depth discussion on the tropospheric correction technique and the results used in this study are contained in our companion paper [Bekaert et al., 2015]. The wrapped and unwrapped interferograms after time series processing and correction for tropospheric errors are included in Bekaert et al. [2015].   [Kampes et al., 2003] and StaMPS (S) [Hooper et al., 2012].

Time Series InSAR Analysis
Once we have corrected the data for the height-dependent tropospheric signal, we model the slow-slip deformation using elastic dislocation theory. As none of the SAR acquisitions were acquired during the event itself, the time series shows a discrete jump at the time of the SSE corresponding to the slow-slip displacement signal d SSE . We assume the secular velocity v sec to be constant, as was found in the GNSS processing [Vergnolle et al., 2010]. In our time series model, d SSE represents the cumulative slow-slip deformation as the deviation from the secular or inter-SSE velocity, which is consistent with the GNSS estimates from Vergnolle et al. [2010] used in this study. In our time series inversion we estimate d SSE and v sec simultaneously with DEM errors, which correlate with perpendicular baseline B perp , and an atmospheric and noise signal common to all interferograms, referred to as the residual master atmospheric signal, a m . The linear model describing the change of the interferometric phase Δ in time, for a single time series point, is where A is the design matrix and x is a vector of unknowns. E{⋅} represents the expectation operator, Δt the time between the master and slave acquisition, n the total number of observations in time (interferograms), and −4 ∕ the conversion from the interferometric phase to line-of-sight displacement. All interferograms up to time t k−1 are prior to the slip event and thus do not directly include an estimate for d SSE , hence the partial zero column in A. All interferograms help in constraining the secular velocity. We estimate the unknown parameters by applying weighted least squares, weighted by the inverse variance-covariance matrix of the data Q Δ . Note that v sec , d SSE , and c topo are independent of the choice of master date; only a m will change if a different master date is selected, under the assumption that the data are unwrapped correctly.
Correcting for the height-dependent tropospheric signal does not remove all the atmospheric signals. While the impact of the residual atmospheric signal is small on the estimation of the DEM errors and the master atmospheric signal, it strongly affects the estimation of the secular velocity and the slow-slip displacement signal. For these we observe large spatial variations near the coast and toward the north of our InSAR study area. We also observe a strong anticorrelation between them ( Figure 3). To account for this in our stochastic model, we define our diagonal variance-covariance matrix Q Δ such that we include the residual noise and slave atmospheric signal. For each interferogram we estimate the data variance, in a nondeforming area (shown by the blue polygon in Figure 1), after correction for an initial estimate of the DEM errors and the master atmospheric signal, estimated assuming uniform weights. The nondeforming area extends over the region where we do not expect line of sight (LOS) displacements associated with the interseismic deformation or the SSE [Vergnolle et al., 2010]. We also excluded Mexico City due to its large subsidence rate.  (Figures 3a-3d), a high spatial correlation is observed between the estimated displacements and the surface topography. After tropospheric correction, Figures 3e-3h, the signal-topography correlation is reduced, as can also be seen in a lower magnitude of the residual master atmosphere, Figure 3c versus Figure 3g. A, B, and C in Figure 3b give the location of the individual points for which the surface displacement time series is shown in Figure 5c.

Error Distribution of the Surface Displacements
We determine the error distribution of our model parameters, used for the slow-slip modeling on the subduction interface, through the percentile bootstrap method [Efron and Tibshirani, 1986]. As opposed to the regular algorithm, where n random samples are taken from n observations, we select n 1 and n 2 random samples respectively from the n 1 and n 2 observations (interferograms) before and after the SSE. Using this approach, the model design matrix remains invertible for each bootstrap simulation, i.e., there will always be at least one interferogram that spans the slip event. Each bootstrap simulation provides us with an estimate for the unknown model parameters (v sec , a m , d SSE , and c topo ) at each point.
Estimating the full variance-covariance matrix directly from bootstrapping is not possible, as outliers in the model strongly bias the covariance terms. Moreover, while the variance terms are not strongly biased by outliers, the variance estimate of a point P i obtained from the bootstrapped data contains, besides its actual variance, also a variance and covariance component of the reference area or point P j , 2 witĥ2 Δij the estimated variance at point P i computed from the bootstrapped data, 2 i and 2 j respectively the actual variance at P i and P j , and C ij (x) the covariance between both P i and P j , which varies with separation distance x. Assuming equal point variance P i and P j , equation (2) can be rewritten to the equation of a semivariogram [Wackernagel, 2003]; see Appendix A for a full derivation. We constructed our noise variance-covariance matrix using the covariance function fitted to the semivariogram of the slow-slip  [Kim et al., 2010] are shown by the green markers. Topography information as well as the GNSS station locations (black square markers) are plotted with respect to the distance to the coastal ACAP GNSS station, located in Acapulco. The interface is modeled in 25 × 25 km dislocation patches, extending 450 km in the strike direction. In total, we use 180 patches. Ocean bottom topography data have been extracted from the GEBCO_08 Grid (version 20091120, http://www.gebco.net), while for the elevations we made use of the SRTM DEM. displacement residuals . The latter is calculated as the difference between the bootstrapped estimates of the slow-slip displacement signal and the weighted least squares estimate, i.e., =d boot SSE −d SSE .

Modeling the 2006 SSE
The position and structure of the Guerrero subduction interface have been previously estimated from a variety of data, including seismicity and focal mechanisms [Singh and Pardo, 1993;Pardo and Suárez, 1995], geologic records, and volcanic features [Ferrari, 2004], and from the analysis of the Meso-America Subduction Experiment (MASE) broadband seismic network [Perez-Campos et al., 2008;Kim et al., 2010]. The Pacific Plate subducts beneath the North America plate about 66 km south-southwest off the coastal city Acapulco with a dip angle of 15 • and becomes subhorizontal at a depth of 40 km, approximately 80 km north-northeast from Acapulco ( Figure 4). The interface remains horizontal until about 230-250 km from the coast, after which it dips steeply below the Trans-Mexican Volcanic Belt, at an angle of 75 • [Kim et al., 2010;Chen and Clayton, 2012]. We use the same subduction geometry as Radiguet et al. [2012] and Cavalie et al. [2013] (red solid line in Figure 4), which correlates with observations from receiver functions [Kim et al., 2010] (green markers).
We model the slip event by dividing the fault into 25 × 25 km dislocation patches and inverting for slip on the interface to fit the slow-slip surface displacements. The relationship between slip on the interface and the surface displacement, assuming a homogeneous elastic half-space, is given by the Green's functions in Okada [1985]. We assume a Poisson's ratio of 0.25 and use the GNSS and our InSAR displacements during the SSE as our observations, resulting in the following slip model: LOS is a vector of the estimated InSAR displacements in the line of sight, d GNSS are vectors of the GNSS displacements (three components), G are matrices of Green's coefficients relating slip to surface displacement, and s is a vector of unknown slip values for each patch. A planar correction is included for the InSAR data ( pl x+ pl y+ pl ) to correct for the long-wavelength orbit errors, where x and y are the InSAR coordinates in a local reference frame. The inversion is weighted by the inverse of the variance-covariance matrix Q, as estimated in section 4. Due to the large number of observations and the long spatial wavelength of the slow-slip deformation signal, we resampled our initial InSAR observations to a 1 km grid before computing the slow-slip surface displacements using weighted least squares, taking care to propagate all errors.
The inverse problem is ill posed and highly sensitive to measurement noise. We therefore regularize the problem by introducing a prior constraint on smoothness, by including . The gradient between the slow-slip deformation signal and topography before and after tropospheric correction (right), estimated in 50 km along track segments using all data in range. Before correction, the signal is strongly correlated in the middle and north of the InSAR study region. After correction, this correlation is decreased, with on average a reduction of 1.1 cm/km. Light markers indicate the region around Mexico City, where subsidence contaminates the tropospheric correction estimates at the selected spatial frequency band of 2-8 km, which was chosen to avoid contamination of the long-wavelength (∼150 km) slow-slip signal. This region is neglected in the slow-slip modeling. Figure S1 in the supporting information gives the same comparison for the secular velocity. Line-of-sight surface displacement time series for individual locations, as shown in Figure 3, (b) before and (c) after the tropospheric correction is applied. Residual master atmosphere and DEM errors are removed. The displacement in time is modeled as an interseismic rate plus a slow-slip step (black solid line). Positive displacements refer to motion toward the satellite, comprising mostly of vertical uplift. Offsets on the y axis are chosen arbitrarily for clarity.
[2009], who select a smoothness coefficient based on a trade-off between the degrees of freedom in the model and the misfit [Correa-Mora et al., 2008]. We apply the method as proposed by Fukuda and Johnson [2008], which solves for smoothness in a Bayesian framework, allowing us to estimate the full probability distribution of each slipping patch, taking into consideration all possible smoothness values. We assume that the distribution of the model residuals d and the Laplacian of the slip ▽s have a normal distribution of zero mean, respectively, ∼ N (0, Q) and ∼ N ( 0, 2 I ) , with 2 an unknown scaling parameter of the Only significant slip, 95% confidence slip >0 cm, is shown. Inverting for the InSAR variance only assumes the estimated signal not to be contaminated from residual atmospheric signals. Inverting for the variance only allows slow slip to extend farther eastward than when using GNSS only. Inclusion of covariance downweights the InSAR toward the GNSS solution, potentially indicating that the InSAR signal is subjected to residual atmospheric signals. Slip extends within the seismogenic zone (10-25 km depths), black solid line, and enters the Guerrero Gap, dashed black line (black dotted line for the northwestern Guerrero Gap). Spatial correlation between the slow-slip region and nonvolcanic tremor  (blue dashed line) correlate with the existence of an ultraslow velocity layer [Song et al., 2009] (green circles). Absence of an ultraslow velocity layer is shown by the cyan circles [Song et al., 2009]. Rake (black arrows) is within ∼ 10 • of the updip direction. The second and third rows show the integrated seismic moment for the longitude at the trench and for depth. A peak at 40 km depth is introduced by the subhorizontal character of the interface at this depth.
Laplacian variance matrix. Unlike Fukuda and Johnson [2008], we assume that our data uncertainties are known, and we do not scale them. The posteriori probability distribution can be written as where K is a normalization constant, N the total number of observations (GNSS and InSAR combined), and M the number of dislocation patches. We treat 2 as unknown and sample its full probability distribution using a Markov chain Monte Carlo method [Mosegaard and Tarantola, 1995]. We do the same for the unknown InSAR plane coefficients ( pl , pl , and pl ) as well as the slip and rake of each slipping patch. We fix the slip rake so that it is within 20 • of the plate convergence direction. We sample our prior distributions using a variable step size which is evaluated every 10,000 model runs [Hooper et al., 2013].

InSAR Processing and InSAR Time Series Estimation
As no GNSS tropospheric estimates are available in Vergnolle et al. [2010], a direct comparison between GNSS and the estimated InSAR tropospheric delays is not possible. We are, however, able to compare the estimated tectonic components from InSAR and GNSS before and after tropospheric correction. Figure 3 gives the estimated time series components before (Figures 3a-3d) and after (Figures 3e-3h) tropospheric correction. Before tropospheric correction, spatial patterns correlating with the topography can be observed in the north and middle of our study regions. For the slow-slip displacements, the variation of the signal toward the north does not correlate with GNSS ( Figure 5a). After correction, the slow-slip signal is consistent with independent estimates from GNSS, and the gradient between the slow-slip displacement signal BEKAERT ET AL.
©2015. The Authors.  and the local topography decreases by 1.1 cm/km. We find that the root-mean-square error (RMSE) between the displacements and those predicted by our time series model decreases after tropospheric correction (Figures 5b and 5c). The reduction increases with distance to the InSAR reference area and has a maximum value of ∼1 cm [Bekaert et al., 2015].
For the secular velocity, we observe the largest displacement rates near the coast (Figure 3e). This signal results from coupling between the subducting slab and the continental plate, causing a buckling of the hanging wall block near the coast. As the horizontal motion is almost perpendicular to the radar line of sight, the signal reflects the vertical velocity almost exclusively. Further inland, the secular displacement rates decrease.
For the slow-slip displacement signal (Figure 3f ), we see the opposite behavior. This is a consequence of the continental plate (hanging wall) slipping back toward the trench, resulting in uplift near the trench and subsidence further inland. Further inland we see the signal level out and diminish. The estimated slow-slip deformation signal near Mexico City is a tropospheric correction artifact, introduced when deformation around Mexico City leaked into the 2-8 km spatial band filter.
For the estimated master atmospheric signal (Figure 3g), we find signals of up to 10 cm which do not show a clear correlation with topography. Toward the northwest of our study region, we observe ∼9 km wavelength mesoscale atmospheric features (wave pattern formation).
The DEM errors reveal a smoothly varying signal in along track with an across-track variation (Figure 3h). These signals are also present in the wrapped interferograms. We believe this signal to be a focusing artifact, BEKAERT ET AL.
©2015. The Authors. The result when inverting GNSS data alone, (middle row) when inverting GNSS together with full covariance information after tropospheric correction, and (bottom row) after tropospheric correction and when ignoring the InSAR covariance information. After tropospheric correction we find the InSAR data to be more consistent with the GNSS observations in the northern region. Near the coast the inclusion of covariance downweights the InSAR data, pushing it toward the GNSS solution. This is due to residual tropospheric atmosphere near the coastal region. Residuals reveal a misfit for the InSAR at the higher spatial resolution, imposed by the smoothness constrained on the slip. Displayed errors are 2 . All shown displacements are in the line of sight (LOS). Locked Guerrero Gap 6.9 (22.4%) 6.9 (18.8%) 6.9 (21.8%) (seismogenic zone) Locked northwest Guerrero Gap 6.7 (13.3%) 6.7 (10.7%) 6.7 (12.1%) (seismogenic zone) a A comparison is made between the moment release in the locked seimogenic zone (10-25 km depth) and the transition zone (25-42 km depth). Also, the moment magnitude release within the locked Guerrero Gap (101.2 • W-99.2 • W) and the northwest portion of the gap (101.2 • W-100.4 • W) are shown. Rigidity is assumed to be 3 × 10 10 Pa.
introduced by the zero-Doppler processing within ROI_PAC. We observed a similar feature in another study region, which disappeared with images focused by ESA. The errors introduced during the focusing stage are limited, as they are correlated with perpendicular baseline and therefore reliably estimated and subtracted.

Slow-Slip Modeling
In total, we performed 2.5 million forward models using Markov chain Monte Carlo sampling to sample the full parameter space for the slip and rake of each slipping patch, the smoothness, and the three InSAR plane coefficients. From all these inversions we kept the maximum a posteriori solution, which represents the maximum of the 364-D histogram (364 parameters) [Hooper et al., 2013]. When inverting GNSS only, this reduces by three, as no plane is included. Figures 6 and 7 give respectively our maximum a posteriori slip solution and the corresponding marginal probability density distributions for the largest slipping patch. The corresponding modeled observations, residuals, and profiles P1-P3 are shown in Figure 8. From the GNSS time series it is unclear if the Oaxaca GNSS sites are solely displaced due to the 2006 Oaxaca SSE or potentially a combination with the Guerrero 2006 slip event. While the displacements associated with the 2006 Guerrero event are expected to be small, we opted to exclude the Oaxaca GNSS sites (HUAT, PINO, and OAXA). Other studies instead force a zero displacement [e.g., Vergnolle et al., 2010;Radiguet et al., 2012;Cavalie et al., 2013].
We compared the solution found when inverting GNSS data only (Figure 6a) with those of a combined GNSS and InSAR inversion. We also compared solutions found when including and excluding InSAR covariance (Figures 6b and 6c).
For the GNSS-only solution, significant slip (slip lower bound of 95% confidence >0 cm) is mostly found within 100 km horizontal distance of the GNSS sites perpendicular to the trench. The slip release corresponds to an equivalent earthquake magnitude of M w 7.3, with the maximum slip (∼17 cm) concentrated at depths of 30-40 km. As summarized in Table 3, we find that most of the slip release (∼76%) occurs on the transition zone (25-40 km depths), with the remaining part almost completely (∼23%) in the seismogenic zone (10-25 km) defined from seismicity [Suárez et al., 1990;Singh and Pardo, 1993] and interseismic studies Larson et al., 2004]. The locked portion of the Guerrero Gap (101.2 • W-99.2 • W) itself accounts for ∼22% of the total release, equivalent to an M w 6.9 earthquake. For the northwest portion of the gap (101.2 • W-100.4 • W, no earthquake since 1911) this is ∼13%, equivalent to an M w 6.7 earthquake.
When including the full InSAR noise variance-covariance information (second column) the InSAR data are downweighted, pushing the solution more toward the GNSS solution on the western side of the InSAR track, while the InSAR still provides a low-level constraint. Including the InSAR data but neglecting covariance (third column) causes the slip solution to extend much more in southeasterly direction. While the maximum slip magnitude and location remain approximately unchanged, we find slip patches of ∼10 cm slip extending about 80 km more southeast than when using GNSS alone. As shown in Table 3, independent of which data are used, the maximum a posteriori slip solution has a similar estimated release of M w 7.3, with ∼19-22% being released in the locked Guerrero Gap (M w ∼6.9) and ∼11-13% in the northwest portion of the gap (M w ∼6.7). Inverting without covariance (but leaving variance) information tends to push the slip solution more downdip, with ∼1.3% more of the total moment released on the transition stage (25-42 km depths).
The maximum a posteriori solution tends toward the smoother side of the marginal probability distributions shown in Figure 7 and does not always coincide with the peaks of the histograms. This can be explained by the fact that the histogram peaks represent the maxima in a two-dimensional parameter space and neglect all other dimensions included in the maximum a posteriori solution. Each slip solution has a different smoothness value, which is also set by the Markov chain Monte Carlo sampling. Figure 9 gives the rougher and smoother slip solutions.

Slow-Slip Modeling Compared With Other Studies
Below, we compare the results when inverting for GNSS alone and the integration of GNSS and InSAR. For the latter we compare the inversions when using the full InSAR variance-covariance matrix, with those ignoring the InSAR covariance, as done by Cavalie et al. [2013]. In all our cases we find slip to enter the seismogenic zone (10-25 km depth) and the Guerrero Gap, reaching depths as shallow as 10 km.

GNSS Inversion
Our slow-slip solution derived using GNSS observations only, Figure 6a, compares well with the overall slip pattern of Correa-Mora et al. [2009], who found a similar spatial extent and magnitude of the slip distribution and a peak slip of ∼19 cm compared to our 17 cm. The solution of Radiguet et al. [2011], where smoothness is chosen based on a trade-off between roughness and misfit, is rougher, with a maximum slip just after the downdip transition on the subhorizontal section and a large updip asperity near the coastal GNSS sites ACAP and COYA. The asperity is likely introduced to fit the vertical displacements at COYO and ACAP, which differ from surrounding GNSS stations [Larson et al., 2007]. Our maximum a posteriori solution fits the GNSS observations well, with an RMSE of 0.3 cm and 0.3 cm for, respectively, the horizontal and vertical components. For the two coastal sites COYA and ACAP, we overestimate the vertical by an average of 0.1 cm, while fitting the horizontal displacements. In our rougher solution, shown in Figure 9, we are capable of estimating both components without significant residuals, illustrating the impact of smoothness. Compared to Vergnolle et al. [2010], we find a slightly smaller maximum slip (17 cm compared to 22.5 cm) and locate it 10 km farther downdip from the trench. Unlike Vergnolle et al. [2010], we find that significant slip extends to depths shallower than 19 km. Inverting with the same interface as Radiguet et al. [2012], we find that our peak slip coincides with their downdip extent of the higher slipping region (∼17 cm). We do not observe the large updip slip of ∼17 cm just north of the coastal ACAP and COYA GNSS sites, potentially introduced due to a rougher solution to decrease misfit of the observations. We find a smaller magnitude release for the total slip distribution than Radiguet et al. [2012], M w 7.3 compared to M w 7.5, as well as the release in the northwestern Guerrero Gap, M w 6.9 compared to M w 7.2. The difference in magnitude could be introduced by the difference in the assumed rigidity, which we assumed to be constant in this study, compared to a layered model in Radiguet et al. [2012].
Within our sampled marginal smoothness probability distribution, we find slip of ∼5 cm reaching depths as shallow as 12 km, which is larger and shallower than previous studies by Correa-Mora et al. [2009] and Vergnolle et al. [2010] but in overall agreement (∼1 cm larger) with the slip magnitudes by Radiguet et al. [2011Radiguet et al. [ , 2012. In agreement with previous studies by Radiguet et al. [2011Radiguet et al. [ , 2012, we find that slow slip occurs in the locked part of the Guerrero Gap.

GNSS and InSAR (With Covariance) Combined Inversion
Our joint GNSS and InSAR slip solution (Figure 6b), when including the InSAR data with full covariance information, tends toward the GNSS solution. The region of maximum slip remains west of the InSAR region toward the 40 km depth subhorizontal transition of the subduction interface. Slip enters the Guerrero Gap, with magnitudes of up to ∼5 cm reaching depths as shallow as 12 km. The InSAR data are not strongly weighted with the inclusion of the covariance information, as the slow-slip signal is correlated on similar length scales to atmospheric noise. Therefore, the Bayesian algorithm assigns a higher probability to keeping the slip smooth than fitting the InSAR data perfectly.

GNSS and InSAR (Neglecting Covariance) Combined Inversion
When ignoring the InSAR error covariance (Figure 6c), the long-wavelength signal is interpreted as slow slip, which requires more slip northward on the subhorizontal, 40 km depth, segment of the interface than the GNSS-only solution, with a maximum slip of ∼17 cm. This is in agreement with Cavalie et al. [2013], who explored the same InSAR track, with a different processing approach. However, unlike Cavalie et al. [2013], we find that the slip distribution is constrained to be parallel to the subducting trench and we do not find any slip extending in a northeast direction. Care needs to be taken with the eastern extent of the slow-slip region. Inverting without covariance assumes implicitly that the residual atmospheric signal is not correlated spatially. If there are correlated residual atmospheric signals, they are assumed to be part of the slow-slip displacements in the inversion. When combining GNSS and InSAR, we find a similar magnitude of slip as for GNSS (∼5 cm), reaching depths as shallow as 12 km. Cavalie et al. [2013] constructed eight interferometric pairs spanning the 2006 slow-slip event, multilooked to a ∼640 m pixel size in order to reduce noise and thus decorrelation. They used a model from GNSS measurements to correct for the interseismic displacement and corrected orbital effects and long-wavelength errors by removing a planar trend utilizing five to six GNSS stations, depending on the interferogram. In this study, we used more SAR acquisitions, in a time series approach, reducing the effects of noise and atmosphere and allowing us to simultaneously invert for the interseismic rates and slow-slip displacement signal. We include the uncertainty in orbit and other long-wavelength errors in our results, as is the case for the smoothness. Cavalie et al. [2013] did not use a dedicated troposphere correction as applied in this  et al., 2010] and with and inter-SSE coupling derived from GNSS by Radiguet et al. [2012]. After removal of the 2006 SSE we find the slip deficit (b) when inverting for GNSS only, (c) in combining GNSS and InSAR with full covariance information, and (d) when combining GNSS and InSAR ignoring covariance. In all cases the slow-slip release coincides with a high inter-SSE coupled region, reducing the slip deficit by more than 50% (Table 4). study; instead, they applied a simple stacking approach [Zebker et al., 1997] using the eight interferograms. In terms of modeling, we assumed the same subduction interface as Cavalie et al. [2013], but rather than imposing an arbitrary scaling factor between InSAR and GNSS, we let the variance of the observations control the weights. We use approximately the same GNSS observations, except that we exclude the OAXACA GNSS stations (PINO, HUAT, and OAXA) rather than forcing them to have zero displacement.

Correlation Between Slow Slip, Ultraslow Velocity Layer and Nonvolcanic Tremor
The coincidence of nonvolcanic tremor (NVT) and SSEs in Japan, Cascadia, known as episodic tremor and slip events, gives weight to the hypothesis that NVT and slow slip originate from the same source or that one initiates the other. We find that our slow-slip solution ( Figure 6) correlates well with the spatial occurrence of the NVT  and the presence of an ultraslow velocity layer [Song et al., 2009;Kim et al., 2010;Song and Kim, 2012]. As the MASE (Meso-American Subduction Experiment) network is a linear array, approximately perpendicular to the trench, the NVT sensitivity in the direction perpendicular to the array is limited. Increased energy release in NVT has been observed during the 2006 SSE [Rubinstein et al., 2010;Husker et al., 2012], with an updip trend in NVT locations as the 2006 SSE progresses . In space, we find that the NVT region occurs just downdip of the maximum slipping region, coinciding with areas slipping by 3-13 cm. The hypothesis of fluid flow and fluid processes at the plate interface generating seismic tremor [Rubinstein et al., 2010] is supported by the correlation between an ultraslow velocity layer [Perez-Campos et al., 2008;Song et al., 2009] and the NVT, low electric resistivity from magnetotelluric [Jodicke et al., 2006], and seismic velocity tomography . The correlation we observe between the slow-slip region and the ultraslow velocity layer supports the hypothesis of slow slip related to high pore pressures.

Earthquake Hazard
We evaluated the slip deficit introduced from the end of the 2001/2002 SSE till the end of the 2006 SSE, Figure 10. This interval of 4.7 years is fairly typical of the time between events, which ranges from BEKAERT ET AL.
3 to 5 years. To calculate the slip deficit at the Guerrero Gap in the absence of any SSE, we multiplied the MORVEL (Mid-Ocean Ridge VELocity) plate convergence rate of 6.1 cm/yr [DeMets et al., 2010] over this period (4.7 years) and subtracted the inter-SSE coupling solution by Radiguet et al. [2012], as derived from 12 years of GNSS observations. We then calculated the slip deficit by differencing the slow-slip estimate from this. As the rake only deviates by up to ∼ 10 • from the plate convergence direction at the Guerrero Gap, we do not account for the rake variations. Irrespective of inverting GNSS data alone, or when combining GNSS with InSAR, we find that slow slip reduces the total moment within the Guerrero Gap by ∼50% (Table 4). This leaves a residual slip deficit, build up over 4.7 years, equivalent to M w an 7.3 earthquake or M w ∼6.8 each year. In the northwestern section of the gap, where no earthquake has been observed since 1911, the slip deficit over the same time period is equivalent to M w 7-7.1. In agreement with Radiguet et al. [2012] we find the slip deficit outside the Guerrero Gap to be higher than inside the gap, where slow slip occurs. Nevertheless, considering the partial release of slow-slip events over time, the slip deficit within the Guerrero Gap can still cause large earthquakes and remains a threat for Mexico City. Assuming that all stress release happens through earthquakes and through slow slip, an equivalent M w ∼8.15 earthquake has been accumulated in the last century within the Guerrero Gap. Recently, two earthquakes occurred on 18 April 2014 (M w ∼7.2) and 8 May 2014 (M w ∼6.4) on the western edge of the Guerrero Gap. Even if they slipped entirely within the locked portion of the Guerrero Gap, they are too small to have a significant impact on the slip deficit.

Conclusions
Whether inverting GNSS observations alone, combining GNSS and InSAR while neglecting the covariance information, or while including the full InSAR covariance information, we find the 2006 Guerrero slow-slip event to have an equivalent earthquake magnitude of M w 7.3. This is in agreement with Correa-Mora et al. [2009] but smaller than previous studies, which estimated a M w ∼7.5. We find larger slip magnitudes (∼5 cm) in the seismogenic zone (10-25 km depths) in the Guerrero Gap than reported previously. Of the total slow-slip moment release we find that 19-25% is released in the locked part of the Guerrero Gap (equivalent to M w ∼6.9). The majority of slip, ∼75% of the total moment, is released in the transition zone (25-40 km depths). We allowed for rake variations within 20 • of the plate convergence direction but find most of the slip to be in the dip direction of the subducting interface, suggesting slip partitioning. The remnant strike-slip component is not released through slow slip. When inverting the combined observations from GNSS and InSAR, we find that the slow slip extends farther east, under the InSAR track. In contrast to a previous InSAR and GNSS study [Cavalie et al., 2013], we find that slip does not extend to the northeast but remains approximately parallel to the trench. Both nonvolcanic tremor and our slow-slip region correlate well with the presence of an ultraslow velocity layer, strengthening the hypothesis that both originate from the same process. We find that slow slip occurs within a highly inter-SSE coupled region. For the period between the end of the 2001/2002 and the end of the 2006 slow-slip event an M w 7.3 earthquake was accumulated or equivalent to M w 6.8 each year. This includes a reduction in total moment by ∼51% introduced by slow slip. In the northwestern portion of the gap, where no earthquake has been observed since 1911, the slip deficit is equivalent to a M w ∼7.1 earthquake. Assuming stress being released by earthquakes and slow slip, and considering an average 4 year slow-slip cycle and absence of a large earthquake in the last 100 years, the Guerrero Gap still has the potential for causing a M w ∼8 earthquake. The recent April and May 2014 Guerrero earthquakes that occurred on the western edge of the Guerrero Gap are too small to decrease this estimate significantly.

Appendix A: InSAR Covariance Matrix From Bootstrapping
We derive the InSAR variance-covariance matrix for the errors by applying the percentile bootstrap method [Efron and Tibshirani, 1986] and using point pair combinations of the errors to estimate a covariance function. As opposed to the regular bootstrap algorithm, where n random samples are taken from n observations, we select n 1 and n 2 random samples respectively from the n 1 and n 2 observations (interferograms) before and after the SSE, such that our model design matrix remains invertible for each bootstrap simulation, equation (1). We estimate the residuals of the slow-slip displacement, , between each bootstrap simulationd boot SSE and the best estimate from weighted least squaresd SSE , and use the residuals for the variance-covariance matrix estimation.
where k represents one of the n bootstrap runs. Estimating the covariance matrix directly from all the residuals is not possible as the covariance terms are biased by outliers. For the variance the influence of outliers is less. Using the bootstrapped residuals, the estimated variancê2 Δij between a point pair P ij follows aŝ2 which in theory corresponds to 2 Δij = 2 i + 2 j − 2C ij (A3) with 2 i and 2 j the variance at P i and P j and C ij covariance between both. By assuming the same variance for point P i and P j , equation (A3) becomes 2 Δij = 2 2 i∕j − 2C ij , which reduces to the definition of a semivariogram [Wackernagel, 2003] with C(0) the variance at point P i and P j and with representing the dissimilarity of the semivariogram depending on the separation distance x between the two points. By binning the experimental variogram and by fitting a Gaussian covariance function in combination with a nugget model for the uncorrelated terms, we were able to derive the full variance-covariance matrix for the slow-slip displacement. We validated the variance-covariance matrix for errors by using its inverse to weight a least squares estimation of the point centroid for each bin and ensuring that each estimated centroid lay within the cloud of points used in the estimation.

Erratum
In the originally published version of this article, the authors omitted an acknowledgment. The following sentence was added: "We thank the National Autonomous University of Mexico (UNAM) and all organizations involved for maintaining the local GPS network and making this data available to the community." This correction has been made to the online version of this article, but not the PDF, and the online version may be considered the authoritative version of record.