Adjoint Slip Inversion under a Constrained Optimization Framework: Revisiting the 2006 Guerrero Slow Slip Event

Understanding the fault behavior through geodetic data has an important impact in our assessment of the seismic hazard. To shed light on the aseismic evolution of a fault, we developed a new slip inversion strategy, the ELADIN (ELastostatic ADjoint INversion) method, that uses the adjoint elastostatic equations to eﬃciently compute the gradient of the cost function. ELADIN is a 2-steps inversion algorithm to better handle the slip constraints. In the ﬁrst step, it ﬁnds the slip that better explain the data without any constraints and the second step reﬁnes the solution imposing the slip constraints through a Gradient Projection Method. In order to get a physical plausible slip distribution and to overcome the poor fault illumination due to scarce data, ELADIN reduces the solution space by means of a von Karman autocorrelation function that controls the wavenumber content of the solution. To estimate the resolution, we propose a mobile checkerboard analysis which allows to measure a lower bound resolution over the fault for an expected slip patch size and an speciﬁc stations deployment. We test ELADIN with synthetic examples and use it to invert the 2006 Guerrero Slow Slip Event (SSE). The later is one of the most studied Mexican SSE that unfortunately was recorded with only 15 stations, so a strong regularization is required. We compared our slip solution with two published slip models and found that our solution preserves the general characteristics observed by the other models such as an updip penetration of the SSE in the Guerrero seismic Gap. Despite this similarity, our resolution analysis indicates that this updip aseismic slip penetration might not be a reliable feature of the 2006 SSE.


Introduction
An elegant and powerful mean to solve geophysical inverse problems is the adjoint method (AM ).Given an objective function, C, measuring the difference between data and a model prediction (i.e., a forward problem), to determine the model parameters that minimize C, the AM allows computing efficiently the derivative of C with respect to the parameters by combining the forward problem and the solution of an adjoint equation (i.e., of an adjoint problem) (Fichtner et al., 2006;Tromp et al., 2005;Tarantola, 1984;Gauthier et al., 1986).Thus, the inverse problem can be solved by using any optimization method that exploits that derivative to find the minimum of C. The AM provides a robust theoretical framework for constrained optimization problems that has been used in many geophysical inverse problems.The AM has been successfully used to solve full-waveform inverse problems in seismology, either to determine the elastic properties of the earth (Tromp et al., 2005;Askan et al., 2007;Fichtner et al., 2010;Krischer et al., 2018) or the kinematic history of earthquake sources (Sánchez-Reyes et al., 2018;Somala et al., 2018).For geodetic data, Kano et al. (2015) used the AM to estimate frictional parameters during the afterslip of an earthquake.
The longterm deformation of the Earth's crust close to the tectonic boundaries may be often explained in terms of the aseismic slip occurring between the plates.Depending on whether the interplate slip rate is larger or smaller than the relative plate motion, the plate interface experiences a relaxing slow slip event (SSE) (Dragert et al., 2001) or a stressing coupling regime (i.e., creeping or full locking) (Simpson et al., 1988), respectively.In the first case, a slip dislocation may predict the associated displacement field.In the later, the crustal deformation could be explained through the backslip formulation (Savage, 1983).The surface displacement is the summation of all contributions from the interface points experiencing either a coupling regime or an SSE.In the present work, to determine the plate interface aseismic slip history in these terms from continuous GPS (or any other geodetic) measurements, we introduce and solve a constrained optimization problem based on the adjoint elastostatic equations with a Tikhonov regularization term (Calvetti et al., 2000;Asnaashari et al., 2013) and a projection operator built with the von Karman autocorrelation function (Mai and Beroza, 2002;Amey et al., 2018).The new method, called ELADIN (ELastostatic ADjoint INversion), simultaneously determines the distribution of the interplate coupling and slow slip from surface displacements.
In those cases where the crustal strain field corresponds to a quasi-static seismotectonic process, the surface displacement is linearly related to the fault slip.However, determining the slip over an extended buried fault from such displacement remains an ill-posed problem.Underdetermination of the model parameters, i.e., the slip distribution, arises from the sparse sampling of the displacement field and the rapidly decreasing sensitivity of displacement to slip with distance from the fault (Nocquet, 2018).One rigorous framework to overcome this problem and to determine the uncertainty of such an inverse problem solution are the Bayesian approaches.The incorporation of prior information through probability density functions (pdf) allows determining the posterior model covariance and pdfs, as well as imposing model restrictions by means of truncated prior pdfs (Tarantola and Valette, 1982;Nocquet, 2018;Minson et al., 2013;Yabuki and Matsu'Ura, 1992;Amey et al., 2018;Nocquet et al., 2014;Nishimura et al., 2004).For instance, Minson et al. (2013) samples the posterior pdf using a Monte Carlo Markov Chain that enables to apply non-negativity constraints and any prior pdf.Dettmer et al. (2014) even showed how the optimal slip parametrization can be estimated from the data.Although Bayesian approaches are widely used and powerful, one important limitation that most have is the large computational load required to determine stochastically the posterior pdfs and thus the uncertainty of the model parameters.Yet, a proper selection of the prior pdf can overcome these issues making the sampling of the posterior pdf much faster through analytical or semi-analytical evaluations (Nocquet, 2018;Benavente et al., 2019).
An alternative to solve the elastostatic inverse problem is by introducing model regularizations and physically consistent restrictions.To prevent unrealistic oscillatory slip distributions, the most common regularization approach is to smooth the solution by applying a Laplacian operator (i.e., penalizing the second derivative of the slip) (McCaffrey et al., 2007;Wallace and Beavan, 2010).Usually, the hyperparameter controlling the smoothing is chosen subjectively by finding a satisfactory compromise between the data fit and the smoothing of the slip distribution.One common strategy to determine the hyperparameter is through an L-curve analysis that looks for an optimal value that keeps the data fitted with the strongest possible regularization (Radiguet et al., 2011).From a statistical approach, the hyperparameter can be determined using objective methods such as the Akaike Bayesian Information criterion (ABIC) (Yabuki and Matsu'Ura, 1992;Miyazaki et al., 2006) or fully Bayesian techniques (Fukuda and Johnson, 2008).Although the Laplacian operator reduces unphysical and rough slip solutions (and thus unreliable large stress drops), this is not the most convenient regularization strategy to preserve the real nature of the fault slip, which has a self-similar spectral signature (Mai and Beroza, 2002).Recently, Amey et al. (2018) proposed to use the von Karman autocorrelation function to build the model covariance matrix such that the penalization term should lead to self-similar slow-slip solutions.
When designing ELADIN, our goal was to introduce a regularization approach that would preserve the above-mentioned nature of faulting and, at the same time, allow a spectral control of the problem solution that guaranties a given resolution criterion.To this purpose we use a von Karman autocorrelation function that reduces the model space to a domain where the wavenumber content of all possible solutions satisfies a minimum slip characteristic length previously determined through robust resolution tests.In our approach, we do not build a model covariance matrix, as proposed by Amey et al. (2018), but we convolve the von Karman correlation function with the slip to project it into a reduced model space.We illustrate the capabilities of the method by inverting GPS data for the 2006 Guerrero SSE, which has been widely investigated in the literature, and describe several benefits that our solution has in comparison with some previous models.Systematic inversion of GPS data along the entire Mexican subduction zone applying the ELADIN method is presented in an associated work (Cruz-Atienza et al., 2020) where we analyzed the aseismic slip history of the plate interface between 2017 and 2019, and its interaction with large earthquakes.

The ELADIN Method
In this section, we first introduce the forward model that allows us to compute the synthetic displacements produced by a slip over the fault.Then, we formulate the inverse problem in a constrained optimization framework, reducing the solution space to control its spectral content with a von Karman autocorrelation function.We also include a Tikhonov term to penalize regions where slip is not expected to occur and to impose slip magnitude constraints.Finally, we present a 2-step algorithm that first solves the inverse problem without slip constraints using the adjoint equations for the gradient computation, and then projects the resulting solution into the feasible solution space which is later improved by means of a Gradient Projection method that imposes the desired physically-consistent slip constraints.

Forward model
The elastostatic representation theorem for the displacement field, u(x), due to a slip, d(ξ), produced at a fault, Σ, is where T i (•, •) is the i-component of the traction vector on the fault computed through the Somigliana tensor, S ij (ξ, x), and the fault normal vector n(ξ).It is a common practice to compute the Somigliana tensor considering a half-space homogeneous medium.However, Williams and Wallace (2015) have shown that neglecting the medium heterogeneities can result in slip overestimates of ∼ 20% for deeper events, and underestimates up to 42% for shallow earthquakes.For this reason, we adopted the AXITRA method (Bouchon and Aki, 1977;Coutant, 1990) for the calculation of the Somigiliana tensor, which allows us to consider heterogeneous layered media.
If the traction and the slip are projected along the plate convergence direction, c-, and the complementary perpendicular direction, p−direction, eq. ( 1) can be written in matrix form as Then, the fault is discretized in M subfaults such that the integral can be approximated as where A i is the i−subfault area.Finaly, if we want to compute the displacement for N receivers, we can order the displacements in a single vector such that the entire computation is reduced to a simple matrix-vector product as where

Inverse problem
The inverse problem consists in recovering the slip at each subfault of a known interface that produces displacements observed at geodetic stations.Due to the linearity of the forward model, eq. ( 4), we construct a quadratic cost function to formulate a convex inverse problem as where U o ∈ R 3N are the displacements observed at the N geodetic stations stored in a single ordered vector, as we did with U in eq. ( 4).Since real data are sparse and may have significant noise, the inverse problem ( 5) is ill-conditioned.In order to face these issues, a problem regularization and realistic physical constraints are introduced next.

Problem Regularization
Most often, the problem regularization is done by means of two elements: a model precision matrix and/or Tikhonov terms.
The model precision matrix is the inverse of the model covariance matrix which controls how sensitive is the slip in a given subfault to the slip on its neighbor subfaults.Radiguet et al. (2011) proposed a subfault correlation that follows a decreasing exponential function according to a defined correlation length.The problem we found with this approach is that the precision matrix for different correlation lengths does not have significantly different effects due to the fast decay of that function.For different types of correlation functions we tested, e.g., gausian and linear correlation functions, the model covariance matrix starts to become ill conditioned when the subfaults size becomes smaller than the correlation length.So, the precision matrices that could be computed were useless.
Tikhonov terms added to the cost function are used to penalize the roughness of the solution.Generally, the penalization is applied to the first or second spatial derivatives of the slip.However, when penalizing the derivatives, the norm of the slip solution is reduced as well.Besides, these two alternatives involve hyperparameters that need to be optimally determined because they control the tradeoff between the misfit of the data and the strength of the regularization.
These inconveniences lead us to propose a new approach that reduces the solution space so that the wavenumber content of the solution (i.e., the minimum characteristic length of the slip patches) can be controlled.The main idea is to apply a filter, projection operator, F , to the slip D.Then, the cost function ( 5) can be formulated as where C d is the data covariance matrix to weight the data according to their quality or proximity to the slip.
Recently, Amey et al. (2018) build a model covariance matrix with the von Karman autocorrelation function and showed that it is a good strategy to guarantee the slip self-similar properties (Mai and Beroza, 2002) that cannot be achieved with a common Laplace regularization.The spatial von Karman autocorrelation function is where H is the Hurst exponent, K H (•) is the modified Bessel function of second kind of order H, r is the correlation length that can be computed as where (s, d) are the coordinates in the along-strike and along-dip directions on the fault, and (a s , a d ) are the correlation lengths in those directions.Instead of using the von Karman autocorrelation function, eq. ( 7), to build a model covariance matrix, we propose to construct a linear operator K which, convolved with the slip D, controls the wavenumber content of the output function along both the strike and dip component.This convolution can be formulated as a matrix-vector product where the projection matrix, F , applies the convolution of the linear operator K to the slip, D, as was done in eq. ( 6) (see appendix A for further details).That is, F projects the slip, D, into a reduced solution space bounded by a chosen wavenumber.

Slip constraints
The model regularization we introduced guarantees that an optimal slip solution can be found.However, this solution may violate some expected physically-consistent restrictions, such as the full-coupling regime limit or slip rakes consistent with the plate convergence direction.Thus, slip constraints need to be imposed according to the available geological information.The cost function (6) can then be reformulated as subject to where β is a hyperparameter, W is a model-weight diagonal matrix that penalizes the slip per subfaults, D p is an a priori slip solution and (D j,l i , D j,u i ) are the lower and upper limits of the i-component of the slip in the j-regime.The slip is either in the SSE regime if its c-component is opposite to the plate convergence direction or in the coupling regime otherwise.If we have an a priori slip solution, D p , we can force the solution to be as close as possible to it by accepting only model changes that improve the data fit.In that case, the weight matrix should be the identity matrix, W = I.If no a priori slip information is available, we simply set D p = 0 and, to obtain the minimum norm solution, we make again W = I.Since we are not interested in getting the minimum norm solution in the present study, we thus set W = 0 everywhere except in the subfaults where we assume free slip (i.e., no coupling or SSE regime).The bigger the weighting value, the bigger is the subfaults slip penalization.The hyperparameter β controls the tradeoff between the fit of the data and the Tikhonov regularization term.Its value only guarantees that the solution does not contain significant slip in the penalized regions.On the other hand, if an a priori slip solution (D p = 0, W = I) is used or a minimal norm solution (D p = 0, W = I) is desired, then β must be determined following an optimal strategy as an L-curve analysis (e.g., Radiguet et al. (2011)) or the ABIC criterion (e.g., Miyazaki et al. (2006)).

Gradient computation: Adjoint method
To solve the inequality-constrained inverse problem (eqs.( 9)-( 11)), first we address the gradient of the cost function without considering the inequality constraints, eq. ( 11).In the framework of constrained inverse problems, the Lagrangian can be computed as where λ are the Lagrange multipliers.The Lagrangian total derivative with respect to the slip, D, is To simplify the computation of the gradient, we follow the adjoint method strategy (Fichtner et al., 2006).We start forcing ∇ λ L = 0 by solving a forward model Ũ = T F D. Then, we use the predicted displacement, Ũ , to compute the adjoint source As a result, the Lagrangian total derivative is the solution of the adjoint problem plus a term related with the slip constraints as Once the gradient of the cost function has been evaluated, we can follow any numerical optimization strategy to find the set of model parameters that minimize that function.

Gradient Projection Method
To avoid dealing with inequality constraints, it is often convenient to bring the current solution into the physically-consistent space after each iteration of the inversion procedure.However, for the slip inversion we realized that such procedure is not convenient because the gradient direction is often orthogonal to the slip constraints making the algorithm to stop.For large scale problems with lower and upper bounds for the variables, Nocedal and Wright (2006) propose the Gradient Projection Method (GPM) as an efficient strategy to deal with inequality restrictions.The GPM consists of two stages per iteration.In the first stage, the steepest descent direction is followed until a bound, i.e., the limit of an inequality constraint, is encountered and needs to be bent to stay feasible.Then, along the resulting piecewise-linear path, a local minimizer, called Cauchy point, is found (see Appendix B for details).For the second stage, a new optimal point is searched in the face of the feasible box on which the Cauchy point lies, i.e., those slip constraints that have reached a limit are changed to equality constraints.It implies that those inequality constraints are now part of the active set.This subproblem is usually not solved exactly since the remaining inequality constraints are usually not considered.
For the slip inversion, we do not follow exactly the GPM to avoid the subproblem of the second stage.This is because we expect that many subfaults in the coupling regime achieve its slip limit and that the number of iterations required was difficult to define.So, after computing the Cauchy point, we directly take it as a new iteration point where the gradient is computed again.Thus, our approach is essentially a steepest descent algorithm that respects the inequality constraints.Our GPM version is slow, so to achieve a fast convergence we then propose an algorithm that is explained in the next section.

2-step inversion algorithm
In order to increase the convergence speed, we developed a 2-step inversion algorithm.The purpose of the first step is to get an optimal initial solution for the GPM.In this step, we solve the unconstrained slip inverse problem using the adjoint method to compute the gradient of the cost function.Once the gradient is obtained, any iterative optimization algorithm can be used to find the optimal solution, e.g., the Conjugate Gradient method, the l-BFGS method, etc.In this work, we use the SEISCOPE optimization toolbox, which is a friendly and powerful optimization library developed in FORTRAN 90 with many available optimization strategies (Métivier and Brossier, 2016).After some performance trials, we decided to use the l-BFGS method.In the second step, we first project the solution into the physically-consistent domain and then we solve the constrained slip inverse problem with a slight modification of the GPM.As explained above, after computing the Cauchy point, instead of reformulating the inverse problem according to the new active set incorporating some inequality constraints, we use it as the new iteration of the slip.This is not a fast strategy, but since we start from a slip distribution that is close to the optimal solution, only a few iterations of the GMP are required (about 200).The pseudcode is described in Algorithm 1.
Step: Unconstrained slip inverse problem (Adjoint method) Data: GPS Data Initialize the slip D 0 = 0; while Convergence is not achieved do 1.Compute a forward problem

Compute the adjoint source
3. Compute the adjoint problem to get the gradient

Resolution
Resolution of our inverse problem essentially depends on the geometry configuration of the problem.This is, on the fault geometry and the distribution of observation sites, i.e., on the displacement field sampling and its sensitivity to dislocations at the fault.For a given problem discretization, synthetic inversions are a powerful mean to quantify how well an inverse method performs.If well-conceived, these tests may lead to very useful resolution information under realistic conditions, i.e., if they include data uncertainties and minimize the dependence on the target model.In the following, we present comprehensive exercises where the restitution of the target model is systematically quantified.To this purpose, for a given slip solution we define the restitution index, r i as where d T i and d I i are the slip of the target and inverted models at the i-subfault, respectively.The slip component used to determine the restitution index can be either the plate convergence or its perpendicular direction.We also introduce the average restitution index, ari, which is the mean of the restitution indexes over the M subfaults that discretize the 3D subduction interface between the Cocos and the North American plates in central Mexico (Cruz-Atienza et al., 2020).r i is one if the inverted slip equals the target slip and zero if the difference between them equals the target value.We have discretized the plate interface with subfaults whose surface projection is a square of 10 × 10 km 2 .To compute the static traction vectors along the interface due to single body forces at the stations, eq. ( 1), we assumed a four-layer 1D structure suitable for the region (Campillo et al., 1996) and used the AXITRA method (Bouchon and Aki, 1977;Coutant, 1990).For the analysis, we have considered all available permanent GPS stations (66 sites) in central Mexico (Cruz-Atienza et al., 2020) and 5 ocean bottom pressure gauges (OBP) deployed in the Guerrero seismic gap since November 2017 (Cruz-Atienza et al., 2018), where only the vertical displacements were considered.

Mobile checkerboard
A widely used strategy to quantify an inverse problem resolution is the checkerboard (CB) test.However, this test is intrinsically linked to the arbitrary choice of the target CB model, which means to the CB unit size, its positions and the absolute modelproperties periodically attributed.For this reason, we performed comprehensive mobile checkerboard (M-CB) tests for different patch sizes (PS).Based on previous GPS data inversions in central Mexico (Radiguet et al., 2012;Cruz-Atienza et al., 2020), we set up the checkerboards by alternating patches with 30 cm of slip (a typical value for SSEs in the region) and -10 cm of back slip (a cumulative value of slip deficit over 20 months in a fully coupled regime assuming a plate convergence rate of 6 cm/year).It should be noted that this type of checkerboard tests, where the resolution of the slip and coupling can be evaluated simultaneously, is not common practice because most of the available inversion methods can only handle these two quantities separately.
Figure 1 shows the inversion results for three CBs with different PS, [60, 80 and 100 km], and the same correlation length, L = 20 km.As we shall see, this value of L maximizes the average restitution index (ari) in these cases where no slip restrictions were imposed (i.e., no gradient projection method was used) and no data uncertainly was considered (i.e., the precision matrix was the identity matrix).Although the data fit is almost perfect in all three cases, it is clear that the target model restitution strongly depends on PS, the slip model characteristic length.As expected, the larger PS the better is the restitution.This is quantified in the right column, where the restitution index, r, is displayed for all subfaults.Besides, two more conclusions stand out: (1) restitution is better in SSE patches than in coupling patches, and (2) the inversion scheme cannot resolve the unrealistic slip discontinuity along the boundary of the CB patches.Both conclusions were expected because the backslip is one third of the positive slip, and because of both the imposed model regularization and the limited sensitivity of displacements with distance to the fault.
Previous results do not provide a reliable estimate of the problem resolution when facing real data because in that case we do not know the actual slip producing the observed displacements.A M-CB test consists in multiple CB inversions so that all possible model positions are explored.Results from the test may be translated into the mobile checkerboard restitution index (mcri) per subfault, which corresponds to the average of the r values estimated for each inversion.The mcri is a quantity that eliminates the resolution dependence on the CB configuration.For a given PS, we performed 6 M-CB tests, one without regularization, L = 0 km and the rest with different correlation lengths, L = 10, 20, 30, 40 and 50 km.Five different PS of 40, 60, 80, 100 and 120 km were considered so each case required a different number of CB inversions to complete the associated M-CB test.Since the horizontal projection of subfaults is 10 km per side and we shifted the CBs with a 20 km jump along the dip and strike directions to complete all possible configurations, the total number of CB inversions in a M-CB test for an given PS in km is (PS/10) 2 .
Figure 2 presents an overview of three M-CB tests for PS of 60, 80 and 100 km (those considered in Figure 1).As expected, in the top row, we see that the mcri increases with the PS, reaching values close to 0.8 in some regions close to the coast where there is the largest density of stations, and where the plate interface is closest to them.In deeper interface regions, between 30 and 50 km depth, mcri falls up to about 0.2 for PS of 60 km and over 0.5 for PS of 100 km along the whole subductions zone.As clearly seen in the right column of Figure 1, the unrealistic slip discontinuities along the patches edges makes the resitution very difficult, so we can considerer the mcri maps of Figure 2 (first row) as a lower resolution bound.Isocontours of these maps for different PSs and optimal correlation lengths thus define reliable fault regions where the inversions should resolve the unknown target slip above the mcri isocontour value (e.g., above 40% of the target slip if mcri equals 0.4).
The M-CB tests also allow to identify the optimum correlation length per subfault that maximizes the ari.This is shown in the second row of Figure 2, where we find that (1) the optimal L decreases for all PSs where the fault is well illuminated (i.e., in regions with high density of stations relative to the interface depth), and that (2) the optimal L increases as the PS decreases in places with sparse stations coverage, as it happens offshore near the subduction trench, for instance.Based on this multiscale analysis we built optimal solutions for the same CBs of Figure 1 by integrating the best inverted slip per subfault that corresponds to the associated optimal correlation length.Resolution improvements for the multiscale models ranged between 10% and 20% as shown in the third row of the figure (compare with the right column of Figure 1).However, something unexpected came out when comparing whole-interface average mcri values for all M-CB tests.Figure 3 shows this metric along with the average data-misfit error (i.e., the L2 norm of the difference between target and inverted displacements) for all tested PSs as a function of L, the correlation length.Although the spatial distribution of the optimal L depends on the slip characteristic length PS, the best average regularization was the same for all PSs and equal to 20 km.Such independency of the average mcri on L for different PSs is determined by the unrealistic slip jumps of the checkerboards slip values that sweep the whole interface no matter the PS.However, as we shall see latter, the optimal regularization length actually increases with PS if both the data uncertainty (the precision matrix) and the slip restrictions (the GPM step) are considered in the inversions.What is remarkable and was indeed expected in Figure 3 is that (1) the maximum restitution values increased with PS, (2) the restitution function for a given PS displayed a concave behavior and (3) the best fitting models are not the best solutions (i.e., those with the highest restitution).Regularization is thus critical to achieve physically acceptable and reliable slip models.

Gaussian slip
The analysis of the previous section did not consider the uncertainty in geodetic measurements that may be significantly large, especially in the vertical component where atmospheric noise and non-tectonic physical signals are particularly present.Nor did the analysis incorporate slip restrictions that are essential to guaranty tectonic expectations in our slip solutions such as backslip smaller than expected for a full-coupled interface regime and slip rake angles close to the plate convergence direction.For this reason, we now study three synthetic exercises in which the target slip corresponds to truncated Gaussian slip distributions (to SSE-like functions) with different spatial supports, surrounded by a fully coupled plate interface.The associated surface displacements (i.e., the inverted data) are strongly and randomly perturbed according to a normal probability distribution given by the data variance per component, which we took as 2.1, 2.5 and 5.1 mm in the north, east and vertical directions, respectively (Radiguet et al., 2011).
Figure 4 shows the target slip models and both, the associated exact displacements (blue arrows) and the perturbed ones (red arrows).The data uncertainty is represented by the gray ellipses at the tips of the perturbed vectors, the semiaxes corresponding to the standard deviation of the normal distribution used to perturb the data per component.The interplate coupling corresponds to three-months cumulative backslip assuming a 6 cm/yr plate convergence (i.e., 1.5 cm), and the geometry and position of the three Gaussian slip patches were inspired by recent SSE solutions found in the region (Cruz-Atienza et al., 2020).Please notice how large are the perturbations.
Inversions for the three Gaussian slip models were done for both the exact and perturbed data.Each set of data was inverted without regularization, F = I, and with correlation lengths of 10, 20, 30, 40, 50 and 60 km.In all cases backlip restrictions were applied by means of the GPM so the interplate coupling could never overcome the value of one.Figure 5A shows some slip solutions for the largest-Gaussian exact data together with the associated restitution maps.Although the data fit is excellent in all cases, acceptable solutions are only retrieved when model regularization is applied.For L = 30 km, the ari is above 0.9 so that the slip solution is almost perfect, except along the Gaussian contour where there is an unrealistic slip discontinuity in the target model (the same situtation as for the checkeboard tests).
When random noise is added to the theoretical observations and the inverse problem is solved by integrating the data uncertainty by means of the precision matrix, the model regularization becomes even more critical to achieve a reliable solution.This can be seen in Figure 5B, where the restitution is very poor around the Gaussian slip area when no regularization is applied as compared with that for L = 40 km, where the ari is also above 0.9 and thus the slip solution is surprisingly good.Also astonishing, results for the other two, smaller Gaussian slip models were very similar (see Appendix C, Figures S1 and S2).A summary of the 42 inversions (14 per Gaussian model) is shown in Figure 6, where we see that although the data-fitting errors for the noisy inversions are roughly four times larger than those obtained from the exact data, the ari in all cases is above 0.9 for the best solutions (i.e., for the optimal L) even for the smallest and circular Gaussian case, which has a slip characteristic length smaller than 80 km centered at 38 km depth (Figure 4A).

The 2006 Guerrero SSE
During the 20 years preceding the devastating 2017 Mw8.2 Tehuantepec earthquake that took place offshore the Oaxaca state, Mexico, long term SSEs in Guerrero occurred almost every four years (i.e., six events between 1998 and 2017) and had remarkably large moment magnitudes (Mw>7.5)(Kostoglodov et al., 2003;Radiguet et al., 2012;Cruz-Atienza et al., 2018).After the earthquake, the regional plate-interface SSE beating has strongly changed so that two other SSEs took place in that state in the next two years (in 2018 and 2019) with much smaller magnitudes (Mw ∼ 7.0) (Cruz-Atienza et al., 2020).Among all Mexican SSEs, the 2006 Guerrero event has been the most investigated despite the poor GPS instrumentation on that time (Kostoglodov et al., 2010;Vergnolle et al., 2010;Radiguet et al., 2011Radiguet et al., , 2012;;Cavalié et al., 2013;Bekaert et al., 2015;Villafuerte and Cruz-Atienza, 2017).One of its most interesting features is that, unlike adjacent subduction segments, the slow slip seems to have penetrated the updip seismogenic region of the plate interface up to 15 km depth in the Guerrero seismic gap.In this section we perform a thorough analysis of the inverse problem resolution for that event and provide what we think are its most reliable features as compared with previous results reported in the literature.

Resolution
In previous sections we found that the problem resolution depends on two main parameters: (1) the slip characteristic length (PS) and ( 2) the inverse-problem correlation length (L).This is true for a given problem geometry (i.e., for a stations array and plate interface geometry).For this reason, we can determine fault regions where resolution (i.e., the restitution index) is high enough for a given L and PS, which means that the inverted slip in those regions is valid within the wavenumber bandwidth associated to the von Karman spectrum for that L. Since only 12 significant GPS stations registered the 2006 SSE, we performed three different M-CB tests considering only the location of these sites and CB periodic c-slip values of -8 and 25 cm.The tests were done for checkerboard unit lengths (PS) of 80, 100 and 120 km, and for L = 0 (no regularization), 10, 20, 30, 40, 50 and 60 km.These resolution exercises assumed reasonable backslip and rake angle restrictions so that the c-slip component ranges within [−8, 32] cm and the angle within [20, −20] • with respect to the c-direction, which implies the admissible range of [−3, 3] cm for the p-slip component.
Plate-interface resolution maps (i.e., for the mcri metric) are shown in Figure 7 as a function of PS and L. As expected, overall mcri values increase with PS for a given L.Although less evident, they also increase with L for a given PS up to a certain correlation value.However, supplementary results not shown reveal that, in the latter case, the high-resolution regions stop expanding for L above 30 km for all three PS cases.The maps show isocontours for mcri = 0.6, which delineate fault regions where the slip solutions are likely to resolve the actual slip within 40% error.As explained previously, these maps represent a lower resolution bound because the M-CB tests assume unrealistically sharp slip discontinuities that strongly penalize the restitution index along the boundaries of the square slip patches (e.g., see Figure 1).For this reason, we expect the resolution within the regions to be higher than the mcri isocontours value.Either way, even in the M-CB test for the maximum PS and L values, the high resolution region does not extend across the whole expected SSE area, as claimed by previous authors using different inversion techniques (Radiguet et al., 2011).Our resolutions maps represent the key piece allowing us to tell something reliable (to some point) about the 2006 SSE.
Figure 8 summaries the results from all M-CB tests in terms of the average mcri and data-misfit L2 error.Although errors are similar for all slip characteristic lengths PS, the maximum average mcri value increase with PS and follow a concave trajectory with L as previously noticed from Figure 7.However, unlike the previous M-CB exercises considering all currently available geodetic stations (Figures 2 and 3), the optimal correlation lengths (i.e., those maximizing the restitution) increase with PS.This remarkable and reasonable result is due to both the slip restrictions and the sparsity of the stations.It tells us that, depending on the characteristic size of the SSE patch we want to solve best, the regularization of the problem must be adapted.For instance, if we are interested in SSE patches with a characteristic length of 80 km, then L = 10 km is the optimal choice.Of course, such small value is detrimental to the extent of the acceptable resolution region, as seen in Figure 7.If L = 20 km, then patches with characteristic length of 100 km will be optimally solved in a larger fault region.

SSE Inversions
The next inversions we present were done using the same GPS data as Radiguet et al. (2011).This means that the displacement time series were carefully pre-processed (Vergnolle et al., 2010) and then corrected for inter-SSE deformations by subtracting the linear trends from the period 2003-2005 per station.Thus, the resulting time series are supposed to show the deviations from the long-term steady motion during the 2006 Guerrero SSE.
Since the inter-SSE displacement trends per station are significantly different in Guerrero (Radiguet et al., 2012), it is important to emphasize that the removal of these individual (usually linear) trends from the data for SSE imaging is an incorrect practice for two simple reasons: (1) the resulting displacement (or velocity) guesstimates no longer correspond to the initial reference system (e.g., fixed tectonic plate), often leading to overestimated deformation values and thus unrealistically high SSE moment magnitudes; and (2) the non-linear transformation implied by the corrections removes the common reference frame between GPS stations, which makes the resulting data set not strictly comparable and therefore its inversion meaningless.The fact that this correction is common practice does not make it acceptable.Either way, for the sake of comparison with previous solutions using this dataset, we have inverted the corrected time series from January 30 (2006) to January 15 (2007) for four different correlations lengths (L = 10, 20, 30 and 40 km) considering slip restrictions, so that the backslip could not overcome the fullcoupling regime in that period and the rake vector could vary +/-20 • from the c-direction.
Figure 9 shows the inversion results for two optimal correlation lengths (L = 20 and 30 km).Since the data is almost perfectly explained in both cases, the preferred solution will depend on both the scale at which we are interested in for interpretations and reasonable physical considerations.Taking the 1 cm slip contour as the effective SSE area, then the moment magnitude of the 2006 event is consistent for both inversions and equal to Mw7.4.For estimating Mw, we considered a typical crustal rigidity µ = 32 × 10 9 Pa.
As shown in the last section, given the poor GPS coverage during the 2006 SSE, the inverse problem regularization plays a critical role to have some confidence in the slip solutions.In the absence of resolution analysis, it is difficult to justify any conclusion, especially between distant stations.For instance, the absence of data along most of the north-west Guerrero seismic gap (NW-GGap) (i.e., between ZIHP and CAYA) (UNAM, 2015) and the Guerrero Costa Chica (i.e., between CPDP and PINO) is unfortunate and obliges us to be cautious in the interpretations.Previous investigations concluded that SSEs behave differently between these two Guerrero subduction segments so that, unlike the Costa Chica, the slow slip in the NW-GGap reaches the updip seismogenic interface zone (i.e., up to 15 km depth) (Radiguet et al., 2011;Cavalié et al., 2013) releasing aseismically a significant part of the accumulated inter-SSE strain (Radiguet et al., 2012;Bekaert et al., 2015).
Figure 10 shows a comparison between our preferred solution (model A) (i.e., for L = 30 km) and two previously published solutions, one from the simultaneous inversion of both GPS and InSAR data (Model B by Cavalié et al. (2013)) and the other from GPS data only (model C by Radiguet et al. (2011)).Our solution is shown together with the associated 60% resolution regions (regions where the average mcri is higher than 0.6), which are taken from Figure 7 according to the optimal solutions of Figure 8. Confidence contours thus delineate the fault regions where solutions should disagree with the actual slip by less than 40% in different wavenumber bandwidths depending on L. The red contours delineate the 60% confidence regions for a slip characteristic length of 80 km, while the green contours depict the same regions for a 120 km characteristic length.Although the three slip solutions are in general consistent, there are clear differences among them.The most visible are (1) the concentration of separated patches in model C (i.e., one of them far from the coast and below 40 km depth, and another one to the east) which may be artifacts due to a lack of regularization (Bekaert et al., 2015), which are consistent north of the CAYA and COYU stations, and (2) the peak slip values that range between 20 and 25 cm.Moment magnitudes are also slightly different (i.e., 7.4 and 7.6 for models A and C, respectively).However, all three models coincide on the updip SSE spreading west of station CAYA, where our model has resolution higher than 60% up to a distance no more than 30 km west of that station.This region is of critical importance because it extends along the NW-GGap, where recent onshore and offshore observations show that slow earthquake indeed happen in a particular way, and thus where the mechanical properties of the plate interface could be different (Cruz-Atienza et al., 2020;Plata-Martínez et al., 2020).Models B and C are remarkably different between stations ZIHP and CAYA, where the InSAR data used for model B does not play any significant role.West of this region, model B predicts a very large shallow penetration of the SSE across the mechanically stable zone where M7+ earthquakes occur every ∼35 years (see past rupture areas) (UNAM, 2015).For this reason, model C, which is consistent with our model A, is the most plausible solution for that zone.Besides, our resolution close to the ZIHP station is higher than 60% as well.In conclusion, our preferred ELADIN solution has the most reliable features of both previously published slip models.

Conclusions
We have introduced the ELADIN method, a new fault-slip inversion technique based on the adjoint elastostatic equations under a constrained optimization framework.The method main characteristics are the projection operator built with the von Karman autocorrelation function to control the spectral content of the solution, i.e, the problem regularization, and the gradient projection method to impose physically-consistent slip restrictions (e.g., interplate coupling smaller than any given value and rake angles consistent with the relative plate motion).To account for the data uncertainty, the method weights the observations according to their individual variances through the precision matrix.Synthetic slip inversions from strongly perturbed data show that the model restitution across the plate interface is surprisingly high (for both SSE and coupled interface regions) when this uncertainty is taken into account.The ELADIN method thus allows determining the aseismic slip on any 3D plate interface (or any fault surface) by simultaneously inverting relaxing slip and coupled fault areas with a spectral control of the problem solution that guaranties a given resolution criterion.We defined this criterion by means of the mobile checkerboard restitution index (mcri), which allows determining fault regions where the resolution (i.e., the slip restitution index) is high enough for a given von Karman autocorrelation length, L. This means that the inverted slip in those regions is valid (to some desired extent) within the wavenumber bandwidth associated to the von Karman spectrum for that L.
After a thorough resolution analysis of the study region, we inverted the 2006 Guerrero SSE.Our preferred slip model (Model A), obtained for L = 30km, was compared with two previously published solutions and found that it retains the most reliable features of these two models.On one hand, our model is consistent with the solution of Cavalié et al. (2013) (Model B) in that it places the maximum slip region above 40 km depth (i.e., downdip from stations CAYA and COYU), where this solution is well constrained by the InSAR data.On the other, although all solutions predict the SSE shallow penetration along a large part of the NW-GGap segment (west of CAYA), our resolution analysis clearly shows that this penetration might not be a reliable feature of the 2006 SSE.However, our Model A is much closer to the solution of Radiguet et al. (2011) (Model C) close to station ZIHP, where only GPS data is available.In this sense and considering also that M7+ earthquakes occur every ∼35 years east from that station (see previous rupture areas in Figure 10), which implies that the plate interface is mechanically unstable, then the extremely large updip SSE penetration predicted by Model B (Cavalié et al., 2013) between stations ZIHP and CAYA seems unrealistic.
A systematic application of the ELADIN method has been made in two associated works (Cruz-Atienza et al., 2020;Villafuerte et al., 2021) to invert recent GPS and InSAR data from the large array shown in Figure 1, in the period 2016-2020, where four major earthquakes and multiple SSEs occurred throughout the Mexican subduction zone.that t i < t i+1 for i ∈ {1, 2, . . ., l −1}.With this set, we construct a set of intervals like {[0, t 1 ], [t 1 , t 2 ], . . ., [t l−1 , t l ]}.Suppose that we have not found the minimizer up to the interval [t j−1 , t j ], then we can model the slip along that interval as where If we substitute eq. ( 24) in the quadratic cost function ( 19), we leave it as a function of ∆t which can be reformulated as where The solution of this problem is Only one of the following three cases can occur (i) If g j−1 > 0 the minimizer is at ∆t * = 0 with t * = t j−1 and p * = p j−1 .
(iii) If ∆t * > t j − t j−1 then try the next interval.
Once the optimal step has been found, ∆t * , the Cauchy point is evaluated as Appendix C. Gaussian slip inversions       4 in terms of the whole-interface average restitution index, ari, and average data-misfit error (red) as a function of the inversions correlation length L. Solid lines correspond to the inversions using the exact data while dashed lines to the inversions with noisy data (see Figure 4).Notice that in all cases the maximum restitutions, ari, are above 0.9.

4.
With the gradient use any iterative optimization algorithm to find an update step ∆D k 5. Update the slip D k+1 = D k + ∆D k .end 2 nd Step: Constrained slip inverse problem (Gradient Projection Method) Data: Optimal solution of 1 st step, D * Project D * into the physically-consistent domain to get the initial solution D 0 ; while Convergence is not achieved do 1.From D k compute the Cauchy point D c k (details in Appendix B) 2. Update the slip D k+1 = D c k .end

Figures
Figures S1 and S2show the synthetic data inversions and restitution indexes with and without noise of the Gaussian-like pulses shown in Figures4A and 4B, respectively.

Figure 1 :
Figure 1: Checkerboard inversions for PS of (A) 60, (B) 80 and (C) 100 km, and correlation length, L, of 20 km.The inverted slip along the plate convergence direction, c-slip, with the surface displacement fits (left column) and the associated restitution index (right column) are displayed on the 3D plate interface (gray contours).Green triangles are the GPS stations.

Figure 2 :Figure 3 :
Figure2: M-CB tests for PS of (A) 60, (B) 80 and (C) 100 km and correlation length, L, of 20 km.Distributions of mcri (first row), the optimal correlation length (second row) and the multiscale assembly of the restitution index (computed from the assembly of the best slip solutions for the CBs shown in Figure1), all of them computed with the c-slips inverted and displayed on the 3D plate interface (gray contours).Green triangles are the GPS stations.

Figure 4 :
Figure 4: Slip models along the c-direction on the plate interface (background colors) and the associated model displacement predictions (arrows) for three Gaussian-like slip patches with different characteristic lengths.Blue and black-solid arrows show the exact surface displacements while red and black-dashed arrows show the same predictions but stochastically perturbed according to the normal distributions given by the data variance per component.

Figure 5 :
Figure 5: Synthetic inversion results for the c-slip model shown in Figure 4C from the exact target displacements (panel A) and from the perturbed (noisy) displacements (panel B).The second row of each panel shows the distribution of the restitution index for the c-slip over the plate interface without regularization and for different values of the correlation length, L.

Figure 6 :
Figure6: Synthetic inversion results for the three Gaussian-like slip functions shown in Figure4in terms of the whole-interface average restitution index, ari, and average data-misfit error (red) as a function of the inversions correlation length L. Solid lines correspond to the inversions using the exact data while dashed lines to the inversions with noisy data (see Figure4).Notice that in all cases the maximum restitutions, ari, are above 0.9.

Figure 7 :Figure 8 :
Figure 7: Plate interface distribution of the mobile checkerboard restitution index, mcri, of the c-slip inverted from M-CB tests corresponding to patch sizes (PS) of 80, 100 and 120 km and correlation lengths L = 10, 20 and 30 km for the 2006 SSE stations configuration.Black contours correspond to mcri values of 0.6 (i.e., slip resolution of 60%).

Figure 9 :
Figure 9: Aseismic slip inversions, in the plate convergence direction, of the 2006 Guerrero SSE for correlation lengths L = 20 km (A) and L = 30 km (B).The plate interface coupling is determined from the ratio between the back slip and the cumulative slip in the inverted period given a plate convergence rate of 6 cm/yr.

Figure 10 :
Figure 10: Comparison of our preferred solution (model A -for L = 30 km, Figure 9) with two previously published models for the 2006 Guerrero SSE, model B from Cavalié et al. (2013) and model C from Radiguet et al. (2011).60% resolution contours for slip-patch (PS) characteristic lengths of 80 and 120 km are shown over model A.

Figure S1 :
Figure S1: Synthetic inversion results along the c-direction for the Gaussian-like slip model shown in Figure 4A from the exact target displacements (panel A) and from the perturbed displacements (panel B).The second row of each panel shows the distribution of the restitution index over the plate interface without regularization and for different values of the correlation length, L.

Figure S2 :
Figure S2: Synthetic inversion results along the c-direction for the Gaussian-like slip model shown in Figure 4B from the exact target displacements (panel A) and from the perturbed displacements (panel B).The second row of each panel shows the distribution of the restitution index over the plate interface without regularization and for different values of the correlation length, L.