# A spectral expansion approach for geodetic slip inversion: implications for the downdip rupture limits of oceanic and continental megathrust earthquakes

A spectral expansion approach for geodetic slip inversion: implications for the downdip rupture... SUMMARY We have developed a data-driven spectral expansion inversion method to place bounds on the downdip rupture depth of large megathrust earthquakes having good InSAR and GPS coverage. This inverse theory approach is used to establish the set of models that are consistent with the observations. In addition, the inverse theory method demonstrates that the spatial resolution of the slip models depends on two factors, the spatial coverage and accuracy of the surface deformation measurements, and the slip depth. Application of this method to the 2010 Mw 8.8 Maule Earthquake shows a slip maximum at 19 km depth tapering to zero at ∼40 km depth. In contrast, the continent–continent megathrust earthquakes of the Himalayas, for example 2015 Mw 7.8 Gorkha Earthquake, shows a slip maximum at 9 km depth tapering to zero at ∼18 km depth. The main question is why is the maximum slip depth of the continental megathrust earthquake only 50 per cent of that observed in oceanic megathrust earthquakes. To understand this difference, we have developed a simple 1-D heat conduction model that includes the effects of uplift and surface erosion. The relatively low erosion rates above the ocean megathrust results in a geotherm where the 450–600 °C transition is centred at ∼40 km depth. In contrast, the relatively high average erosion rates in the Himalayas of ∼1 mm yr–1 results in a geotherm where the 450–600 °C transition is centred at ∼20 km. Based on these new observations and models, we suggest that the effect of erosion rate on temperature explains the difference in the maximum depth of the seismogenic zone between Chile and the Himalayas. Satellite geodesy, Inverse theory, Intra-plate processes 1 INTRODUCTION Over the past few decades, advances in space geodesy such as Interferometric Synthetic Aperture Radar (InSAR) and global positioning system (GPS) have facilitated detailed mapping of surface deformation. These high quality and regional datasets have enabled scientists to build 3-D deformation fields for many large earthquakes. Usually, the fault geometry can be reasonably approximated from aftershocks, geophysical imaging and geological field studies. Elastic models that relate slip at depth to surface deformation can then be constructed by discretizing this fault into small patches before estimating the slip on each patch. While this method is straightforward, it hardly provides an estimate of the spatial resolution nor unbiased uncertainties associated with the slip model. It is well known that the equations which relate slip at depth to surface deformation contain an upward continuation term (Steketee 1958), which exponentially attenuates the rupture signal at wavelengths smaller than the depth, therefore making this inverse problem inherently non-unique (Parker 1994). This is true even if the data coverage is complete and the data have no errors. The non-uniqueness becomes evident in parametrized inversions when the discretized patch size is smaller than the depth; the inversions are poorly conditioned resulting in wild oscillations in slip distribution. In order to stabilize the solution, smoothness regularization combined with a non-negativity constraint is applied to the inversion, which acts as a strong prior knowledge to the model (Jonsson et al. 2002; Simons et al. 2002; Fialko 2004; Tarantola 2005). In some studies, a checkerboard test is used to illustrate the resolving power of the data (Fialko 2004; Tong et al. 2010). Tong et al. (2010) performed this test for the Mw 8.8 2010 Maule earthquake and found the model can only recover a 20-km checker up to about 35 km depth, but can recover a 40-km checker to 60 km depth. These results are consistent with the inverse relation between model resolvability and rupture depth. For several recent megathrust earthquakes, surface deformation has been well recorded and many studies have incorporated these data into earthquake rupture inversions using the standard parameter estimation approach. Inversions for the 2010 Mw 8.8 Maule earthquake (oceanic plate subducting) suggest a tapering of slip towards zero at ∼40 km depth (Tong et al. 2010; Lorito et al. 2011). This depth is similar to many other oceanic subduction slip inversions [e.g. the 2011 Mw 9.0 Tohoku-oki (Simons et al. 2011; Minson et al. 2014); 2004 Mw 9.2 Sumatra–Andaman (Ammon et al. 2005; Lay et al. 2005; Chlieh et al. 2007)]. By contrast, inversions for the 2015 Mw 7.8 Nepal earthquake (continental plate subducting) suggest a tapering of slip towards zero at less than 20 km (Galetzka et al. 2015; Wang & Fialko 2015). This depth is consistent with the maximum depth of strong interseismic coupling resolved by geodetic models along the Himalayan front (Stevens & Avouac 2015). This difference raises two questions. First, is the difference in downdip rupture limit well resolved by the observations? If it is well resolved, what physical mechanisms may explain such a large difference? This requires a real uncertainty estimate as well as an analysis of resolving power for inversion models. To address the first question, we adopt the approach of linear inverse theory (Gilbert 1971; Parker 1977). As noted by Parker (1977): ‘The quality that distinguishes inverse theory from the parameter estimation problem of statistics is the unknowns are functions, not merely a handful of real numbers. This means that the solution contains in principle an infinite number of variables, and therefore with real data the problem is as underdetermined as it can be.’ The approach begins with the well-known Green's functions relating slip at depth on a prescribed fault plane to surface deformation. Following the spectral expansion method, we then construct a set of up to N orthonormal slip functions that span the set of N surface observations. The solution and its uncertainty are constructed from a finite number M ≤ N of these functions, or kernels, that match the data to within their uncertainty. The low order kernels are generally smoothest and well constrained by the data while the higher order kernels have more oscillations and are more poorly constrained. The solution is constructed from the well-resolved kernels. We use this approach to provide bounds on the depth distribution of slip during the Mw 8.8 Maule (oceanic) and Mw 7.8 Gorka (continental) megathrust areas, and show that the factor-of-two difference in slip depth is indeed well resolved by the geodetic data. To address the second question, the pronounced difference in the slip depth, we develop a simple 1-D thermal model that includes the effects of vertical advection of heat due to erosion-induced rebound (Royden 1993). We find that an erosion rate of 1 mm yr–1 decreases the depth of the 600 °C isotherm by a factor of two, which can explain the relatively shallow rupture depth of the Gorka earthquake. 2 THE SPECTRAL EXPANSION APPROACH Here we give a brief summary of the spectral expansion theory from Parker (1977, 1994) as applied to static slip inversion. The theory is originally from Gilbert (1971). The main difference between our formulation and the original formulation is that our model is 2-D vector slip (strike and dip) as a function of the two dimensions on the fault surface m(x). Note that the fault does not need to be a planar surface. The data are surface deformation sampled at a finite number of points on the surface of the Earth. GPS can provide 3 components of deformation while InSAR provides 1 or 2 line of sight components. Each component of each data point is di and has an uncertainty σi. We use the elastic half-space Green's function (Okada 1985) to relate vector slip at depth to each component of surface deformation. We note that Green's functions for a more complicated elastic structure could be used. The convolution of the model and the Green's function is given by   $${d_i} = \mathop \int \nolimits_{\rm{\Sigma }}^{} {G_i}\left( {\boldsymbol{x}} \right)\ {\boldsymbol{m}}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }},\$$ (1)where Σ represents integration over the entire fault plane. After dividing both sides of eq. (1) by the uncertainty σi the equation can be rewritten as   $$d{'_i} = \mathop \int \nolimits_{\rm{\Sigma }}^{} {G'_i}\left( {\boldsymbol{x}} \right){\boldsymbol{m}}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }}.\$$ (2) Now the data d'i as well as the Green's functions G'i are dimensionless and the data have unit standard deviation. Following Parker (1977) we form the Gram matrix $${{\boldsymbol \Gamma }}$$  $${{\rm{\Gamma }}_{ij}} = \mathop \int \nolimits_{\rm{\Sigma }}^{} {G'_i}\left( {\boldsymbol{x}} \right){G'_j}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }}.\$$ (3) Thus, $${{\boldsymbol \Gamma }}$$ is positive-definite and symmetric, and can be diagonalized with an orthonormal matrix O of eigenvectors   \begin{eqnarray} {{{\bf O}}^{\rm{T}}} {\boldsymbol \Gamma} {\bf O}\ &=& {\boldsymbol{\ }}{{\boldsymbol \Lambda }},\quad {\rm{where\ }}{{\boldsymbol \Lambda }}\ = {\rm{\ diag}} \{ {{{\rm{\lambda }}_1},{{\rm{\lambda }}_2}, \ldots ,{{\rm{\lambda }}_{N\ }}} \}\nonumber\\ &&\qquad{\rm{\ and\ }}{{\rm{\lambda }}_1} \ge {{\rm{\lambda }}_2} \ge {{\rm{\lambda }}_3} \ge \ldots {{\rm{\lambda }}_{N\ }} > 0, \end{eqnarray} (4)where λi are the eigenvalues; the full set of eigenvalues is called the spectrum. It's not difficult to show that   \begin{eqnarray} {\lambda _i} \!=\! {O_i}^{\rm{T}}\ {{\boldsymbol \Gamma }}\ {O_i} &=& \mathop \sum \limits_j {O_{ji}}\mathop \int \nolimits_{\rm{\Sigma }}^{} {{\bf G}}{{{\bf G}}^{\rm{T}}}{\rm{d\Sigma }}{O_{ij}}\nonumber\\ &=& \mathop \int \nolimits_{\rm{\Sigma }}^{} \mathop \sum \limits_j {O_{ji}}{{\bf G}}{{{\bf G}}^{\rm{T}}}{O_{ij}}{\rm{d\Sigma }},\ \end{eqnarray} (5)where G is a column vector of normalized Green's functions and the Gram matrix $${{\boldsymbol \Gamma }}$$ can be treated as a dyad matrix. We rewrite the equation as following:   $${{\boldsymbol \Lambda }}\ = \mathop \int \nolimits_{\rm{\Sigma }}^{} {{{\bf O}}^{\rm{T}}}{{\bf G}}{{{\bf G}}^{\rm{T}}}{{\bf O}}\, {\rm{d\Sigma }}.\$$ (6) Now consider the kernel functions defined by   $${\psi _i}\left( {\boldsymbol{x}} \right) = {\lambda _i}^{ - \frac{1}{2}}\ \mathop \sum \limits_j {O_{ji}}{G'_j}\left( {\boldsymbol{x}} \right).$$ (7) It can be easily verified from the eq. (5) that ψi(x) are an orthonormal set   $$\mathop \int \nolimits_{\rm{\Sigma }}^{} {\psi _i}\left( {\boldsymbol{x}} \right){\psi _j}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }} = {\delta _{ij}}.\$$ (8) Therefore, we consider the expansion of the model m(x) in terms of these kernels   $${\boldsymbol {m}}( x) = \mathop \sum \limits_j {a_i}\ {\psi _i}\left( {\boldsymbol{x}} \right) + {\psi _*},$$ (9)where ψ* is the annihilator (i.e. models that produce no surface deformation at the positions of the data points). The coefficients of the expansion can be derived   \begin{equation*} {a_i} = \mathop \int \nolimits_{\rm{\Sigma }}^{} {\psi _i}\left( {\boldsymbol{x}} \right){\boldsymbol{m}}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }} \end{equation*}   \begin{equation*} \quad= \mathop \int \limits_{\rm{\Sigma }}^{} {\lambda _i}^{ - \frac{1}{2}}\mathop \sum \limits_j {O_{ji}}{G'_j}\left( {\boldsymbol{x}} \right){\boldsymbol{m}}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }}\ \end{equation*}   \begin{equation*} \quad= {\lambda _i}^{ - \frac{1}{2}}\ \mathop \sum \limits_j {O_{ji}}\mathop \int \nolimits_{\rm{\Sigma }}^{} {G'_i}\left( {\boldsymbol{x}} \right){\boldsymbol{m}}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }} \end{equation*}   $$\quad= {\lambda _i}^{ - \frac{1}{2}}\ \mathop \sum \limits_j {O_{ji}}d{'_j}.$$ (10) Since all the d'j have unit variance and the matrix O is orthonormal, the error estimates of ai are simply λi − 1/2 and ai are statistically independent. When the inverse problem is handled with this approach, it is easy to isolate the parts that are better determined, that is those have smaller error or larger eigenvalue. 3 APPLICATION TO 2-D AND 3-D TEST CASES We first test this approach using an infinitely-long vertical strike-slip fault. Two input slip distributions are used for testing. The first is uniform slip (1 m) with depth between the surface and 12 km (Fig. 1c, black dashed). The second is a Gaussian-shaped slip centred at 8 km depth (Fig. 1, black dashed). We then use the 2-D half-space Green's function (Cohen 1999) to forward generate surface observations from –50 to 50 km distance with 200 m sampling increment and add a random noise of 2 per cent of the maximum shear deformation (Fialko 2004). Next, we normalize the Green's function, compute the Gram matrix analytically and decompose it to get the eigenvalues (Fig. 1a) as well as the orthonormal kernels (Fig. 1b). Note that the orthonormal kernels have oscillations as a function of depth with more oscillations for higher kernel number. After this step, we calculate the coefficient for each kernel and combine them to get the inverted slip (Fig. 1c, blue line). The uncertainty of the inverted slip can be acquired by propagating the error of each coefficient to the model (Fig. 1c, magenta dash--dotted line). Figure 1. View large Download slide 2-D infinite fault test case. (a) Eigenvalues versus order, with the blue circle denoting the order selected for the step case (c) and blue star for the Gaussian case (d) when using noisy data. (b) Some examples of orthonormal kernels versus depth. (c) Recovered slip versus depth plot for input slip being a step function and (d) recovered slip versus depth plot for input slip being a Gaussian function. For (c) and (d), the black dashed lines represent the input slip, blue lines are slip inverted using noisy data with magenta dash--dotted lines being the uncertainty bounds, light blue lines are the slip inverted using noise-free data, transparent red and green lines are inverted slip using noise-free data with more kernels. The forward model and fitting can be found in Fig. S1. Figure 1. View large Download slide 2-D infinite fault test case. (a) Eigenvalues versus order, with the blue circle denoting the order selected for the step case (c) and blue star for the Gaussian case (d) when using noisy data. (b) Some examples of orthonormal kernels versus depth. (c) Recovered slip versus depth plot for input slip being a step function and (d) recovered slip versus depth plot for input slip being a Gaussian function. For (c) and (d), the black dashed lines represent the input slip, blue lines are slip inverted using noisy data with magenta dash--dotted lines being the uncertainty bounds, light blue lines are the slip inverted using noise-free data, transparent red and green lines are inverted slip using noise-free data with more kernels. The forward model and fitting can be found in Fig. S1. An additional test is done to explore the recoverability by inverting the synthetic data with no noise. The light blue line in Fig. 1(c) shows the solution using the exact same number of kernels, while the transparent red and green lines are solutions using more. It's obvious that with more kernels, the solution will get closer to the step function we input, despite the fact that due to the upward continuation, a full recovery of the step function is not expected. Also, the rapid reduction in the eigenvalue causes the error from the data to be exponentially magnified when it comes to higher order. To further examine this approach, we set a Gaussian input slip centring 8 km depth with the maximum slip being 1 m (Fig. 1d). The conclusion remains the same. In both cases, the surface deformation is well fit with the selected kernels determined from the misfit function $${\chi ^2} = \frac{1}{N}\ \mathop \sum \nolimits_{i\ = \ 1}^N {(\frac{{{d_i} - d_i^{\rm{pred}}}}{{{\sigma _i}}})^2}\$$ (Fig. S1). Note that from the theory (Section 2), if the same fault geometry is determined, once the locations of observations are set, every computation step will be the same until the calculation of coefficients for the kernels. In other words, both cases shown in Figs 1(c) and (d) are just different recompositions of the kernels shown in Fig. 1(b), based on different surface observations. After considering the 2-D test, we expanded the formulation to a 3-D test case. We use a 350-km-long 250-km-wide planar fault with N30°E strike and 15° dip towards east. The input dip-slip is a Gaussian function shown in Fig. 2(a), with the corresponding moment magnitude being 8.0. Then we generate a forward model using the half-space Green's function (Okada 1985) and project the vector surface deformation into descending and ascending line-of-sight (LOS) directions to simulate the InSAR observations. Random noise (2 per cent of maximum amplitude) (Fialko 2004) is added and then the data are down sampled with a quad-tree algorithm (Jonsson et al.2002) based on the curvature of the deformation field. Due to the complexity of the 3-D half-space Green's function (Okada 1985), the computation of the Gram matrix is done with numerical integration for both strike and dip slip. We then decompose the Gram matrix to recover the eigenvalues and orthonormal kernels (Figs S3i and S6). We found the 3-D case has a much slower fall-off in misfit with increasing number of kernel functions than the 2-D case. To determine how many kernels to use, we examine the misfit function (Fig. 2c) and chose the smallest number with reasonable misfit reduction (e.g. 100 kernels), in which case, more kernels will result in an explosion in model's uncertainty but no significant improvement on data fitting. While the along-strike cumulative dip slip shows a good fit to the input (Fig. 2b) with the recovered moment magnitude being 7.99, a small component of strike slip was also ‘recovered’, potentially due to incomplete coverage from down-sampled data and the noise we added. Besides, we also discovered some short wavelength undulations in the inverted slip, which we believe is due to the nature of this decomposition method. Figure 2. View large Download slide 3-D test case. (a) Input dip slip with positive defined as hanging wall moving up-dip. The black dashed line is the fault trace. (b) Cumulative dip slip $$p( z )\ = \mathop \int \nolimits_0^L m( {x,z} ){\rm{d}}x\ \$$ (Simons et al. 2002) versus depth plot with black dashed line being the input dip slip, blue line representing the inverted dip slip and magenta dash-dotted lines denoting the uncertainty. (c) Misfit versus number of kernels, with the red dot being the order selected. (d) Inverted dip slip and (f) inverted strike slip with hanging wall moving north-eastward defined as positive. Resolving wavelength for dip slip (e) and strike slip (g) constructed by the algorithm shown in Fig. S2. Data fitting and eigenvalues can be found in Fig. S3. Figure 2. View large Download slide 3-D test case. (a) Input dip slip with positive defined as hanging wall moving up-dip. The black dashed line is the fault trace. (b) Cumulative dip slip $$p( z )\ = \mathop \int \nolimits_0^L m( {x,z} ){\rm{d}}x\ \$$ (Simons et al. 2002) versus depth plot with black dashed line being the input dip slip, blue line representing the inverted dip slip and magenta dash-dotted lines denoting the uncertainty. (c) Misfit versus number of kernels, with the red dot being the order selected. (d) Inverted dip slip and (f) inverted strike slip with hanging wall moving north-eastward defined as positive. Resolving wavelength for dip slip (e) and strike slip (g) constructed by the algorithm shown in Fig. S2. Data fitting and eigenvalues can be found in Fig. S3. To further understand the spatial resolution of the inversion, we developed a technique to qualitatively represent the resolving power of the slip model (Fig. S2). For each kernel used in the construction of the model, we pick out the peaks and troughs, together with the corners and sides, to create a set of nodes. Next, we apply Delaunay triangulation (Delaunay 1934) to these nodes and set their values to the average of the distances between each node and its adjacent nodes. Then we fit these distance measures at each vertex using splines in tension (Smith & Wessel 1990) and filter the result to get the resolving wavelength map for the kernel. When applying this to all the kernels used in the spectral expansion inversion, the minimum value among all the kernels at every position on the fault plane is considered as the resolving wavelength for the model (Figs 2e and g). We also performed further tests for dip slip and a mixture of dip and strike slip scenarios, the results of which are presented in supplementary information (Figs S7 and S8). The two test cases reveal some interesting features of this spectral expansion approach. First, the eigenvalues derived from the Gram matrix decrease slower in the 3-D case than the 2-D case. We believe this is because the data in the 3-D case requires extra constraints for one more dimension, thus more eigenvalues are needed to achieve the same level of trade-offs. Second, the resolving power of the model depends a lot on data sampling. It's obvious, even from the equations in Section 2, that the shapes of kernels are completely determined by the model geometry and distribution of the data samples. The derived resolving wavelength could also act as a guide for the digitization of parameter estimation models. Third, the 2-D case shows an increase in uncertainty at shallower depth since it has complete surface coverage. The 3-D case shows a similar feature, but suddenly the uncertainty becomes smaller at the surface, which we relate to the lack of near-fault data. The truth is that the data can never be dense enough when it comes close to the fault, as the Green's function will change dramatically over a short distance. This inversion approach relies heavily on the similarity of Green's function between data points, that is the redundancy of data. Thus, in this case, the quad-tree algorithm may not be optimal for down sampling, as the redundancy of data is dependent on fault depth. Fourth, although the computation is defined as the inner product of the Green’ functions over the model space, that is how similar their Green's functions are, the Gram matrix could potentially represent the covariance (after scaling by the data's uncertainty) between data points, if we assume the data having zero mean. In other words, if the measurement of one data point changes away from zero due to an earthquake, the measurement of a nearby data point may also change correspondingly. Finally, the fault could be in any shape as long as the integration shown in eq. (3) can be done properly, which gives a strong adaptability to this approach. 4 APPLICATION TO THE 2010 MAULE AND THE 2015 GORKHA EARTHQUAKES We first use this spectral expansion approach to model the slip for the Mw 8.8 Maule, Chile earthquake. This earthquake was a megathrust event that ruptured a mature seismic gap (Moreno et al.2010) and generated a strong tsunami throughout the Pacific Ocean (Fritz et al.2011). Previous study shows coseismic rupture extending downdip to a depth of ∼40 km, which correlates well with interseismic coupling along the Andean subduction zone in south-central Chile (Lorito et al.2011; Métois et al. 2012). Here we apply the spectral expansion approach to the processed InSAR and GPS data presented in (Tong et al.2010) and analyse the uncertainty and resolving power of the model. We use the same fault geometry as in their study, with the fault plane extending 680 km along strike N16.5°E and 250 km downdip 15° towards east. From the misfit function analysis (Fig. 3e), we determined that 150 kernels were needed to represent the coseismic slip. The solution (Figs 3a and b) shows a total slip moment corresponding to moment magnitude 8.79 assuming an average shear modulus of 40 GPa, with the thrust component being 1.86  × 1022Nm and the strike component being 5.85  × 1021 Nm. This is close to the seismic moment reported by (NEIC 2010) (Mw 8.8). From the resolving wavelength analysis (Figs 3c and d), the model has the best resolution at intermediate depth (20–40 km), with resolution decreasing dramatically outside of this range. More interestingly, at shallower depth, due to the lack of offshore data, the resolution is actually worse than at deeper depth. This brings into question the very small uncertainty estimates near the fault surface. Are we so sure the uncertainty of the slip estimate is almost zero close to fault trace? The answer is yes, with a sacrifice in resolution, the statistical precision is almost intrinsically guaranteed. The cumulative dip slip (Fig. 5a) increases from 0 to 20 km depth and decreases to zero at roughly 45 km depth. This is in good agreement with previous studies and the uncertainty of our model together with the resolving power analysis show that the reduction in slip amplitudes is well determined by observations. Figure 3. View large Download slide Inverted slip for the 2010 Mw 8.8 Maule Earthquake. (a) Inverted dip slip and (b) inverted strike slip, with same positivity definition as Figs 2(d) and (f). The black dashed line in (a) represents the fault trace. Resolving wavelength for dip slip (c) and strike slip (d). (e) Misfit versus the number of kernels used. The data fitting, uncertainty estimate and eigenvalues can be found in Fig. S4. Figure 3. View large Download slide Inverted slip for the 2010 Mw 8.8 Maule Earthquake. (a) Inverted dip slip and (b) inverted strike slip, with same positivity definition as Figs 2(d) and (f). The black dashed line in (a) represents the fault trace. Resolving wavelength for dip slip (c) and strike slip (d). (e) Misfit versus the number of kernels used. The data fitting, uncertainty estimate and eigenvalues can be found in Fig. S4. Figure 5. View largeDownload slide Cumulative dip slip versus depth plots for (a) the 2010 Mw 8.8 Maule Earthquake and (b) is the 2015 Mw 7.8 Gorkha Earthquake, with blue lines denoting the cumulative dip slip and magenta dash-dotted lines representing the uncertainty bounds. Figure 5. View largeDownload slide Cumulative dip slip versus depth plots for (a) the 2010 Mw 8.8 Maule Earthquake and (b) is the 2015 Mw 7.8 Gorkha Earthquake, with blue lines denoting the cumulative dip slip and magenta dash-dotted lines representing the uncertainty bounds. We then reinvestigate the 2015 Mw 7.8 Gorkha earthquake in Nepal. This earthquake took place on a tectonic boundary due to the Indian-Eurasian collision and ruptured 140 km on the Main Himalayan Fault (Galetzka et al. 2015). In contrast to megathrust events at subduction zones, previous studies all show a very limited downdip rupture extent with co-seismic slip typically terminating between 15 and 25 km depth (Wang & Fialko 2015). Here we apply the spectral expansion inversion to the InSAR data from (Galetzka et al. 2015; Lindsey et al.2015) and GPS data from (Galetzka et al. 2015). We down-sample the InSAR data using the quad-tree algorithm with increased sampling density closer to the fault in order to increase the resolution at shallow depth. The sampling stops 20 km away from fault in order to maintain numerical stability and save computation time. We then combine the InSAR and GPS data to perform the spectral expansion inversion. The model is created following the main frontal thrust geometry from (Ader et al.2012), with the fault plane extending 200 km along strike N71.5°W and 150 km along dip 7° towards north (Wang & Fialko 2015). By examining the misfit function, we determined that 250 kernels were needed to represent the coseismic slip (Fig. 4e). The estimated total moment for the coseismic slip (Figs 4a and b) is 6.36  × 1021 Nm, with the corresponding moment magnitude being 7.80, assuming an average shear modulus of 30 GPa, which agrees with the seismic moment reported by (NEIC 2015) (Mw 7.8). From the resolving wavelength analysis (Figs 4c and d), the resolving power generally decays with depth except in areas with denser samples due to the quad-tree sampling. The cumulative dip slip with depth plot (Fig. 5b) shows slip rapidly increasing at about 7 km depth, reaches its maximum around 10 km, and then decreases towards zero gradually. The slip versus depth profile is in agreement with the previous studies (Wang & Fialko 2015). From the uncertainty estimates and the resolving power analysis, this decay is also well resolved by the data. Figure 4. View large Download slide Inverted slip for the 2015 Mw 7.8 Gorkha Earthquake. (a) Inverted dip slip with same positivity definition as Fig. 2(d) and (b) is the inverted strike slip with positive defined as the hanging wall moving north-westward. The black dashed line in (a) represents the fault trace. Resolving wavelength for dip slip (c) and strike slip (d). (e) Misfit function versus the number of kernels used. The data fitting, uncertainty estimate and eigenvalues can be found in Fig. S5. Figure 4. View large Download slide Inverted slip for the 2015 Mw 7.8 Gorkha Earthquake. (a) Inverted dip slip with same positivity definition as Fig. 2(d) and (b) is the inverted strike slip with positive defined as the hanging wall moving north-westward. The black dashed line in (a) represents the fault trace. Resolving wavelength for dip slip (c) and strike slip (d). (e) Misfit function versus the number of kernels used. The data fitting, uncertainty estimate and eigenvalues can be found in Fig. S5. For these two cases, the InSAR data's uncertainty is taken from previous studies, with 10 cm for all InSAR data in the Maule case (Tong et al.2010) and 2.3, 5.4 and 4.1 cm for track 47, 48 and 157 in the Nepal case (Wang & Fialko 2015). One interesting feature from the spectral expansion is that the eigenvalues computed for the Maule case (Fig. S4i) decreases faster than for the Nepal case (Fig. S5i), potentially because the Nepal earthquake has a shallower dipping angle or more data close to trench. 5 DISCUSSIONS 5.1 The influence of erosion rate on thermal structure The inversion approach we developed in this study confirms that the contrasting downdip rupture extents of the two earthquakes are indeed well determined by observations. This raises the next question: Why is the downdip rupture extent of Himalayan megathrust earthquakes (15–20 km) a factor of 2 shallower than the downdip limit (>40 km) of megathrust earthquakes in oceanic subduction zones? Temperature exerts a strong influence on both dehydration reactions and deformation mechanisms. In particular, the transition from friction to intracrystalline plasticity as the dominant deformation mechanism in quartzofeldspathic rocks occurs at temperatures of 325–350 °C, and is often proposed to define the upper temperature bound for seismic behaviour in crustal rocks (Hyndman & Wang 1993, 1995; Hyndman et al.1997; Oleskevich et al. 1999). It has been also proposed that earthquake ruptures initiated below this temperature bound may propagate with decreasing slip to where the temperature is 450 °C (Hyndman et al. 1997; Hyndman 2007). Due to the existence of subduction, the temperature is generally higher on the side of the overriding plate, which is roughly around 600 °C (Klingelhoefer et al.2010). The rapid subduction of cooler oceanic lithosphere will depress the isotherms on the megathrust interface, but this can account for only a 5 km deepening when subduction rate changes from 2 to 4 cm yr–1 (Klingelhoefer et al.2010). To further understand what is causing the factor of 2 difference between the Himalayian rupture depth and the Maule rupture depth, we follow Royden (1993) and developed a 1-D heat conduction model that accounts for the effect of erosion rate (Appendix). We assume that the surface is at a quasi-steady stage where the uplift from accretion and rebound is balanced by the effect of erosion. We then adopt the parameters in Royden (1993) who assumes the erosion rate to be 1 mm yr–1 and vary the remaining variables to calculate the corresponding temperature profiles (Fig. 6). From our analysis, both the heat production in the upper crust and the mantle temperature has little effect on depth of the 600 °C isotherm, while the erosion rate influences the temperature profile significantly. This can be understood by considering the effects of accretion or erosion/rebound; the hotter materials are directly brought up towards the surface, which is far more efficient than heat conduction. Recent studies show that in the region of the Gorkha earthquake, the surface erosion rate increases from about 0.5 mm yr–1 in the up-dip region to 3.5 mm yr–1 downdip (Whipple et al.2016), while for south-central Chile, this number is quite small (<0.3 mm yr–1) (Aguilar et al.2011; Carretier et al.2013). Fig. 6(a) shows that an average 1 mm yr–1 erosion rate is sufficient to bring the 600 °C isotherm from a depth of 40–17 km, which is consistent with our findings from spectral expansion inversion. Figure 6. View largeDownload slide Temperature versus depth for 1-D heat conduction model. (a) Models with different accretion/erosion rate. (b) Models with different upper crust heat generation rate. (c) Models with different mantle temperature. Figure 6. View largeDownload slide Temperature versus depth for 1-D heat conduction model. (a) Models with different accretion/erosion rate. (b) Models with different upper crust heat generation rate. (c) Models with different mantle temperature. 5.2 Correlations between forearc ridges and the downdip limit of the seismogenic zone Bassett & Watts (2015a,b) have developed a spectral averaging method of suppressing the large-amplitude, long-wavelength, trench-normal topographic and gravimetric expression of subduction zones. The global application of this technique revealed that in at-least five circum-Pacific subduction zones with a history of large megathrust earthquakes, the downdip limit of the seismogenic zone was approximately correlated with the location of a trench-parallel forearc ridge (TPFR). The longest TPFR parallels the coastline of Chile for >2000 km and is coincident with the Chilean coastal Cordillera. As shown in Fig. 7, the TPFR in southern Chile is well correlated with the maximum depth of earthquakes with Global Centroid Moment Tensors (GCMT; Dziewonski et al. 1981; Ekström et al.2012) consistent with megathrust slip, the maximum depth of co-seismic slip in the 2010 Mw 8.8 Maule earthquake (Tong et al.2010; Delouis et al.2010; this study), and the maximum depth of interseismic coupling as constrained by campaign GPS measurements (Moreno et al.2010). A similar correlation has been recognized in northern Chile (Béjar-Pizarro et al.2013). Figure 7. View largeDownload slide (a) Residual free-air gravity anomaly at the 2010 Chile Mw 8.8 earthquake region with grey dashed line denotes the trench axis and the black dashed line represents the TPFR. Grey pluses mark GCMT thrust events with Mw < 6.5. GCMT events with Mw ≥ 6.5 are plotted with beach ball, scaled for magnitude, coloured for depth. (b) Same as (a) with grey solid contours are the inverted slip pattern with each line representing 3 m slip. (c) Fraction of Plate Convergence from (Moreno et al.2010). The grey dashed line and black dashed line are the same as in (a) and (b). Figure 7. View largeDownload slide (a) Residual free-air gravity anomaly at the 2010 Chile Mw 8.8 earthquake region with grey dashed line denotes the trench axis and the black dashed line represents the TPFR. Grey pluses mark GCMT thrust events with Mw < 6.5. GCMT events with Mw ≥ 6.5 are plotted with beach ball, scaled for magnitude, coloured for depth. (b) Same as (a) with grey solid contours are the inverted slip pattern with each line representing 3 m slip. (c) Fraction of Plate Convergence from (Moreno et al.2010). The grey dashed line and black dashed line are the same as in (a) and (b). To investigate whether a similar relation is observed for the Gorka earthquake, we have applied this technique of spectral averaging to calculate residual topographic anomalies along the Himalayas. We have analysed the global SRTM15 topographic grid, which has a resolution of 15 arcsec (roughly 500 m). Residual topographic anomalies along the Himalayas are shown in Fig. 8. The Gorka earthquake occurred immediately up-dip of an elongated, strike-parallel region of positive residual topographic anomalies, which correlates remarkably well with the 3500 m elevation contour of the Himalayas (Avouac 2003; Stevens & Avouac 2015). This elevated band may be analogous to the TPFR previously observed in oceanic subduction zones. Black contours show the distribution of interseismic coupling as constrained by campaign GPS measurements (Stevens & Avouac 2015). It is intriguing that we observe a similar spatial association between the downdip transition from an interseismically locked to creeping megathrust, with the elevated residual topography of the high Himalayas. Figure 8. View largeDownload slide Residual topography for the 2015 Mw 7.8 Nepal earthquake region. The grey dashed line denotes the Main Frontal Thrust (MFT). The grey lines are the locking contours from (Stevens & Avouac 2015). The grey lines in the subplot at upper-right corner are 1-m slip contours. The bottom-left plot is the average elevation at this region. Figure 8. View largeDownload slide Residual topography for the 2015 Mw 7.8 Nepal earthquake region. The grey dashed line denotes the Main Frontal Thrust (MFT). The grey lines are the locking contours from (Stevens & Avouac 2015). The grey lines in the subplot at upper-right corner are 1-m slip contours. The bottom-left plot is the average elevation at this region. In subduction zones, the broad pattern of interseismic deformation can be modelled by placing an edge dislocation at the downdip end of the locked fault with its Burger's vector parallel to fault dip (Savage 1983). This model predicts a broad interseismic subsidence over the majority of the locked fault, with interseismic uplift located near the downdip limit of locking. Over many earthquake cycles, the incomplete recovery of interseismic elastic deformation by earthquakes and post-seismic slip may lead to permanent deformation of the forearc that spatially mirrors the sign of interseismic deformation. This model has been proposed to explain long-term coastal uplift in oceanic subduction zones in Cascadia (Kelsey et al. 1994) and Mexico (Ramirez-Herrera et al.2004) and we suggest unrecovered interseismic deformation may also be contributing to the elevation and growth of the high Himalaya in eastern Nepal, which is in agreement with previous studies (Meade 2010). 6 CONCLUSIONS In this study, we developed and tested a data-driven spectral expansion approach that possess the advantages of (1) mathematical robustness, (2) smoothness and stableness without regularization, (3) capability of adopting complex fault surface, (4) giving real uncertainty estimates and (5) providing resolving power analysis. With this approach, we confirmed that the downdip ruptures of the 2010 Mw 8.8 Maule earthquake and the 2015 Mw 7.8 Gorkha earthquake are well resolved by the data. The Maule earthquake rupture terminates around 45 km depth and the Gorkha earthquake rupture stops shallower than 20 km. Both these depths are well correlated with the location of trench-parallel forearc ridges, which may reflect unrecovered interseismic elastic deformation. To understand the difference in maximum rupture depth between the two cases, we constructed 1-D thermal conduction models with different erosion rate, upper-crust heat production and mantle temperature. By comparing different models, we conclude that the erosion rate has a major influence on the downdip rupture limit. ACKNOWLEDGEMENTS We want to thank the two reviewers for their valuable suggestions. This study was funded by the NASA Earth Surface and Interior Program (NNX16AK93G) and the Southern California Earthquake Center (SCEC). SCEC is funded by the NSF cooperative Agreement EAR-1033462 and USGS Cooperative Agreement G12AC20038. Special thanks to Prof Catherine Constable for sharing with us her knowledge on Geophysical Inverse Theory. We also want to thank the authors from Tong et al. (2010), Lindsey et al. (2015) and Galetzka et al. (2015) for sharing their data to us. REFERENCES Ader T. et al.  , 2012. Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: implications for seismic hazard, J. geophys. Res. , 117, B04403, doi:10.1029/2011JB009071. Google Scholar CrossRef Search ADS   Aguilar G. et al.  , 2011. Variability in erosion rates related to the state of landscape transience in the semi‐arid Chilean Andes, Earth Surface Processes and Landforms , 36( 13), 1736– 1748. Ammon C.J. et al.  , 2005. Rupture process of the 2004 Sumatra-Andaman earthquake, Science , 308( 5725), 1133– 1139. Google Scholar CrossRef Search ADS PubMed  Avouac J.P., 2003. Mountain building, erosion, and the seismic cycle in the Nepal Himalaya, Adv. Geophys. , 46, 1– 80. Google Scholar CrossRef Search ADS   Bassett D., Watts A.B., 2015a. Gravity anomalies, crustal structure, and seismicity at subduction zones: 1. Seafloor roughness and subducting relief, Geochem. Geophys. Geosyst. , 16( 5), 1508– 1540. Google Scholar CrossRef Search ADS   Bassett D., Watts A.B., 2015b. Gravity anomalies, crustal structure, and seismicity at subduction zones: 2. Interrelationships between fore‐arc structure and seismogenic behavior, Geochem. Geophys. Geosyst. , 16( 5), 1541– 1576. Google Scholar CrossRef Search ADS   Béjar-Pizarro M., Socquet A., Armijo R., Carrizo D., Genrich J., Simons M., 2013. Andean structural control on interseismic coupling in the North Chile subduction zone, Nat. Geosci. , 6( 6), 462– 467. Google Scholar CrossRef Search ADS   Carretier S. et al.  , 2013. Slope and climate variability control of erosion in the Andes of central Chile, Geology , 41( 2), 195– 198. Google Scholar CrossRef Search ADS   Chlieh M. et al.  , 2007. Coseismic slip and afterslip of the great Mw 9.15 Sumatra–Andaman earthquake of 2004, Bull. seism. Soc. Am. , 97( 1A), S152– S173. Google Scholar CrossRef Search ADS   Cohen S.C., 1999. Numerical models of crustal deformation in seismic zones, Adv. Geophys. , 41, 133– 231. Google Scholar CrossRef Search ADS   Delaunay B., 1934. Sur la sphere vide. Izv. Akad. Nauk SSSR, Otdelenie Matematicheskii i Estestvennyka Nauk , 7( 793-800), 1– 2. Delouis B., Nocquet J.-M., Vallée M., 2010. Slip distribution of the February 27, 2010 Mw = 8.8 Maule earthquake, central Chile, from static and high‐rate GPS, InSAR, and broadband teleseismic data, Geophys. Res. Lett. , 37, L17305, doi:10.1029/2010GL043899. Google Scholar CrossRef Search ADS   Dziewonski A.M., Chou T.A., Woodhouse J.H., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, J. geophys. Res. , 86( B4), 2825– 2852. Google Scholar CrossRef Search ADS   Ekström G., Nettles M., Dziewoński A.M., 2012. The global CMT project 2004–2010: Centroid-moment tensors for 13,017 earthquakes, Phys. Earth planet. Inter. , 200, 1– 9. Google Scholar CrossRef Search ADS   Fialko Y., 2004. Probing the mechanical properties of seismically active crust with space geodesy: study of the coseismic deformation due to the 1992 Mw7.3 Landers (southern California) earthquake, J. geophys. Res. , 109, B03307, doi:10.1029/2003JB002756. Fritz H.M. et al.  , 2011. Field survey of the 27 February 2010 Chile tsunami, Pure appl. Geophys. , 168( 11), 1989– 2010. Google Scholar CrossRef Search ADS   Galetzka J. et al.  , 2015. Slip pulse and resonance of the Kathmandu basin during the 2015 Gorkha earthquake, Nepal, Science , 349( 6252), 1091– 1095. Google Scholar CrossRef Search ADS PubMed  Gilbert F., 1971. Ranking and winnowing gross Earth data for inversion and resolution, Geophys. J. Int. , 23( 1), 125– 128. Google Scholar CrossRef Search ADS   Hyndman R.D., 2007. The Seismogenic Zone of Subduction Thrust Faults , Columbia Univ. Press, pp. 15– 40. Hyndman R.D., Wang K., 1993. Thermal constraints on the zone of major thrust earthquake failure: the Cascadia subduction zone, J. geophys. Res. , 98( B2), 2039– 2060. Google Scholar CrossRef Search ADS   Hyndman R.D., Yamano M., Oleskevich D.A., 1997. The seismogenic zone of subduction thrust faults, Island Arc , 6( 3), 244– 260. Google Scholar CrossRef Search ADS   Hyndman R.D., Wang K., Yamano M., 1995. Thermal constraints on the seismogenic portion of the southwestern Japan subduction thrust, J. geophys. Res. , 100( B8), 15 373–15 392. Google Scholar CrossRef Search ADS   Jónsson S., Zebker H., Segall P., Amelung F., 2002. Fault slip distribution of the 1999 Mw 7.1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements, Bull. seism. Soc. Am. , 92( 4), 1377– 1389. Google Scholar CrossRef Search ADS   Kelsey H.M., Engebretson D.C., Mitchell C.E., Ticknor R.L., 1994. Topographic form of the Coast Ranges of the Cascadia margin in relation to coastal uplift rates and plate subduction, J. geophys. Res. , 99( B6), 12 245–12 255. Google Scholar CrossRef Search ADS   Klingelhoefer F. et al.  , 2010. Limits of the seismogenic zone in the epicentral region of the 26 December 2004 great Sumatra‐Andaman earthquake: results from seismic refraction and wide‐angle reflection surveys and thermal modeling, J. geophys. Res. , 115, B01304, doi:10.1029/2009JB006569. Google Scholar CrossRef Search ADS   Lay T. et al.  , 2005. The great Sumatra-Andaman earthquake of 26 December 2004, Science , 308( 5725), 1127– 1133. Google Scholar CrossRef Search ADS PubMed  Lindsey E.O., Natsuaki R., Xu X., Shimada M., Hashimoto M., Melgar D., Sandwell D.T., 2015. Line‐of‐sight displacement from ALOS‐2 interferometry: Mw 7.8 Gorkha Earthquake and Mw 7.3 aftershock, Geophys. Res. Lett. , 42( 16), 6655– 6661. Google Scholar CrossRef Search ADS   Lorito S. et al., 2011. Limited overlap between the seismic gap and coseismic slip of the great 2010 Chile earthquake, Nat. Geosci. , 4( 3), 173– 177. Google Scholar CrossRef Search ADS   Meade B.J., 2010. The signature of an unbalanced earthquake cycle in Himalayan topography?, Geology , 38( 11), 987– 990. Google Scholar CrossRef Search ADS   Métois M., Socquet A., Vigny C., 2012. Interseismic coupling, segmentation and mechanical behavior of the central Chile subduction zone, J. geophys. Res. , 117, B03406, doi:10.1029/2011JB008736. Google Scholar CrossRef Search ADS   Minson S.E. et al.  , 2014. Bayesian inversion for finite fault earthquake source models—II: The 2011 great Tohoku-oki, Japan earthquake, Geophys. J. Int. , 198( 2), 922– 940. Google Scholar CrossRef Search ADS   Moreno M., Rosenau M., Oncken O., 2010. 2010 Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone, Nature , 467( 7312), 198– 202. Google Scholar CrossRef Search ADS PubMed  National Earthquake Information Center (NEIC), 2010. Magnitude 8.8-Offshore Maule, Chile, U.S. Geol. Surv., Denver, CO. Available at: http://earthquake.usgs.gov/earthquakes/eventpage/official20100227063411530_30#executive. National Earthquake Information Center (NEIC), 2015. Magnitude 8.8-Offshore Maule, Chile, U.S. Geol. Surv., Denver, CO. Available at: http://earthquake.usgs.gov/earthquakes/eventpage/official20100227063411530_30#executive. Okada Y., 1985. Surface deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am. , 75( 4), 1135– 1154. Oleskevich D.A., Hyndman R.D., Wang K., 1999. The updip and downdip limits to great subduction earthquakes: Thermal and structural models of Cascadia, south Alaska, SW Japan, and Chile, J. geophys. Res. , 104( B7), 14 965–14 991. Google Scholar CrossRef Search ADS   Parker R.L., 1977. Understanding inverse theory, Ann. Rev. Earth planet. Sci. , 5, 35– 64. Google Scholar CrossRef Search ADS   Parker R.L., 1994. Geophysical Inverse Theory , Princeton Univ. Press. Ramırez-Herrera M.T., Kostoglodov V., Urrutia-Fucugauchi J., 2004. Holocene-emerged notches and tectonic uplift along the Jalisco coast, Southwest Mexico, Geomorphology , 58( 1), 291– 304. Google Scholar CrossRef Search ADS   Royden L.H., 1993. The steady state thermal structure of eroding orogenic belts and accretionary prisms, J. geophys. Res. , 98, 4487– 4507. Google Scholar CrossRef Search ADS   Savage J.C., 1983. A dislocation model of strain accumulation and release at a subduction zone, J. geophys. Res. , 88( B6), 4984– 4996. Google Scholar CrossRef Search ADS   Simons M., Fialko Y., Rivera L., 2002. Coseismic deformation from the 1999 Mw 7.1 Hector Mine, California, earthquake as inferred from InSAR and GPS observations, Bull. seism. Soc. Am. , 92( 4), 1390– 1402. Google Scholar CrossRef Search ADS   Simons M. et al.  , 2011. The 2011 magnitude 9.0 Tohoku-Oki earthquake: mosaicking the megathrust from seconds to centuries, Science , 332( 6036), 1421– 1425. Google Scholar CrossRef Search ADS PubMed  Smith W.H.F., Wessel P., 1990. Gridding with continuous curvature splines in tension, Geophysics , 55( 3), 293– 305. Google Scholar CrossRef Search ADS   Steketee J.A., 1958. On Volterra's dislocations in a semi-infinite elastic medium, Can. J. Phys. , 36( 2), 192– 205. Google Scholar CrossRef Search ADS   Stevens V.L., Avouac J.P., 2015. Interseismic coupling on the main Himalayan thrust, Geophys. Res. Lett. , 42( 14), 5828– 5837. Google Scholar CrossRef Search ADS   Tarantola A., 2005. Inverse Problem Theory and Methods for Model Parameter Estimation , Chap. 1, pp. 24– 40, SIAM. Tong X. et al., 2010. The 2010 Maule, Chile earthquake: downdip rupture limit revealed by space geodesy, Geophys. Res. Lett. , 37, L24311, doi:10.1029/2010GL045805. Google Scholar CrossRef Search ADS   Wang K., Fialko Y., 2015. Slip model of the 2015 Mw 7.8 Gorkha (Nepal) earthquake from inversions of ALOS‐2 and GPS data, Geophys. Res. Lett. , 42( 18), 7452– 7458. Google Scholar CrossRef Search ADS   Whipple K.X., Shirzaei M., Hodges K.V., Arrowsmith J.R., 2016. Active shortening within the Himalayan orogenic wedge implied by the 2015 Gorkha earthquake, Nat. Geosci. , 9( 9), 711– 716. Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. (a) Comparison between surface deformation derived from the input slip (red dots, noise added) and inverted slip (7 kernels, black line) corresponding to Fig. 2(c). (b) Misfit function versus the number of kernels used for Fig. 2(c). (c) Comparison between surface deformation corresponding to Fig. 2(d) (11 kernels). (d) Misfit function plot versus number of kernels used for Fig. 2(d). Figure S2. (a) An example 2-D kernel for the 2010 Maule Earthquake, with black circles denoting the peaks and troughs of the kernel. (b) Delaunay triangulation for the nodes (circles) in (a) (adding 8 nodes at corners and middle of sides, value of each node is the average distance from the node to its adjacent nodes), with colour representing the average from the triangle's three nodes. (c) Map generated applying a surface technique and filtering to the nodes in (b). Figure S3. Plots (a), (b) and (c) are the subsampled data, forwarded model, and misfit for the simulated noisy ascending data. Plots (d), (e) and (f) are the subsampled data, forwarded model, and misfit for the simulated noisy descending data. (g) Uncertainty map for dip slip in Fig. 3(d). (h) Uncertainty map for strike slip in Fig. 3(f). (i) Eigenvalues versus order. Figure S4. Plots (a), (b) and (c) are the subsampled data, forwarded model, and misfit for the ascending data for the 2010 Mw 8.8 Maule Earthquake. Plots (d), (e) and (f) are the subsampled data, forwarded model, and misfit for the descending data. (g) Uncertainty map for dip slip in Fig. 4a. (h) Uncertainty map for strike slip in Fig. 4b. (i) Eigenvalues versus order. Figure S5. Plots (a), (b) and (c) are the subsampled data, forwarded model, and misfit for the data from ALOS-2 ascending track 157 for the 2015 Mw 7.8 Nepal earthquake. Plots (d), (e) and (f) are the subsampled data, forwarded model, and misfit for the data from ALOS-2 descending track 47. Plots (g), (h) and (i) are the subsampled data, forwarded model, and misfit for the data from ALOS-2 descending track 48. (j) Uncertainty map for dip slip in Fig. 5a. (k) is the uncertainty map for strike slip in Fig. 5b. (l) Eigenvalues versus order. Figure S6. Some examples of kernels for the 3-D test. The upper row shows the dip slip kernels and the lower row shows the strike slip kernels. The wavelength of the kernel decreases as the order increases. Figure S7. 3-D strike-slip test case. (a) Input strikeslip with the black dashed line denoting the fault trace. (b) Cumulative strike slip $$p( z ) = \mathop \int \nolimits_0^L m( {x,z} ){\rm{d}}x$$ versus depthwith black dashed line being the input strike slip, blue line representing the inverted strike slip, and magenta dash-dot lines denoting the uncertainty. (c) Misfit versus order with the red dot being the order selected. (d) Inverted dip slip and (f) inverted strike slip. Resolving wavelength for dip slip (e) and strike slip (g) constructed by the algorithm shown in Fig. S2. Figure S8. 3-D mixture of dip and strike sliptest case. (a) Input strike and dip slip (same amplitude) with the black dashed line denoting the fault trace. (b) Cumulative slip $$p( z ) = \mathop \int \nolimits_0^L m( {x,z} ){\rm{d}}x$$ versus depth plot with black dashed line being the input dip/strike slip, blue line representing the inverted dip slip and cyan dash-dot lines denoting the uncertainty, red line representing the inverted strike slip, and magenta dash-dot lines denoting the uncertainty. (c) Misfit versus order with the red dot being the order selected. (d) Inverted dip slip and (f) inverted strike slip. Resolving wavelength for dip slip (e) and strike slip (g) constructed by the algorithm shown in Fig. S2. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX: APPENDIX: DERIVATION OF THE 1-D HEAT CONDUCTION MODEL WITH ACCRETION/EROSION Following (Royden 1993), we simplified and modified the heater conduction equation as follows, assuming the accretion and erosion rates are in balance.   \begin{equation*} \frac{{{\rm d^2}T}}{{{\rm{d}}{z^2}}} + \frac{v}{\alpha }\frac{{{\rm{d}}T}}{{{\rm{d}}z}} + \frac{{{A_0}}}{K}\ {{\rm{e}}^{ - z/h}} = \ 0, \end{equation*} where T is the temperature, z is the depth, v is the accretion rate, α is the thermal diffusivity, A0 is the heat production rate at the surface, K is the thermal conductivity and h is the characteristic decaying depth for heat production. We simplify the equation a little more by substituting some of the variables as follows:   \begin{equation*} \frac{{{{\rm{d}}^2}T}}{{{\rm{d}}{z^2}}} + a\frac{{{\rm{d}}T}}{{{\rm{d}}z}} + b\ {{\rm{e}}^{ - cz}} = \ 0 \end{equation*} Then further simplify the equation form a 2nd order ODE to a 1st order ODE by letting dT/dz = y. Then we have the equation as   \begin{equation*} \frac{{{\rm{d}}y}}{{{\rm{d}}z}} + ay + b\ {{\rm{e}}^{ - cz}} = \ 0 \end{equation*} As for the type of 1st order ODE like y' + P(x)y = Q(x), the general solution is   \begin{equation*} y\ = {e^{ - \mathop \int \nolimits_{{x_0}}^x P\left( t \right){\rm{d}}t}}\ \int Q\left( x \right){{\rm{e}}^{\mathop \int \nolimits_{{x_0}}^x P\left( t \right){\rm{d}}t}}{\rm{d}}x. \end{equation*} Thus, we have the solution for y as   \begin{equation*} y\ = \ - \frac{b}{{a - c}}\left( {{e^{ - cz}} + {C_1}{e^{ - az}}} \right) \end{equation*} Substitute y with dT/dz and solve the ODE, we’ll get   \begin{equation*} T\ = \frac{b}{{c\left( {a - c} \right)}}\ {e^{ - cz}} + {C_1}\frac{b}{{a\left( {a - c} \right)}}{e^{ - az}} + {C_2} \end{equation*} Consider the boundary condition   \begin{equation*} z\ = \ 0,\quad \ T\ = {T_0}{\rm{\ }} \end{equation*}   \begin{equation*} z\ = \ {\rm{H}},\quad \ T\ = {T_m}\ \end{equation*} We can get the coefficients   \begin{equation*} \ {C_1} = \ - \left( {{T_m}\frac{{a\left( {a - c} \right)}}{b} + \frac{a}{c}\left( {1 - {e^{ - cH}}} \right)} \right)\bigg/\left( {1 - {e^{ - aH}}} \right) \end{equation*}   \begin{equation*} \ {C_2} = \left( {{T_m} + \frac{b}{{c\left( {a - c} \right)}}\left( {{e^{ - aH}} - {e^{ - cH}}} \right)} \right)\ \bigg/\left( {1 - {e^{ - aH}}} \right) \end{equation*} © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# A spectral expansion approach for geodetic slip inversion: implications for the downdip rupture limits of oceanic and continental megathrust earthquakes

, Volume 212 (1) – Jan 1, 2018
12 pages

/lp/ou_press/a-spectral-expansion-approach-for-geodetic-slip-inversion-implications-ZkrV5K9Lk4
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx408
Publisher site
See Article on Publisher Site

### Abstract

SUMMARY We have developed a data-driven spectral expansion inversion method to place bounds on the downdip rupture depth of large megathrust earthquakes having good InSAR and GPS coverage. This inverse theory approach is used to establish the set of models that are consistent with the observations. In addition, the inverse theory method demonstrates that the spatial resolution of the slip models depends on two factors, the spatial coverage and accuracy of the surface deformation measurements, and the slip depth. Application of this method to the 2010 Mw 8.8 Maule Earthquake shows a slip maximum at 19 km depth tapering to zero at ∼40 km depth. In contrast, the continent–continent megathrust earthquakes of the Himalayas, for example 2015 Mw 7.8 Gorkha Earthquake, shows a slip maximum at 9 km depth tapering to zero at ∼18 km depth. The main question is why is the maximum slip depth of the continental megathrust earthquake only 50 per cent of that observed in oceanic megathrust earthquakes. To understand this difference, we have developed a simple 1-D heat conduction model that includes the effects of uplift and surface erosion. The relatively low erosion rates above the ocean megathrust results in a geotherm where the 450–600 °C transition is centred at ∼40 km depth. In contrast, the relatively high average erosion rates in the Himalayas of ∼1 mm yr–1 results in a geotherm where the 450–600 °C transition is centred at ∼20 km. Based on these new observations and models, we suggest that the effect of erosion rate on temperature explains the difference in the maximum depth of the seismogenic zone between Chile and the Himalayas. Satellite geodesy, Inverse theory, Intra-plate processes 1 INTRODUCTION Over the past few decades, advances in space geodesy such as Interferometric Synthetic Aperture Radar (InSAR) and global positioning system (GPS) have facilitated detailed mapping of surface deformation. These high quality and regional datasets have enabled scientists to build 3-D deformation fields for many large earthquakes. Usually, the fault geometry can be reasonably approximated from aftershocks, geophysical imaging and geological field studies. Elastic models that relate slip at depth to surface deformation can then be constructed by discretizing this fault into small patches before estimating the slip on each patch. While this method is straightforward, it hardly provides an estimate of the spatial resolution nor unbiased uncertainties associated with the slip model. It is well known that the equations which relate slip at depth to surface deformation contain an upward continuation term (Steketee 1958), which exponentially attenuates the rupture signal at wavelengths smaller than the depth, therefore making this inverse problem inherently non-unique (Parker 1994). This is true even if the data coverage is complete and the data have no errors. The non-uniqueness becomes evident in parametrized inversions when the discretized patch size is smaller than the depth; the inversions are poorly conditioned resulting in wild oscillations in slip distribution. In order to stabilize the solution, smoothness regularization combined with a non-negativity constraint is applied to the inversion, which acts as a strong prior knowledge to the model (Jonsson et al. 2002; Simons et al. 2002; Fialko 2004; Tarantola 2005). In some studies, a checkerboard test is used to illustrate the resolving power of the data (Fialko 2004; Tong et al. 2010). Tong et al. (2010) performed this test for the Mw 8.8 2010 Maule earthquake and found the model can only recover a 20-km checker up to about 35 km depth, but can recover a 40-km checker to 60 km depth. These results are consistent with the inverse relation between model resolvability and rupture depth. For several recent megathrust earthquakes, surface deformation has been well recorded and many studies have incorporated these data into earthquake rupture inversions using the standard parameter estimation approach. Inversions for the 2010 Mw 8.8 Maule earthquake (oceanic plate subducting) suggest a tapering of slip towards zero at ∼40 km depth (Tong et al. 2010; Lorito et al. 2011). This depth is similar to many other oceanic subduction slip inversions [e.g. the 2011 Mw 9.0 Tohoku-oki (Simons et al. 2011; Minson et al. 2014); 2004 Mw 9.2 Sumatra–Andaman (Ammon et al. 2005; Lay et al. 2005; Chlieh et al. 2007)]. By contrast, inversions for the 2015 Mw 7.8 Nepal earthquake (continental plate subducting) suggest a tapering of slip towards zero at less than 20 km (Galetzka et al. 2015; Wang & Fialko 2015). This depth is consistent with the maximum depth of strong interseismic coupling resolved by geodetic models along the Himalayan front (Stevens & Avouac 2015). This difference raises two questions. First, is the difference in downdip rupture limit well resolved by the observations? If it is well resolved, what physical mechanisms may explain such a large difference? This requires a real uncertainty estimate as well as an analysis of resolving power for inversion models. To address the first question, we adopt the approach of linear inverse theory (Gilbert 1971; Parker 1977). As noted by Parker (1977): ‘The quality that distinguishes inverse theory from the parameter estimation problem of statistics is the unknowns are functions, not merely a handful of real numbers. This means that the solution contains in principle an infinite number of variables, and therefore with real data the problem is as underdetermined as it can be.’ The approach begins with the well-known Green's functions relating slip at depth on a prescribed fault plane to surface deformation. Following the spectral expansion method, we then construct a set of up to N orthonormal slip functions that span the set of N surface observations. The solution and its uncertainty are constructed from a finite number M ≤ N of these functions, or kernels, that match the data to within their uncertainty. The low order kernels are generally smoothest and well constrained by the data while the higher order kernels have more oscillations and are more poorly constrained. The solution is constructed from the well-resolved kernels. We use this approach to provide bounds on the depth distribution of slip during the Mw 8.8 Maule (oceanic) and Mw 7.8 Gorka (continental) megathrust areas, and show that the factor-of-two difference in slip depth is indeed well resolved by the geodetic data. To address the second question, the pronounced difference in the slip depth, we develop a simple 1-D thermal model that includes the effects of vertical advection of heat due to erosion-induced rebound (Royden 1993). We find that an erosion rate of 1 mm yr–1 decreases the depth of the 600 °C isotherm by a factor of two, which can explain the relatively shallow rupture depth of the Gorka earthquake. 2 THE SPECTRAL EXPANSION APPROACH Here we give a brief summary of the spectral expansion theory from Parker (1977, 1994) as applied to static slip inversion. The theory is originally from Gilbert (1971). The main difference between our formulation and the original formulation is that our model is 2-D vector slip (strike and dip) as a function of the two dimensions on the fault surface m(x). Note that the fault does not need to be a planar surface. The data are surface deformation sampled at a finite number of points on the surface of the Earth. GPS can provide 3 components of deformation while InSAR provides 1 or 2 line of sight components. Each component of each data point is di and has an uncertainty σi. We use the elastic half-space Green's function (Okada 1985) to relate vector slip at depth to each component of surface deformation. We note that Green's functions for a more complicated elastic structure could be used. The convolution of the model and the Green's function is given by   $${d_i} = \mathop \int \nolimits_{\rm{\Sigma }}^{} {G_i}\left( {\boldsymbol{x}} \right)\ {\boldsymbol{m}}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }},\$$ (1)where Σ represents integration over the entire fault plane. After dividing both sides of eq. (1) by the uncertainty σi the equation can be rewritten as   $$d{'_i} = \mathop \int \nolimits_{\rm{\Sigma }}^{} {G'_i}\left( {\boldsymbol{x}} \right){\boldsymbol{m}}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }}.\$$ (2) Now the data d'i as well as the Green's functions G'i are dimensionless and the data have unit standard deviation. Following Parker (1977) we form the Gram matrix $${{\boldsymbol \Gamma }}$$  $${{\rm{\Gamma }}_{ij}} = \mathop \int \nolimits_{\rm{\Sigma }}^{} {G'_i}\left( {\boldsymbol{x}} \right){G'_j}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }}.\$$ (3) Thus, $${{\boldsymbol \Gamma }}$$ is positive-definite and symmetric, and can be diagonalized with an orthonormal matrix O of eigenvectors   \begin{eqnarray} {{{\bf O}}^{\rm{T}}} {\boldsymbol \Gamma} {\bf O}\ &=& {\boldsymbol{\ }}{{\boldsymbol \Lambda }},\quad {\rm{where\ }}{{\boldsymbol \Lambda }}\ = {\rm{\ diag}} \{ {{{\rm{\lambda }}_1},{{\rm{\lambda }}_2}, \ldots ,{{\rm{\lambda }}_{N\ }}} \}\nonumber\\ &&\qquad{\rm{\ and\ }}{{\rm{\lambda }}_1} \ge {{\rm{\lambda }}_2} \ge {{\rm{\lambda }}_3} \ge \ldots {{\rm{\lambda }}_{N\ }} > 0, \end{eqnarray} (4)where λi are the eigenvalues; the full set of eigenvalues is called the spectrum. It's not difficult to show that   \begin{eqnarray} {\lambda _i} \!=\! {O_i}^{\rm{T}}\ {{\boldsymbol \Gamma }}\ {O_i} &=& \mathop \sum \limits_j {O_{ji}}\mathop \int \nolimits_{\rm{\Sigma }}^{} {{\bf G}}{{{\bf G}}^{\rm{T}}}{\rm{d\Sigma }}{O_{ij}}\nonumber\\ &=& \mathop \int \nolimits_{\rm{\Sigma }}^{} \mathop \sum \limits_j {O_{ji}}{{\bf G}}{{{\bf G}}^{\rm{T}}}{O_{ij}}{\rm{d\Sigma }},\ \end{eqnarray} (5)where G is a column vector of normalized Green's functions and the Gram matrix $${{\boldsymbol \Gamma }}$$ can be treated as a dyad matrix. We rewrite the equation as following:   $${{\boldsymbol \Lambda }}\ = \mathop \int \nolimits_{\rm{\Sigma }}^{} {{{\bf O}}^{\rm{T}}}{{\bf G}}{{{\bf G}}^{\rm{T}}}{{\bf O}}\, {\rm{d\Sigma }}.\$$ (6) Now consider the kernel functions defined by   $${\psi _i}\left( {\boldsymbol{x}} \right) = {\lambda _i}^{ - \frac{1}{2}}\ \mathop \sum \limits_j {O_{ji}}{G'_j}\left( {\boldsymbol{x}} \right).$$ (7) It can be easily verified from the eq. (5) that ψi(x) are an orthonormal set   $$\mathop \int \nolimits_{\rm{\Sigma }}^{} {\psi _i}\left( {\boldsymbol{x}} \right){\psi _j}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }} = {\delta _{ij}}.\$$ (8) Therefore, we consider the expansion of the model m(x) in terms of these kernels   $${\boldsymbol {m}}( x) = \mathop \sum \limits_j {a_i}\ {\psi _i}\left( {\boldsymbol{x}} \right) + {\psi _*},$$ (9)where ψ* is the annihilator (i.e. models that produce no surface deformation at the positions of the data points). The coefficients of the expansion can be derived   \begin{equation*} {a_i} = \mathop \int \nolimits_{\rm{\Sigma }}^{} {\psi _i}\left( {\boldsymbol{x}} \right){\boldsymbol{m}}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }} \end{equation*}   \begin{equation*} \quad= \mathop \int \limits_{\rm{\Sigma }}^{} {\lambda _i}^{ - \frac{1}{2}}\mathop \sum \limits_j {O_{ji}}{G'_j}\left( {\boldsymbol{x}} \right){\boldsymbol{m}}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }}\ \end{equation*}   \begin{equation*} \quad= {\lambda _i}^{ - \frac{1}{2}}\ \mathop \sum \limits_j {O_{ji}}\mathop \int \nolimits_{\rm{\Sigma }}^{} {G'_i}\left( {\boldsymbol{x}} \right){\boldsymbol{m}}\left( {\boldsymbol{x}} \right){\rm{d\Sigma }} \end{equation*}   $$\quad= {\lambda _i}^{ - \frac{1}{2}}\ \mathop \sum \limits_j {O_{ji}}d{'_j}.$$ (10) Since all the d'j have unit variance and the matrix O is orthonormal, the error estimates of ai are simply λi − 1/2 and ai are statistically independent. When the inverse problem is handled with this approach, it is easy to isolate the parts that are better determined, that is those have smaller error or larger eigenvalue. 3 APPLICATION TO 2-D AND 3-D TEST CASES We first test this approach using an infinitely-long vertical strike-slip fault. Two input slip distributions are used for testing. The first is uniform slip (1 m) with depth between the surface and 12 km (Fig. 1c, black dashed). The second is a Gaussian-shaped slip centred at 8 km depth (Fig. 1, black dashed). We then use the 2-D half-space Green's function (Cohen 1999) to forward generate surface observations from –50 to 50 km distance with 200 m sampling increment and add a random noise of 2 per cent of the maximum shear deformation (Fialko 2004). Next, we normalize the Green's function, compute the Gram matrix analytically and decompose it to get the eigenvalues (Fig. 1a) as well as the orthonormal kernels (Fig. 1b). Note that the orthonormal kernels have oscillations as a function of depth with more oscillations for higher kernel number. After this step, we calculate the coefficient for each kernel and combine them to get the inverted slip (Fig. 1c, blue line). The uncertainty of the inverted slip can be acquired by propagating the error of each coefficient to the model (Fig. 1c, magenta dash--dotted line). Figure 1. View large Download slide 2-D infinite fault test case. (a) Eigenvalues versus order, with the blue circle denoting the order selected for the step case (c) and blue star for the Gaussian case (d) when using noisy data. (b) Some examples of orthonormal kernels versus depth. (c) Recovered slip versus depth plot for input slip being a step function and (d) recovered slip versus depth plot for input slip being a Gaussian function. For (c) and (d), the black dashed lines represent the input slip, blue lines are slip inverted using noisy data with magenta dash--dotted lines being the uncertainty bounds, light blue lines are the slip inverted using noise-free data, transparent red and green lines are inverted slip using noise-free data with more kernels. The forward model and fitting can be found in Fig. S1. Figure 1. View large Download slide 2-D infinite fault test case. (a) Eigenvalues versus order, with the blue circle denoting the order selected for the step case (c) and blue star for the Gaussian case (d) when using noisy data. (b) Some examples of orthonormal kernels versus depth. (c) Recovered slip versus depth plot for input slip being a step function and (d) recovered slip versus depth plot for input slip being a Gaussian function. For (c) and (d), the black dashed lines represent the input slip, blue lines are slip inverted using noisy data with magenta dash--dotted lines being the uncertainty bounds, light blue lines are the slip inverted using noise-free data, transparent red and green lines are inverted slip using noise-free data with more kernels. The forward model and fitting can be found in Fig. S1. An additional test is done to explore the recoverability by inverting the synthetic data with no noise. The light blue line in Fig. 1(c) shows the solution using the exact same number of kernels, while the transparent red and green lines are solutions using more. It's obvious that with more kernels, the solution will get closer to the step function we input, despite the fact that due to the upward continuation, a full recovery of the step function is not expected. Also, the rapid reduction in the eigenvalue causes the error from the data to be exponentially magnified when it comes to higher order. To further examine this approach, we set a Gaussian input slip centring 8 km depth with the maximum slip being 1 m (Fig. 1d). The conclusion remains the same. In both cases, the surface deformation is well fit with the selected kernels determined from the misfit function $${\chi ^2} = \frac{1}{N}\ \mathop \sum \nolimits_{i\ = \ 1}^N {(\frac{{{d_i} - d_i^{\rm{pred}}}}{{{\sigma _i}}})^2}\$$ (Fig. S1). Note that from the theory (Section 2), if the same fault geometry is determined, once the locations of observations are set, every computation step will be the same until the calculation of coefficients for the kernels. In other words, both cases shown in Figs 1(c) and (d) are just different recompositions of the kernels shown in Fig. 1(b), based on different surface observations. After considering the 2-D test, we expanded the formulation to a 3-D test case. We use a 350-km-long 250-km-wide planar fault with N30°E strike and 15° dip towards east. The input dip-slip is a Gaussian function shown in Fig. 2(a), with the corresponding moment magnitude being 8.0. Then we generate a forward model using the half-space Green's function (Okada 1985) and project the vector surface deformation into descending and ascending line-of-sight (LOS) directions to simulate the InSAR observations. Random noise (2 per cent of maximum amplitude) (Fialko 2004) is added and then the data are down sampled with a quad-tree algorithm (Jonsson et al.2002) based on the curvature of the deformation field. Due to the complexity of the 3-D half-space Green's function (Okada 1985), the computation of the Gram matrix is done with numerical integration for both strike and dip slip. We then decompose the Gram matrix to recover the eigenvalues and orthonormal kernels (Figs S3i and S6). We found the 3-D case has a much slower fall-off in misfit with increasing number of kernel functions than the 2-D case. To determine how many kernels to use, we examine the misfit function (Fig. 2c) and chose the smallest number with reasonable misfit reduction (e.g. 100 kernels), in which case, more kernels will result in an explosion in model's uncertainty but no significant improvement on data fitting. While the along-strike cumulative dip slip shows a good fit to the input (Fig. 2b) with the recovered moment magnitude being 7.99, a small component of strike slip was also ‘recovered’, potentially due to incomplete coverage from down-sampled data and the noise we added. Besides, we also discovered some short wavelength undulations in the inverted slip, which we believe is due to the nature of this decomposition method. Figure 2. View large Download slide 3-D test case. (a) Input dip slip with positive defined as hanging wall moving up-dip. The black dashed line is the fault trace. (b) Cumulative dip slip $$p( z )\ = \mathop \int \nolimits_0^L m( {x,z} ){\rm{d}}x\ \$$ (Simons et al. 2002) versus depth plot with black dashed line being the input dip slip, blue line representing the inverted dip slip and magenta dash-dotted lines denoting the uncertainty. (c) Misfit versus number of kernels, with the red dot being the order selected. (d) Inverted dip slip and (f) inverted strike slip with hanging wall moving north-eastward defined as positive. Resolving wavelength for dip slip (e) and strike slip (g) constructed by the algorithm shown in Fig. S2. Data fitting and eigenvalues can be found in Fig. S3. Figure 2. View large Download slide 3-D test case. (a) Input dip slip with positive defined as hanging wall moving up-dip. The black dashed line is the fault trace. (b) Cumulative dip slip $$p( z )\ = \mathop \int \nolimits_0^L m( {x,z} ){\rm{d}}x\ \$$ (Simons et al. 2002) versus depth plot with black dashed line being the input dip slip, blue line representing the inverted dip slip and magenta dash-dotted lines denoting the uncertainty. (c) Misfit versus number of kernels, with the red dot being the order selected. (d) Inverted dip slip and (f) inverted strike slip with hanging wall moving north-eastward defined as positive. Resolving wavelength for dip slip (e) and strike slip (g) constructed by the algorithm shown in Fig. S2. Data fitting and eigenvalues can be found in Fig. S3. To further understand the spatial resolution of the inversion, we developed a technique to qualitatively represent the resolving power of the slip model (Fig. S2). For each kernel used in the construction of the model, we pick out the peaks and troughs, together with the corners and sides, to create a set of nodes. Next, we apply Delaunay triangulation (Delaunay 1934) to these nodes and set their values to the average of the distances between each node and its adjacent nodes. Then we fit these distance measures at each vertex using splines in tension (Smith & Wessel 1990) and filter the result to get the resolving wavelength map for the kernel. When applying this to all the kernels used in the spectral expansion inversion, the minimum value among all the kernels at every position on the fault plane is considered as the resolving wavelength for the model (Figs 2e and g). We also performed further tests for dip slip and a mixture of dip and strike slip scenarios, the results of which are presented in supplementary information (Figs S7 and S8). The two test cases reveal some interesting features of this spectral expansion approach. First, the eigenvalues derived from the Gram matrix decrease slower in the 3-D case than the 2-D case. We believe this is because the data in the 3-D case requires extra constraints for one more dimension, thus more eigenvalues are needed to achieve the same level of trade-offs. Second, the resolving power of the model depends a lot on data sampling. It's obvious, even from the equations in Section 2, that the shapes of kernels are completely determined by the model geometry and distribution of the data samples. The derived resolving wavelength could also act as a guide for the digitization of parameter estimation models. Third, the 2-D case shows an increase in uncertainty at shallower depth since it has complete surface coverage. The 3-D case shows a similar feature, but suddenly the uncertainty becomes smaller at the surface, which we relate to the lack of near-fault data. The truth is that the data can never be dense enough when it comes close to the fault, as the Green's function will change dramatically over a short distance. This inversion approach relies heavily on the similarity of Green's function between data points, that is the redundancy of data. Thus, in this case, the quad-tree algorithm may not be optimal for down sampling, as the redundancy of data is dependent on fault depth. Fourth, although the computation is defined as the inner product of the Green’ functions over the model space, that is how similar their Green's functions are, the Gram matrix could potentially represent the covariance (after scaling by the data's uncertainty) between data points, if we assume the data having zero mean. In other words, if the measurement of one data point changes away from zero due to an earthquake, the measurement of a nearby data point may also change correspondingly. Finally, the fault could be in any shape as long as the integration shown in eq. (3) can be done properly, which gives a strong adaptability to this approach. 4 APPLICATION TO THE 2010 MAULE AND THE 2015 GORKHA EARTHQUAKES We first use this spectral expansion approach to model the slip for the Mw 8.8 Maule, Chile earthquake. This earthquake was a megathrust event that ruptured a mature seismic gap (Moreno et al.2010) and generated a strong tsunami throughout the Pacific Ocean (Fritz et al.2011). Previous study shows coseismic rupture extending downdip to a depth of ∼40 km, which correlates well with interseismic coupling along the Andean subduction zone in south-central Chile (Lorito et al.2011; Métois et al. 2012). Here we apply the spectral expansion approach to the processed InSAR and GPS data presented in (Tong et al.2010) and analyse the uncertainty and resolving power of the model. We use the same fault geometry as in their study, with the fault plane extending 680 km along strike N16.5°E and 250 km downdip 15° towards east. From the misfit function analysis (Fig. 3e), we determined that 150 kernels were needed to represent the coseismic slip. The solution (Figs 3a and b) shows a total slip moment corresponding to moment magnitude 8.79 assuming an average shear modulus of 40 GPa, with the thrust component being 1.86  × 1022Nm and the strike component being 5.85  × 1021 Nm. This is close to the seismic moment reported by (NEIC 2010) (Mw 8.8). From the resolving wavelength analysis (Figs 3c and d), the model has the best resolution at intermediate depth (20–40 km), with resolution decreasing dramatically outside of this range. More interestingly, at shallower depth, due to the lack of offshore data, the resolution is actually worse than at deeper depth. This brings into question the very small uncertainty estimates near the fault surface. Are we so sure the uncertainty of the slip estimate is almost zero close to fault trace? The answer is yes, with a sacrifice in resolution, the statistical precision is almost intrinsically guaranteed. The cumulative dip slip (Fig. 5a) increases from 0 to 20 km depth and decreases to zero at roughly 45 km depth. This is in good agreement with previous studies and the uncertainty of our model together with the resolving power analysis show that the reduction in slip amplitudes is well determined by observations. Figure 3. View large Download slide Inverted slip for the 2010 Mw 8.8 Maule Earthquake. (a) Inverted dip slip and (b) inverted strike slip, with same positivity definition as Figs 2(d) and (f). The black dashed line in (a) represents the fault trace. Resolving wavelength for dip slip (c) and strike slip (d). (e) Misfit versus the number of kernels used. The data fitting, uncertainty estimate and eigenvalues can be found in Fig. S4. Figure 3. View large Download slide Inverted slip for the 2010 Mw 8.8 Maule Earthquake. (a) Inverted dip slip and (b) inverted strike slip, with same positivity definition as Figs 2(d) and (f). The black dashed line in (a) represents the fault trace. Resolving wavelength for dip slip (c) and strike slip (d). (e) Misfit versus the number of kernels used. The data fitting, uncertainty estimate and eigenvalues can be found in Fig. S4. Figure 5. View largeDownload slide Cumulative dip slip versus depth plots for (a) the 2010 Mw 8.8 Maule Earthquake and (b) is the 2015 Mw 7.8 Gorkha Earthquake, with blue lines denoting the cumulative dip slip and magenta dash-dotted lines representing the uncertainty bounds. Figure 5. View largeDownload slide Cumulative dip slip versus depth plots for (a) the 2010 Mw 8.8 Maule Earthquake and (b) is the 2015 Mw 7.8 Gorkha Earthquake, with blue lines denoting the cumulative dip slip and magenta dash-dotted lines representing the uncertainty bounds. We then reinvestigate the 2015 Mw 7.8 Gorkha earthquake in Nepal. This earthquake took place on a tectonic boundary due to the Indian-Eurasian collision and ruptured 140 km on the Main Himalayan Fault (Galetzka et al. 2015). In contrast to megathrust events at subduction zones, previous studies all show a very limited downdip rupture extent with co-seismic slip typically terminating between 15 and 25 km depth (Wang & Fialko 2015). Here we apply the spectral expansion inversion to the InSAR data from (Galetzka et al. 2015; Lindsey et al.2015) and GPS data from (Galetzka et al. 2015). We down-sample the InSAR data using the quad-tree algorithm with increased sampling density closer to the fault in order to increase the resolution at shallow depth. The sampling stops 20 km away from fault in order to maintain numerical stability and save computation time. We then combine the InSAR and GPS data to perform the spectral expansion inversion. The model is created following the main frontal thrust geometry from (Ader et al.2012), with the fault plane extending 200 km along strike N71.5°W and 150 km along dip 7° towards north (Wang & Fialko 2015). By examining the misfit function, we determined that 250 kernels were needed to represent the coseismic slip (Fig. 4e). The estimated total moment for the coseismic slip (Figs 4a and b) is 6.36  × 1021 Nm, with the corresponding moment magnitude being 7.80, assuming an average shear modulus of 30 GPa, which agrees with the seismic moment reported by (NEIC 2015) (Mw 7.8). From the resolving wavelength analysis (Figs 4c and d), the resolving power generally decays with depth except in areas with denser samples due to the quad-tree sampling. The cumulative dip slip with depth plot (Fig. 5b) shows slip rapidly increasing at about 7 km depth, reaches its maximum around 10 km, and then decreases towards zero gradually. The slip versus depth profile is in agreement with the previous studies (Wang & Fialko 2015). From the uncertainty estimates and the resolving power analysis, this decay is also well resolved by the data. Figure 4. View large Download slide Inverted slip for the 2015 Mw 7.8 Gorkha Earthquake. (a) Inverted dip slip with same positivity definition as Fig. 2(d) and (b) is the inverted strike slip with positive defined as the hanging wall moving north-westward. The black dashed line in (a) represents the fault trace. Resolving wavelength for dip slip (c) and strike slip (d). (e) Misfit function versus the number of kernels used. The data fitting, uncertainty estimate and eigenvalues can be found in Fig. S5. Figure 4. View large Download slide Inverted slip for the 2015 Mw 7.8 Gorkha Earthquake. (a) Inverted dip slip with same positivity definition as Fig. 2(d) and (b) is the inverted strike slip with positive defined as the hanging wall moving north-westward. The black dashed line in (a) represents the fault trace. Resolving wavelength for dip slip (c) and strike slip (d). (e) Misfit function versus the number of kernels used. The data fitting, uncertainty estimate and eigenvalues can be found in Fig. S5. For these two cases, the InSAR data's uncertainty is taken from previous studies, with 10 cm for all InSAR data in the Maule case (Tong et al.2010) and 2.3, 5.4 and 4.1 cm for track 47, 48 and 157 in the Nepal case (Wang & Fialko 2015). One interesting feature from the spectral expansion is that the eigenvalues computed for the Maule case (Fig. S4i) decreases faster than for the Nepal case (Fig. S5i), potentially because the Nepal earthquake has a shallower dipping angle or more data close to trench. 5 DISCUSSIONS 5.1 The influence of erosion rate on thermal structure The inversion approach we developed in this study confirms that the contrasting downdip rupture extents of the two earthquakes are indeed well determined by observations. This raises the next question: Why is the downdip rupture extent of Himalayan megathrust earthquakes (15–20 km) a factor of 2 shallower than the downdip limit (>40 km) of megathrust earthquakes in oceanic subduction zones? Temperature exerts a strong influence on both dehydration reactions and deformation mechanisms. In particular, the transition from friction to intracrystalline plasticity as the dominant deformation mechanism in quartzofeldspathic rocks occurs at temperatures of 325–350 °C, and is often proposed to define the upper temperature bound for seismic behaviour in crustal rocks (Hyndman & Wang 1993, 1995; Hyndman et al.1997; Oleskevich et al. 1999). It has been also proposed that earthquake ruptures initiated below this temperature bound may propagate with decreasing slip to where the temperature is 450 °C (Hyndman et al. 1997; Hyndman 2007). Due to the existence of subduction, the temperature is generally higher on the side of the overriding plate, which is roughly around 600 °C (Klingelhoefer et al.2010). The rapid subduction of cooler oceanic lithosphere will depress the isotherms on the megathrust interface, but this can account for only a 5 km deepening when subduction rate changes from 2 to 4 cm yr–1 (Klingelhoefer et al.2010). To further understand what is causing the factor of 2 difference between the Himalayian rupture depth and the Maule rupture depth, we follow Royden (1993) and developed a 1-D heat conduction model that accounts for the effect of erosion rate (Appendix). We assume that the surface is at a quasi-steady stage where the uplift from accretion and rebound is balanced by the effect of erosion. We then adopt the parameters in Royden (1993) who assumes the erosion rate to be 1 mm yr–1 and vary the remaining variables to calculate the corresponding temperature profiles (Fig. 6). From our analysis, both the heat production in the upper crust and the mantle temperature has little effect on depth of the 600 °C isotherm, while the erosion rate influences the temperature profile significantly. This can be understood by considering the effects of accretion or erosion/rebound; the hotter materials are directly brought up towards the surface, which is far more efficient than heat conduction. Recent studies show that in the region of the Gorkha earthquake, the surface erosion rate increases from about 0.5 mm yr–1 in the up-dip region to 3.5 mm yr–1 downdip (Whipple et al.2016), while for south-central Chile, this number is quite small (<0.3 mm yr–1) (Aguilar et al.2011; Carretier et al.2013). Fig. 6(a) shows that an average 1 mm yr–1 erosion rate is sufficient to bring the 600 °C isotherm from a depth of 40–17 km, which is consistent with our findings from spectral expansion inversion. Figure 6. View largeDownload slide Temperature versus depth for 1-D heat conduction model. (a) Models with different accretion/erosion rate. (b) Models with different upper crust heat generation rate. (c) Models with different mantle temperature. Figure 6. View largeDownload slide Temperature versus depth for 1-D heat conduction model. (a) Models with different accretion/erosion rate. (b) Models with different upper crust heat generation rate. (c) Models with different mantle temperature. 5.2 Correlations between forearc ridges and the downdip limit of the seismogenic zone Bassett & Watts (2015a,b) have developed a spectral averaging method of suppressing the large-amplitude, long-wavelength, trench-normal topographic and gravimetric expression of subduction zones. The global application of this technique revealed that in at-least five circum-Pacific subduction zones with a history of large megathrust earthquakes, the downdip limit of the seismogenic zone was approximately correlated with the location of a trench-parallel forearc ridge (TPFR). The longest TPFR parallels the coastline of Chile for >2000 km and is coincident with the Chilean coastal Cordillera. As shown in Fig. 7, the TPFR in southern Chile is well correlated with the maximum depth of earthquakes with Global Centroid Moment Tensors (GCMT; Dziewonski et al. 1981; Ekström et al.2012) consistent with megathrust slip, the maximum depth of co-seismic slip in the 2010 Mw 8.8 Maule earthquake (Tong et al.2010; Delouis et al.2010; this study), and the maximum depth of interseismic coupling as constrained by campaign GPS measurements (Moreno et al.2010). A similar correlation has been recognized in northern Chile (Béjar-Pizarro et al.2013). Figure 7. View largeDownload slide (a) Residual free-air gravity anomaly at the 2010 Chile Mw 8.8 earthquake region with grey dashed line denotes the trench axis and the black dashed line represents the TPFR. Grey pluses mark GCMT thrust events with Mw < 6.5. GCMT events with Mw ≥ 6.5 are plotted with beach ball, scaled for magnitude, coloured for depth. (b) Same as (a) with grey solid contours are the inverted slip pattern with each line representing 3 m slip. (c) Fraction of Plate Convergence from (Moreno et al.2010). The grey dashed line and black dashed line are the same as in (a) and (b). Figure 7. View largeDownload slide (a) Residual free-air gravity anomaly at the 2010 Chile Mw 8.8 earthquake region with grey dashed line denotes the trench axis and the black dashed line represents the TPFR. Grey pluses mark GCMT thrust events with Mw < 6.5. GCMT events with Mw ≥ 6.5 are plotted with beach ball, scaled for magnitude, coloured for depth. (b) Same as (a) with grey solid contours are the inverted slip pattern with each line representing 3 m slip. (c) Fraction of Plate Convergence from (Moreno et al.2010). The grey dashed line and black dashed line are the same as in (a) and (b). To investigate whether a similar relation is observed for the Gorka earthquake, we have applied this technique of spectral averaging to calculate residual topographic anomalies along the Himalayas. We have analysed the global SRTM15 topographic grid, which has a resolution of 15 arcsec (roughly 500 m). Residual topographic anomalies along the Himalayas are shown in Fig. 8. The Gorka earthquake occurred immediately up-dip of an elongated, strike-parallel region of positive residual topographic anomalies, which correlates remarkably well with the 3500 m elevation contour of the Himalayas (Avouac 2003; Stevens & Avouac 2015). This elevated band may be analogous to the TPFR previously observed in oceanic subduction zones. Black contours show the distribution of interseismic coupling as constrained by campaign GPS measurements (Stevens & Avouac 2015). It is intriguing that we observe a similar spatial association between the downdip transition from an interseismically locked to creeping megathrust, with the elevated residual topography of the high Himalayas. Figure 8. View largeDownload slide Residual topography for the 2015 Mw 7.8 Nepal earthquake region. The grey dashed line denotes the Main Frontal Thrust (MFT). The grey lines are the locking contours from (Stevens & Avouac 2015). The grey lines in the subplot at upper-right corner are 1-m slip contours. The bottom-left plot is the average elevation at this region. Figure 8. View largeDownload slide Residual topography for the 2015 Mw 7.8 Nepal earthquake region. The grey dashed line denotes the Main Frontal Thrust (MFT). The grey lines are the locking contours from (Stevens & Avouac 2015). The grey lines in the subplot at upper-right corner are 1-m slip contours. The bottom-left plot is the average elevation at this region. In subduction zones, the broad pattern of interseismic deformation can be modelled by placing an edge dislocation at the downdip end of the locked fault with its Burger's vector parallel to fault dip (Savage 1983). This model predicts a broad interseismic subsidence over the majority of the locked fault, with interseismic uplift located near the downdip limit of locking. Over many earthquake cycles, the incomplete recovery of interseismic elastic deformation by earthquakes and post-seismic slip may lead to permanent deformation of the forearc that spatially mirrors the sign of interseismic deformation. This model has been proposed to explain long-term coastal uplift in oceanic subduction zones in Cascadia (Kelsey et al. 1994) and Mexico (Ramirez-Herrera et al.2004) and we suggest unrecovered interseismic deformation may also be contributing to the elevation and growth of the high Himalaya in eastern Nepal, which is in agreement with previous studies (Meade 2010). 6 CONCLUSIONS In this study, we developed and tested a data-driven spectral expansion approach that possess the advantages of (1) mathematical robustness, (2) smoothness and stableness without regularization, (3) capability of adopting complex fault surface, (4) giving real uncertainty estimates and (5) providing resolving power analysis. With this approach, we confirmed that the downdip ruptures of the 2010 Mw 8.8 Maule earthquake and the 2015 Mw 7.8 Gorkha earthquake are well resolved by the data. The Maule earthquake rupture terminates around 45 km depth and the Gorkha earthquake rupture stops shallower than 20 km. Both these depths are well correlated with the location of trench-parallel forearc ridges, which may reflect unrecovered interseismic elastic deformation. To understand the difference in maximum rupture depth between the two cases, we constructed 1-D thermal conduction models with different erosion rate, upper-crust heat production and mantle temperature. By comparing different models, we conclude that the erosion rate has a major influence on the downdip rupture limit. ACKNOWLEDGEMENTS We want to thank the two reviewers for their valuable suggestions. This study was funded by the NASA Earth Surface and Interior Program (NNX16AK93G) and the Southern California Earthquake Center (SCEC). SCEC is funded by the NSF cooperative Agreement EAR-1033462 and USGS Cooperative Agreement G12AC20038. Special thanks to Prof Catherine Constable for sharing with us her knowledge on Geophysical Inverse Theory. We also want to thank the authors from Tong et al. (2010), Lindsey et al. (2015) and Galetzka et al. (2015) for sharing their data to us. REFERENCES Ader T. et al.  , 2012. Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: implications for seismic hazard, J. geophys. Res. , 117, B04403, doi:10.1029/2011JB009071. Google Scholar CrossRef Search ADS   Aguilar G. et al.  , 2011. Variability in erosion rates related to the state of landscape transience in the semi‐arid Chilean Andes, Earth Surface Processes and Landforms , 36( 13), 1736– 1748. Ammon C.J. et al.  , 2005. Rupture process of the 2004 Sumatra-Andaman earthquake, Science , 308( 5725), 1133– 1139. Google Scholar CrossRef Search ADS PubMed  Avouac J.P., 2003. Mountain building, erosion, and the seismic cycle in the Nepal Himalaya, Adv. Geophys. , 46, 1– 80. Google Scholar CrossRef Search ADS   Bassett D., Watts A.B., 2015a. Gravity anomalies, crustal structure, and seismicity at subduction zones: 1. Seafloor roughness and subducting relief, Geochem. Geophys. Geosyst. , 16( 5), 1508– 1540. Google Scholar CrossRef Search ADS   Bassett D., Watts A.B., 2015b. Gravity anomalies, crustal structure, and seismicity at subduction zones: 2. Interrelationships between fore‐arc structure and seismogenic behavior, Geochem. Geophys. Geosyst. , 16( 5), 1541– 1576. Google Scholar CrossRef Search ADS   Béjar-Pizarro M., Socquet A., Armijo R., Carrizo D., Genrich J., Simons M., 2013. Andean structural control on interseismic coupling in the North Chile subduction zone, Nat. Geosci. , 6( 6), 462– 467. Google Scholar CrossRef Search ADS   Carretier S. et al.  , 2013. Slope and climate variability control of erosion in the Andes of central Chile, Geology , 41( 2), 195– 198. Google Scholar CrossRef Search ADS   Chlieh M. et al.  , 2007. Coseismic slip and afterslip of the great Mw 9.15 Sumatra–Andaman earthquake of 2004, Bull. seism. Soc. Am. , 97( 1A), S152– S173. Google Scholar CrossRef Search ADS   Cohen S.C., 1999. Numerical models of crustal deformation in seismic zones, Adv. Geophys. , 41, 133– 231. Google Scholar CrossRef Search ADS   Delaunay B., 1934. Sur la sphere vide. Izv. Akad. Nauk SSSR, Otdelenie Matematicheskii i Estestvennyka Nauk , 7( 793-800), 1– 2. Delouis B., Nocquet J.-M., Vallée M., 2010. Slip distribution of the February 27, 2010 Mw = 8.8 Maule earthquake, central Chile, from static and high‐rate GPS, InSAR, and broadband teleseismic data, Geophys. Res. Lett. , 37, L17305, doi:10.1029/2010GL043899. Google Scholar CrossRef Search ADS   Dziewonski A.M., Chou T.A., Woodhouse J.H., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, J. geophys. Res. , 86( B4), 2825– 2852. Google Scholar CrossRef Search ADS   Ekström G., Nettles M., Dziewoński A.M., 2012. The global CMT project 2004–2010: Centroid-moment tensors for 13,017 earthquakes, Phys. Earth planet. Inter. , 200, 1– 9. Google Scholar CrossRef Search ADS   Fialko Y., 2004. Probing the mechanical properties of seismically active crust with space geodesy: study of the coseismic deformation due to the 1992 Mw7.3 Landers (southern California) earthquake, J. geophys. Res. , 109, B03307, doi:10.1029/2003JB002756. Fritz H.M. et al.  , 2011. Field survey of the 27 February 2010 Chile tsunami, Pure appl. Geophys. , 168( 11), 1989– 2010. Google Scholar CrossRef Search ADS   Galetzka J. et al.  , 2015. Slip pulse and resonance of the Kathmandu basin during the 2015 Gorkha earthquake, Nepal, Science , 349( 6252), 1091– 1095. Google Scholar CrossRef Search ADS PubMed  Gilbert F., 1971. Ranking and winnowing gross Earth data for inversion and resolution, Geophys. J. Int. , 23( 1), 125– 128. Google Scholar CrossRef Search ADS   Hyndman R.D., 2007. The Seismogenic Zone of Subduction Thrust Faults , Columbia Univ. Press, pp. 15– 40. Hyndman R.D., Wang K., 1993. Thermal constraints on the zone of major thrust earthquake failure: the Cascadia subduction zone, J. geophys. Res. , 98( B2), 2039– 2060. Google Scholar CrossRef Search ADS   Hyndman R.D., Yamano M., Oleskevich D.A., 1997. The seismogenic zone of subduction thrust faults, Island Arc , 6( 3), 244– 260. Google Scholar CrossRef Search ADS   Hyndman R.D., Wang K., Yamano M., 1995. Thermal constraints on the seismogenic portion of the southwestern Japan subduction thrust, J. geophys. Res. , 100( B8), 15 373–15 392. Google Scholar CrossRef Search ADS   Jónsson S., Zebker H., Segall P., Amelung F., 2002. Fault slip distribution of the 1999 Mw 7.1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements, Bull. seism. Soc. Am. , 92( 4), 1377– 1389. Google Scholar CrossRef Search ADS   Kelsey H.M., Engebretson D.C., Mitchell C.E., Ticknor R.L., 1994. Topographic form of the Coast Ranges of the Cascadia margin in relation to coastal uplift rates and plate subduction, J. geophys. Res. , 99( B6), 12 245–12 255. Google Scholar CrossRef Search ADS   Klingelhoefer F. et al.  , 2010. Limits of the seismogenic zone in the epicentral region of the 26 December 2004 great Sumatra‐Andaman earthquake: results from seismic refraction and wide‐angle reflection surveys and thermal modeling, J. geophys. Res. , 115, B01304, doi:10.1029/2009JB006569. Google Scholar CrossRef Search ADS   Lay T. et al.  , 2005. The great Sumatra-Andaman earthquake of 26 December 2004, Science , 308( 5725), 1127– 1133. Google Scholar CrossRef Search ADS PubMed  Lindsey E.O., Natsuaki R., Xu X., Shimada M., Hashimoto M., Melgar D., Sandwell D.T., 2015. Line‐of‐sight displacement from ALOS‐2 interferometry: Mw 7.8 Gorkha Earthquake and Mw 7.3 aftershock, Geophys. Res. Lett. , 42( 16), 6655– 6661. Google Scholar CrossRef Search ADS   Lorito S. et al., 2011. Limited overlap between the seismic gap and coseismic slip of the great 2010 Chile earthquake, Nat. Geosci. , 4( 3), 173– 177. Google Scholar CrossRef Search ADS   Meade B.J., 2010. The signature of an unbalanced earthquake cycle in Himalayan topography?, Geology , 38( 11), 987– 990. Google Scholar CrossRef Search ADS   Métois M., Socquet A., Vigny C., 2012. Interseismic coupling, segmentation and mechanical behavior of the central Chile subduction zone, J. geophys. Res. , 117, B03406, doi:10.1029/2011JB008736. Google Scholar CrossRef Search ADS   Minson S.E. et al.  , 2014. Bayesian inversion for finite fault earthquake source models—II: The 2011 great Tohoku-oki, Japan earthquake, Geophys. J. Int. , 198( 2), 922– 940. Google Scholar CrossRef Search ADS   Moreno M., Rosenau M., Oncken O., 2010. 2010 Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone, Nature , 467( 7312), 198– 202. Google Scholar CrossRef Search ADS PubMed  National Earthquake Information Center (NEIC), 2010. Magnitude 8.8-Offshore Maule, Chile, U.S. Geol. Surv., Denver, CO. Available at: http://earthquake.usgs.gov/earthquakes/eventpage/official20100227063411530_30#executive. National Earthquake Information Center (NEIC), 2015. Magnitude 8.8-Offshore Maule, Chile, U.S. Geol. Surv., Denver, CO. Available at: http://earthquake.usgs.gov/earthquakes/eventpage/official20100227063411530_30#executive. Okada Y., 1985. Surface deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am. , 75( 4), 1135– 1154. Oleskevich D.A., Hyndman R.D., Wang K., 1999. The updip and downdip limits to great subduction earthquakes: Thermal and structural models of Cascadia, south Alaska, SW Japan, and Chile, J. geophys. Res. , 104( B7), 14 965–14 991. Google Scholar CrossRef Search ADS   Parker R.L., 1977. Understanding inverse theory, Ann. Rev. Earth planet. Sci. , 5, 35– 64. Google Scholar CrossRef Search ADS   Parker R.L., 1994. Geophysical Inverse Theory , Princeton Univ. Press. Ramırez-Herrera M.T., Kostoglodov V., Urrutia-Fucugauchi J., 2004. Holocene-emerged notches and tectonic uplift along the Jalisco coast, Southwest Mexico, Geomorphology , 58( 1), 291– 304. Google Scholar CrossRef Search ADS   Royden L.H., 1993. The steady state thermal structure of eroding orogenic belts and accretionary prisms, J. geophys. Res. , 98, 4487– 4507. Google Scholar CrossRef Search ADS   Savage J.C., 1983. A dislocation model of strain accumulation and release at a subduction zone, J. geophys. Res. , 88( B6), 4984– 4996. Google Scholar CrossRef Search ADS   Simons M., Fialko Y., Rivera L., 2002. Coseismic deformation from the 1999 Mw 7.1 Hector Mine, California, earthquake as inferred from InSAR and GPS observations, Bull. seism. Soc. Am. , 92( 4), 1390– 1402. Google Scholar CrossRef Search ADS   Simons M. et al.  , 2011. The 2011 magnitude 9.0 Tohoku-Oki earthquake: mosaicking the megathrust from seconds to centuries, Science , 332( 6036), 1421– 1425. Google Scholar CrossRef Search ADS PubMed  Smith W.H.F., Wessel P., 1990. Gridding with continuous curvature splines in tension, Geophysics , 55( 3), 293– 305. Google Scholar CrossRef Search ADS   Steketee J.A., 1958. On Volterra's dislocations in a semi-infinite elastic medium, Can. J. Phys. , 36( 2), 192– 205. Google Scholar CrossRef Search ADS   Stevens V.L., Avouac J.P., 2015. Interseismic coupling on the main Himalayan thrust, Geophys. Res. Lett. , 42( 14), 5828– 5837. Google Scholar CrossRef Search ADS   Tarantola A., 2005. Inverse Problem Theory and Methods for Model Parameter Estimation , Chap. 1, pp. 24– 40, SIAM. Tong X. et al., 2010. The 2010 Maule, Chile earthquake: downdip rupture limit revealed by space geodesy, Geophys. Res. Lett. , 37, L24311, doi:10.1029/2010GL045805. Google Scholar CrossRef Search ADS   Wang K., Fialko Y., 2015. Slip model of the 2015 Mw 7.8 Gorkha (Nepal) earthquake from inversions of ALOS‐2 and GPS data, Geophys. Res. Lett. , 42( 18), 7452– 7458. Google Scholar CrossRef Search ADS   Whipple K.X., Shirzaei M., Hodges K.V., Arrowsmith J.R., 2016. Active shortening within the Himalayan orogenic wedge implied by the 2015 Gorkha earthquake, Nat. Geosci. , 9( 9), 711– 716. Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. (a) Comparison between surface deformation derived from the input slip (red dots, noise added) and inverted slip (7 kernels, black line) corresponding to Fig. 2(c). (b) Misfit function versus the number of kernels used for Fig. 2(c). (c) Comparison between surface deformation corresponding to Fig. 2(d) (11 kernels). (d) Misfit function plot versus number of kernels used for Fig. 2(d). Figure S2. (a) An example 2-D kernel for the 2010 Maule Earthquake, with black circles denoting the peaks and troughs of the kernel. (b) Delaunay triangulation for the nodes (circles) in (a) (adding 8 nodes at corners and middle of sides, value of each node is the average distance from the node to its adjacent nodes), with colour representing the average from the triangle's three nodes. (c) Map generated applying a surface technique and filtering to the nodes in (b). Figure S3. Plots (a), (b) and (c) are the subsampled data, forwarded model, and misfit for the simulated noisy ascending data. Plots (d), (e) and (f) are the subsampled data, forwarded model, and misfit for the simulated noisy descending data. (g) Uncertainty map for dip slip in Fig. 3(d). (h) Uncertainty map for strike slip in Fig. 3(f). (i) Eigenvalues versus order. Figure S4. Plots (a), (b) and (c) are the subsampled data, forwarded model, and misfit for the ascending data for the 2010 Mw 8.8 Maule Earthquake. Plots (d), (e) and (f) are the subsampled data, forwarded model, and misfit for the descending data. (g) Uncertainty map for dip slip in Fig. 4a. (h) Uncertainty map for strike slip in Fig. 4b. (i) Eigenvalues versus order. Figure S5. Plots (a), (b) and (c) are the subsampled data, forwarded model, and misfit for the data from ALOS-2 ascending track 157 for the 2015 Mw 7.8 Nepal earthquake. Plots (d), (e) and (f) are the subsampled data, forwarded model, and misfit for the data from ALOS-2 descending track 47. Plots (g), (h) and (i) are the subsampled data, forwarded model, and misfit for the data from ALOS-2 descending track 48. (j) Uncertainty map for dip slip in Fig. 5a. (k) is the uncertainty map for strike slip in Fig. 5b. (l) Eigenvalues versus order. Figure S6. Some examples of kernels for the 3-D test. The upper row shows the dip slip kernels and the lower row shows the strike slip kernels. The wavelength of the kernel decreases as the order increases. Figure S7. 3-D strike-slip test case. (a) Input strikeslip with the black dashed line denoting the fault trace. (b) Cumulative strike slip $$p( z ) = \mathop \int \nolimits_0^L m( {x,z} ){\rm{d}}x$$ versus depthwith black dashed line being the input strike slip, blue line representing the inverted strike slip, and magenta dash-dot lines denoting the uncertainty. (c) Misfit versus order with the red dot being the order selected. (d) Inverted dip slip and (f) inverted strike slip. Resolving wavelength for dip slip (e) and strike slip (g) constructed by the algorithm shown in Fig. S2. Figure S8. 3-D mixture of dip and strike sliptest case. (a) Input strike and dip slip (same amplitude) with the black dashed line denoting the fault trace. (b) Cumulative slip $$p( z ) = \mathop \int \nolimits_0^L m( {x,z} ){\rm{d}}x$$ versus depth plot with black dashed line being the input dip/strike slip, blue line representing the inverted dip slip and cyan dash-dot lines denoting the uncertainty, red line representing the inverted strike slip, and magenta dash-dot lines denoting the uncertainty. (c) Misfit versus order with the red dot being the order selected. (d) Inverted dip slip and (f) inverted strike slip. Resolving wavelength for dip slip (e) and strike slip (g) constructed by the algorithm shown in Fig. S2. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX: APPENDIX: DERIVATION OF THE 1-D HEAT CONDUCTION MODEL WITH ACCRETION/EROSION Following (Royden 1993), we simplified and modified the heater conduction equation as follows, assuming the accretion and erosion rates are in balance.   \begin{equation*} \frac{{{\rm d^2}T}}{{{\rm{d}}{z^2}}} + \frac{v}{\alpha }\frac{{{\rm{d}}T}}{{{\rm{d}}z}} + \frac{{{A_0}}}{K}\ {{\rm{e}}^{ - z/h}} = \ 0, \end{equation*} where T is the temperature, z is the depth, v is the accretion rate, α is the thermal diffusivity, A0 is the heat production rate at the surface, K is the thermal conductivity and h is the characteristic decaying depth for heat production. We simplify the equation a little more by substituting some of the variables as follows:   \begin{equation*} \frac{{{{\rm{d}}^2}T}}{{{\rm{d}}{z^2}}} + a\frac{{{\rm{d}}T}}{{{\rm{d}}z}} + b\ {{\rm{e}}^{ - cz}} = \ 0 \end{equation*} Then further simplify the equation form a 2nd order ODE to a 1st order ODE by letting dT/dz = y. Then we have the equation as   \begin{equation*} \frac{{{\rm{d}}y}}{{{\rm{d}}z}} + ay + b\ {{\rm{e}}^{ - cz}} = \ 0 \end{equation*} As for the type of 1st order ODE like y' + P(x)y = Q(x), the general solution is   \begin{equation*} y\ = {e^{ - \mathop \int \nolimits_{{x_0}}^x P\left( t \right){\rm{d}}t}}\ \int Q\left( x \right){{\rm{e}}^{\mathop \int \nolimits_{{x_0}}^x P\left( t \right){\rm{d}}t}}{\rm{d}}x. \end{equation*} Thus, we have the solution for y as   \begin{equation*} y\ = \ - \frac{b}{{a - c}}\left( {{e^{ - cz}} + {C_1}{e^{ - az}}} \right) \end{equation*} Substitute y with dT/dz and solve the ODE, we’ll get   \begin{equation*} T\ = \frac{b}{{c\left( {a - c} \right)}}\ {e^{ - cz}} + {C_1}\frac{b}{{a\left( {a - c} \right)}}{e^{ - az}} + {C_2} \end{equation*} Consider the boundary condition   \begin{equation*} z\ = \ 0,\quad \ T\ = {T_0}{\rm{\ }} \end{equation*}   \begin{equation*} z\ = \ {\rm{H}},\quad \ T\ = {T_m}\ \end{equation*} We can get the coefficients   \begin{equation*} \ {C_1} = \ - \left( {{T_m}\frac{{a\left( {a - c} \right)}}{b} + \frac{a}{c}\left( {1 - {e^{ - cH}}} \right)} \right)\bigg/\left( {1 - {e^{ - aH}}} \right) \end{equation*}   \begin{equation*} \ {C_2} = \left( {{T_m} + \frac{b}{{c\left( {a - c} \right)}}\left( {{e^{ - aH}} - {e^{ - cH}}} \right)} \right)\ \bigg/\left( {1 - {e^{ - aH}}} \right) \end{equation*} © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations