SUMMARY
The slip and static stress drop of the 2015 Mw 7.8 Gorkha, Nepal earthquake has been studied using its excellent near-source static observations and a newly developed finite fault inversion algorithm with an average stress drop constraint. A series of nonlinear inversions with different target stress drops are conducted to search for the solution that not only most accurately fits the geodetic data, but also has an energy-based stress drop, |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$|, matching the target value. Our result reveals that the misfit to the geodetic data gradually decreases when the |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| of the inverted slip model increases from 2 to 7 MPa; but becomes nearly constant when |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| further increases. Hence, only the lower bound of |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$|, that is, |${\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}} ^{\mathrm{ min}}}$|(∼7 MPa), of the Gorkha earthquake can be constrained with near-fault geodetic measurements, consistent with our previous study using just far-field seismic data. The upper bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| can be constrained with an extra constraint on the maximum local stress drop. An artificial upper bound can also be reached if using a large subfault size to represent the source. The lower bound of |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| leads to the lower bound of the apparent available energy and the upper bound of the radiation efficiency (|${\eta _\mathrm{ R}}$|, 0.09–0.15), though the latter is also sensitive to the determination of the radiated seismic energy. We find that the lower |${\eta _\mathrm{ R}}$| and the reported high rupture velocity (80–93 per cent shear wave speed) can be reconciled by considering the aspect ratio of the dominant slip patch. We recommend using |${\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}} ^{\mathrm{ min}}}$| to replace various slip smoothing constraints to stabilize the finite fault inversion because of its clear physical meaning.
1 INTRODUCTION
Second, in order to explore the evolution of slip on faults, the source process is parametrized to be as realistic as possible (e.g. Hartzell & Helmberger 1982; Olson & Apsel 1982; Ide et al.1996; Ji et al.2002b). Hundreds of free parameters are often used, despite the fact that not all of them are well constrained by the available observations. Consequently, regularizations, such as smoothing constraints applied to the fault slip, are implemented to stabilize the solutions. Thus, the uncertainties of individual parameters are often difficult to access because they are not only caused by the limited data, but also these regularizations. Different constraints inevitably result in differences in inverted slip model distributions (e.g. Hartzell et al.1991), though the long-wavelength features tend to be well preserved.
As mentioned above, various regulations have to be applied to finite fault studies to stabilize the inversions. These regulations inevitably affect the inverted solutions but how to appropriately select them is a challenging task for source modellers. For instance, slip-smoothing constraints are used in most inversion algorithms (e.g. Hartzell & Heaton 1983). As pointed out by Yabuki & Matsu'ura (1992), the stress field in and around the fault region will be finite. The finiteness of stress thus requires the smooth variation of slip along the fault surface. Subsequently, the true target of these finite fault inversion algorithms becomes ‘the most spatially smooth slip model that still explains the observations well’. The definition of the slip model roughness, which is concurrently minimized during the inversions, is, however, not unique. Hartzell & Heaton (1983) favoured a simple smoothing scheme, that is, |${D_i} - {D_j} \cong 0$|, where |${D_i}$| and |${D_j}$| are the slip on adjacent subfaults. Du et al. (1992) defined the fault roughness as the Laplacian of the inverted slip distribution along the fault surface. If the finite fault inversion is carried out with the singular value decomposition method, the smoothing of the fault slip might also be achieved by properly truncating the small eigenvalues (e.g. Hartzell & Heaton 1983). These methods impact the inverted slip distribution differently (Adams et al.2017), and there is no physical reason to select one approach over the others.
Adams et al. (2017) noted the good correlation between the Laplacian of the slip distribution and energy-based average stress drop |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| and suggested that the |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| can be used as an alternative way to smooth the fault slip. Unlike the constraints mentioned above, this regulation, which we name as ‘minimum average stress drop’ constraint, has a clear physical meaning. Because of the aforementioned relationship in eq. (2), and the fact that the seismic potency of a large earthquake can often be well constrained, searching for the lower bound of |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| is equivalent to searching for the lower bound of the available energy (|${\rm{\Delta }}{W_0}$|) of this earthquake. In other words, we can replace the previous target of finite fault inversions, the smoothest slip distribution that can explain the observations well, with a more physical and well-defined new target: the slip model with the lowest available energy that can explain the observations well.
The ability to estimate the lower bound of the available energy is important in studying the physics of the source process. For example, a lower bound for the available energy infers an upper bound for the seismic radiation efficiency (|${\eta _\mathrm{ R}})$|, which is defined as the ratio of the radiated seismic energy to the available energy, |${\eta _\mathrm{ R}} = \frac{{{E_\mathrm{ R}}}}{{{\rm{\Delta }}{W_0}}}$| (e.g. Husseini 1977; Venkataraman & Kanamori 2004). Husseini (1977) and Venkataraman & Kanamori (2004) showed that |${\eta _\mathrm{ R}}$| is positively correlated with the average rupture velocity. Thus, a rupture with a small radiation efficieny implies more energy is dissipated as fracture energy to create the rupture plane. Ye et al. (2016) systematically analysed 110 Mw > 7 earthquakes in the world and reported the average value of |${\eta _\mathrm{ R}}$| to be 0.34 (logarithmic) or 0.49 (linear).
Here we study the slip distribution and static stress drop of the 2015 Mw7.8 Gorkha, Nepal earthquake, which has some of the best near-field geodetic observations for earthquakes with similar or larger magnitudes (e.g. Lindsey et al.2015, Wang & Fialko 2015; Mcnamara et al.2016; Yue et al.2016). The Gorkha earthquake has previously been extensively studied with seismic and geodetic data individually or combined together (e.g. Avouac et al.2015; Denolle et al.2015; Lindsey et al.2015; Lay et al.2016; Mcnamara et al.2016; Yue et al.2016). These studies revealed that this earthquake has a fast average rupture velocity, varying from 2.8 to 3.2 km s−1 (Avouac et al.2015; Fan & Shearer 2015; Zhang et al.2016), 83–92 per cent of the shear wave speed (3.5 km s−1) in the depth of coseismic rupture. In this study, we will first explore the lower and upper bounds of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| during the 2015 Nepal earthquake that are consistent with this excellent near-field static data. Second, a classic L-curve method is introduced to quantitatively define the lower bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$|, that is, |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$|. The lower bound of |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| leads to the upper bound of the radiation efficiency, that is, |${\eta _\mathrm{ R}}^{\mathrm{ max}}$|, which is in fact much smaller than the previous reported mean values of global large earthquakes (Ye et al.2016). We subsequently show that the low |${\eta _\mathrm{ R}}$| and the high rupture velocity can be reconciled by considering the shape of the slip patch. In the end, we recommend using |${\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}} ^{\mathrm{ min}}}$| to replace various slip smoothing constraints in order to stabilize the finite fault inversion.
2 TECTONIC BACKGROUND AND OBSERVATIONS
The 2015 April 25 MW 7.8 Gorkha earthquake occurred at 06:11:25 (UTC) on the Main Himalayan Thrust (MHT) as a result of thrust faulting on or near the main thrust interface between the subducting India plate and the overriding Eurasia plate (Fig. 1; Avouac et al.2015). Previous geodetic studies revealed that the MHT is the primary fault interface accommodating 20 mm yr−1 of convergence between these two plates (Molnar 1988; Larson et al.1999; Avouac 2003). Using interseismic GPS observations, Ader et al. (2012) concluded that the MHT is locked from the surface to approximately 20 km depth. The segment of the MHT where the Gorkha earthquake occurred previously ruptured in 1833, with a rupture length of 100 km (Bilham et al.2001), comparable to the 2015 Gorkha earthquake.
Tectonic setting of the 2015 MW 7.9 Gorkha earthquake with the GCMT solution and focal mechanism shown. The star denotes the USGS epicentre and the blue box shows the surface projection of the fault plane used in the inverse models.
We use coseismic geodetic data distributed to the public domain in this analysis. The unwrapped InSAR data were downloaded from http://topex.ucsd.edu/nepal/, which counts the LOS (line-of-sight) surface deformation between 2015 February 22 and 2015 May 3 measured by Advanced Land Observing Satellite 2 (ALOS-2, operated by Japanese Aerospace Exploration Agency; Lindsey et al.2015). The data have been subsampled to 678 points using the QuadTree resampling method (Lohman & Simons 2005). The coseismic displacements at four near-fault GPS stations deployed by the Caltech Tectonics Observatory (e.g. Avouac et al.2015) were processed by Advanced Rapid Imaging and Analysis Center for Natural Hazards at Jet Propulsion Laboratory (JPL) using GIPSY-OASIS software (Galetzka et al.2015). Fig. 2(a) illustrates the distribution of the LOS displacement measurements, with the locations of the InSAR and GPS points. The colour maps the LOS displacements, which change from −70.9 to 107.6 cm. The southern end shows a larger and positive coseismic LOS offset compared to the northern region and implies a greater uplift compared to the smaller subsidence in the north.
(a) Distribution of the resampled points of InSAR LOS measurements (black triangles) and GPS stations (black circles), accompanied with the image of the unwrapped ALOS-2 LOS displacement field. The colour bar denotes the motion amplitude. The black and red arrows show the GPS observed and synthetic vertical offsets, respectively (see the text for details). (b) Similar to (a) but synthetic image of the LOS displacement field. The black and red arrows show the GPS observed and synthetic horizontal offsets, respectively. (c) Distribution of the resampled points (triangles) with the colours showing the amplitude of the residuals (observed synthetic).
3 METHOD
where Djk and λjk are the dislocation amplitude and rake angle, respectively. |$Y_{jk}^1$| and |$Y_{jk}^2$| are the subfault Green's functions for unit slip in the strike and downdip direction, respectively. Ordinarily, during the finite fault inversion with geodetic data, the fault model m, described by Djk and λjk, is determined by minimizing the difference between the synthetic and observed static displacements.
In this study, we simply define Efit(m) as the root mean square (RMS) of the difference between the observed and synthetic static displacements. Emoment is a function of the discrepancy between the inverted seismic moment and the value based on other methods, such as the GCMT solution. The weight λ1 is fixed at 0.1 in this study for all inversions. However as we discuss later, changing this to a negligible value (10−4) does not affect our conclusions. Finally, the last term is called the ‘fixed average stress drop constraint’, and it is zero when the |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}( m )$| of the inverted slip model m is equal to the pre-assigned |${\tau _{\mathrm{ target}}}$|. The logarithm function is selected here to limit the range of variation due to different |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}( m )$|. Note that as |${\overline {{\rm{\Delta }}\tau } _E}( m )$| approaches |${\tau _{\mathrm{ target}}}$|, |$| {ln( {\frac{{{{\overline {{\rm{\Delta }}\tau } }_E}( \mathrm{ m} )}}{{{\tau _{\mathrm{ target}}}}}} )} |$|converges to zero as |$\frac{{| {{{\overline {{\rm{\Delta }}\tau } }_\mathrm{ E}}( m ) - {\tau _{\mathrm{ target}}}} |}}{{{\tau _{\mathrm{ target}}}}}$|. λτ is the coefficient balancing the effects of Efit(m) and the ‘fixed average stress drop constraint’. It can be selected using a grid search approach (Adams et al.2017).
During this study, we let |${\tau _{\mathrm{ target}}}$| vary within a large range and a simulated annealing algorithm (Ji et al.2002a) is used to search for the minimum of E(m) for each given |${\tau _{\mathrm{ target}}}$| (Adams et al.2017).
4 FAULT GEOMETRY AND VELOCITY STRUCTURE
We approximate the fault with a single rectangular fault plane inferred from the GCMT solution (http://www.globalcmt.org), the USGS epicentre (28.231°N, 84.731°E) and a slightly shallow hypocentre depth of 12 km. It is oriented with a strike of 293° and dips 7° to the north. Similar fault geometry is also used in previous studies (e.g. Yue et al.2016). After some preliminary tests, the final fault plane used in this study extends 176 km along-strike and 80 km downdip, spanning a depth range of 8.3 to 18.1 km. The fault plane is divided into 880 subfaults measuring 4.0 km (along-strike) and 4.0 km (downdip). The fault slip and rake angle of each subevent are inverted, leading to a total of 1760 free parameters. To test the impact of the subfault size, we also conduct inversions using a subfault size of 8 km by 8 km, resulting in a total of 440 free parameters. It is important to recognize that the 4 × 4 km subfault case is underdetermined as we have fewer data points (678) than free parameters, whereas the 8 × 8 km case is not. During the inversion, the slip on individual subfaults is allowed to change from 0.0 to 15 m, and the rake angle is allowed to vary from 80° to 140°. We also limit the slip at the edges of the fault plane to less than 0.01 m.
We use a 1-D layered velocity model (Fig. 3) interpolated from the Global CRUST 2.0 model (Bassin et al.2000) to approximate the earth structure in the source region. It is important to note that from a depth of 0.6 to 21.5 km, which encompasses the depth region of coseismic rupture, the S-wave speed is 3.5 km s, according to this model. The corresponding rigidity is then 3.31 × 104 MPa. While the surface static displacements of the layered structure are computed using the method of Xie & Yao (1989), the constant velocity and density in the source depth range allow us to quickly calculate the on-fault stress change using the analytic solution for a half-space earth (Okada 1992).
1-D Velocity model used in all inversions. The P-wave velocity (km s−1), S-wave velocity (km s−1) and the density (g cm−3) are shown in blue, red and green, respectively.
5 RESULTS
We have conducted three sets of inversions to explore the upper and lower bounds of |${\overline {{\rm{\Delta \tau }}} _{\rm{E}}}$| and the potential uncertainties associated with subfault size and the maximum local stress drop. We will first investigate the results using 4 km by 4 km subfaults and then discuss the uncertainties.
5.1 |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| of the 2015 Gorkha earthquake
We have conducted inversions with 22 different target stress drops (|${\tau _{\mathrm{ target}}})$| ranging from 2 to 50 MPa. The subfault size of all fault models is 4 km by 4 km. Inversions are allowed a variety of different initial models, but the number of iterations is fixed to 400 after preliminary tests. For each |${\tau _{\mathrm{ target}}}$|, we conduct multiple inversions with the coefficient λτ gradually increasing from a negligible value of 10−4 until 0.1 and select the model with the smallest λτ value that still could guarantee the difference between |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| of the inverted model and |${\tau _{\mathrm{ target}}}$| to be less than 10 per cent (Adams et al.2017). Fig. 4 shows the relationship between |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| of the inverted models and the value of the corresponding misfit function |${E_{\mathrm{ fit}}}$|. It can be seen that |${E_{\mathrm{ fit}}}$| generally decreases quickly as |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| increases from 2 to 7 MPa, and then remains essentially constant when |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| increases from 7 to 50 MPa, albeit with small fluctuations. The pattern is identical to our previous work using teleseismic data (Adams et al.2017).
(a) Energy-based average stress drop (|$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$|, MPa) versus RMS misfit (cm) for inversions using a subfault size of 4 × 4 km. Only the fixed average stress drop constraint is used (eq. 6). The misfit first decreases as |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| increases, suggesting it is possible to define a clear lower bound. The average misfit for |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| ≥ 7 MPa is plotted as the black dashed line, and ±2σ of the misfit range due to the choice of the initial model is plotted as the red lines. (b) Comparison of the relationships of the energy-based average stress drop (|$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$|, MPa) and misfit (cm) for inversions using a subfault size of 8 × 8 km (in red) and 4 × 4 km (in blue). Only the fixed average stress drop constraint is used. The misfit in both cases first decreases as |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| increases, suggesting it is possible to define a clear lower bound. However, only the misfit in the tests using 8 × 8 km subfaults increases significantly again as |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| > 25 MPa. (c) Comparison of the relationships of the energy-based average stress drop (|$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| in MPa) and misfit (cm) for inversions using the additional fixed maximum local stress constraint (in blue), and inversions without (in red), both with 4 × 4 km subfaults and the fixed average stress drop constraint. The misfit in both cases first decreases as |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| increases, suggesting that it is possible to define a clear lower bound. The local stress constraint imposes a maximum local stress drop of 25 MPa. Thus, the blue curve cannot converge to any values greater than 25 MPa. The blue curve also increases rapidly as |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| increases and approaches the local maximum.
In principle the simulated annealing algorithm should be independent of the initial models, but in practice, because of the limitations such as a fixed maximum iteration number, small dependences are often observed. We conduct 10 additional inversions with different initial models for the same |${\tau _{\mathrm{ target}}}$| of 7 MPa. The |${E_{\mathrm{ fit}}}$| values of these models have a mean of 2.694 cm and a standard deviation (|$\sigma )$| of 0.011. The two red dashed lines in Fig. 4(a) define the |$\pm 2\sigma$| misfit range due to the choice of the initial model. Note that all |${E_{\mathrm{ fit}}}$| values of models with |${\rm{\Delta }}{\tau _\mathrm{ E}}$|>7 MPa fall within this range. Hence, this excellent geodetic data set allows us to define the lower bound of |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| but fails to constrain its upper bound, similar to our previous conclusion using teleseismic data (Adams et al.2017).
All inverted models with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}} \ge 7\ \mathrm{ M}\mathrm{ Pa}$| explain the data well. Fig. 2(c) shows the distribution of the residuals (observed-synthetics) of the worst scenario, that is, the solution with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}} = 7\ \mathrm{ M}\mathrm{ Pa}$|. The misfits for the majority of the InSAR points across the actual fault plane are near zero. However, relatively large discrepancies are found in the points far to the south and near the bottom left-hand corner of the assigned fault plane. As large misfits for measurements in similar regions can also be seen in other studies, the misfits are likely caused by errors unrelated to the coseismic rupture.
The exact lower bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$|, or |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$|, depends on the errors and limitations of observations, the synthetic earth response and the inversion procedure. If only considering the error caused by the choice of the initial model, we might argue that, with 95 per cent confidence, the |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| of 2015 Nepal earthquake is about 7 MPa (upper red dashed line in Fig. 4a; Adams et al.2017). Alternatively, Adams et al. (2017) found that the |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| of the inverted model changes monotonically with the roughness of the inverted slip model measured as the Laplacian of fault slip. One of the classic methods that is widely used to deal with the trade-off between data misfit and the roughness of the inverted slip model is the L-curve method (e.g. Hansen & O'Leary 1993; Mendoza & Hartzell 2013). We adopt it here to deal with the balance between |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| and |${E_{\mathrm{ fit}}}$|. Analogously, |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| is defined as the ‘corner’ or the point with the maximum curvature of the log(|${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$|)–log(|${E_{\mathrm{ fit}}})$| plot. Following Hansen & O'Leary (1993), a local smoothing polynomial function is applied before a 2-D spline curve is used to fit the set of points that make up the L-curve: which are the average stress drops of the model and their respective inversion misfit values. Then, the point of maximum curvature on this spline curve is computed and the ‘corner’ becomes the data point closest to this calculated point. After applying this algorithm, the corner, that is, the lower bound, of the misfit versus stress drop curve is at 7 MPa (Fig. 5a), which is consistent with our earlier results.
(a) L-curve illustrating the trade-off between the average stress drop (in MPa) and the misfit (cm). The black dots denote the models used which are treated as data. A fourth order 2-D spline curve used to fit the data is in blue. The corner, or the point of maximum curvature of the blue line, is located near |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| = 7 MPa (see the text for details). (b) Fault roughness (Laplacian of slip) versus RMS misfit (cm). The dashed lines represent the values of the Laplacian and the RMS of the model with |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| = 7 MPa.
As pointed out by Adams et al. (2017), there is a positive correlation between |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| and the fault roughness defined using the Laplacian of the inverted slip distribution. For a comparison, in Fig. 5(b) we show the RMS misfit and fault roughness of every model used in Fig. 5(a). The two dashed lines denote the values of these two parameters associated with the model with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$|, respectively. Visually these models also fall in an L-shape curve and the location of the model with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| is near the point of maximum curvature of this inferred L-shape curve (Fig. 5b). Hence, even using the slip roughness as a constraint, the model with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| is still one of its optimal solutions. This similarity implies that the solutions using conventional finite fault methods will have a similar |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| that is close to |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| if all of them use similar data sets. This is consistent with the previous analysis of the published slip models (Causse et al.2014).
Fig. 6(a) shows the slip distributions of the models with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}} = {\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$|, hereafter referred to as ‘Meb-7MPa’. It has a cumulative seismic moment of 7.55 × 1020 Nm (Mw 7.8). The inverted peak slip is 8.67 m. The weighted average fault slip amplitude |$\bar{D}$|, defined as |$\bar{D} = \frac{{\mathop \sum \nolimits_j \mathop \sum \nolimits_k D_{jk}^2}}{{\mathop \sum \nolimits_j \mathop \sum \nolimits_k {D_{jk}}}}\ $|, where Djk is the slip amplitude of subfault (j,k) (Ji et al.2002b), is 4.3 m (Fig. 7b). The slip distribution of this model is compact, as only 43 per cent of all subfaults have significant slip, which is here defined as slip larger than 0.87 m, 10 per cent of the peak slip. Most of the subfaults with significant slip are within the largest slip patch with an inverse triangular shape, which is outlined by the 10 per cent peak slip contour in white (Fig. 6a). Several small slip patches can also be seen at the edges of the fault plane, which are likely caused by errors in the observations.
Slip distributions for inversions with just the fixed average stress drop constraint, (a) and (b), and also with the additional maximum local stress constraint, (c) and (d). Panels (a) and (c) have |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| equal to 7 MPa, and (b) and (d) have |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| = 20 MPa. Contours outline the areas with slip (m) greater than 10 per cent of the peak slip value. Visible greater degree of heterogeneity in tests done with a higher |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$|, though the overall shape of the large slip patch is well resolved in both distributions.
(a) Map projection of the difference between the slip distributions of Meb-20 MPa and Meb-7 MPa. (b) Variation of weighted average slip (m) for subfaults with displacement greater than 10 per cent of the peak slip value versus stress drop (MPa). The vertical dashed line marks |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| = 7 MPa and the horizontal line shows the corresponding average slip value.
According to the values of the objective function, all models with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}} \ge$|7 MPa cannot be distinguished simply by their fit to the observations (Fig. 2a); all of them also have a similar seismic moment. Fig. 6(b) shows the slip distribution of one such representative: the model with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| equal to 20 MPa, which we will refer to as ‘Meb-20MPa’. Note that the edges of the asperities are again outlined using the same criterion as above, by contours of 10 per cent of the peak slip. Comparing Fig. 6(a) with Fig. 6(b), we can see that as |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| increases, the fine scale roughness intensifies. The peak slip of Meb-20 MPa is ∼15 m, and its weighted average fault slip |$\bar{D}$| is 5.7 m, 25 per cent larger than that of model Meb-7 MPa. On the other hand, the spatial extent and location of the main slip asperity remains relatively intact, despite the increased heterogeneity. For example, it can be seen that the shape of the large triangular asperity in model Meb-20 MPa spans nearly the same area as that of Meb-7 MPa (white contour line; Figs 6a and b). Fig. 7(a) illustrates the difference between Meb-7 and Meb-20 MPa. On the individual subfaults, the differences range from −5.85 to 9.32 m, which is notable considering the peak slip of model Meb-7 MPa is only 8.67 m. The mean of the differences is 0.002 m and the standard deviation is 1.67 m. Over the main slip asperity, regions of the slip distributions of the two models generally differ in the range of −4 to +4 m. Fig. 7(a) then highlights the difficulty to constrain the slip of individual subfaults using available geodetic data. Then, as shown in Fig. 7(b), the weighted average fault slip |$\bar{D}$| increases almost linearly from 4.3 to 8.5 m as |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| increases from 7 to 50 MPa.
Figs 8(a) and (b) show the static stress drop distributions of Meb-7 MPa and Meb-20 MPa, respectively. The heterogeneity of the slip distributions is clearly reflected in the stress distributions, with the most variability in the model with the higher stress drop. The subfaults with the large stress drop are also correlated with the subfaults with high slip. The peak stress drop of Meb-20 MPa is about two times greater than that of Meb-7 MPa.
Stress distributions for inversions using just the fixed average stress drop constraint, (a) and (b), and with the additional maximum local stress constraint, (c) and (d). Panels (a) and (c) have |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| equal to 7 MPa, and (b) and (d) have |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| = 20 MPa. Shear stress drop in the overall slip direction (rake angle of 96°) is shown. Visible greater degree of heterogeneity in the tests done with a higher |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$|, though the overall regions of high stress are well resolved in both distributions. Clearly, imposing a local maximum allowable stress drop forces a generally smoother and more widely distributed stress distribution, (c) and (d). In all cases, the respective slip distribution (shown in Fig. 6) is contoured and superimposed on the stress distributions here.
Furthermore, tests were conducted that relaxed the moment constraint, that is, |${\lambda _1}$| was set to 0.0001 rather than 0.1, and the inverted seismic moment increased by 8 per cent (6.98 × 1020 Nm versus 7.55 × 1020 Nm). The misfit value also decreased slightly from 2.71 to 2.69. The only notable difference was that when the moment constraint was less strict, the peak slip was slightly more: 9.3 versus 8.7 m.
5.2 Can we constrain |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ max}}$| of the 2015 Gorkha earthquake?
We have explored the possibility of constraining the maximum of |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| and found that under two circumstances |${E_{\mathrm{ fit}}}$| will increase again when |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| becomes large enough. First, we find that the pattern of the |${E_{\mathrm{ fit}}}$| versus |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| curve shown in Fig. 4(a) depends on the subfault size. For comparison, we have repeated the above analysis with models using 8 km by 8 km subfaults, that is, four times larger in area. The number of inverted free parameters is subsequently reduced by a factor of 4. As shown in Fig. 4(b), the value of the misfit function |${E_{\mathrm{ fit}}}$| decreases quickly as |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| increases from 2 MPa and reaches the minimum when |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| equals 7–8 MPa, consistent with the above analyses using 4 km by 4 km subfaults. However, as expected the minimum value of |${E_{\mathrm{ fit}}}$| is 3.395, larger than the analysis using 4 km by 4 km subfaults. When |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| further increases, the value of |${E_{\mathrm{ fit}}}$| starts to increase as well, particularly when |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| exceeds 18 MPa (Fig. 4b). Hence, when using large subfaults, the |${E_{\mathrm{ fit}}}$| versus |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| curve seems to suggest that the upper bound of |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| can also be constrained using geodetic data. Note that using a larger subfault size limits the shortest wavelength of the spatial variation of slip allowed, and subsequently, the shortest wavelength of the local stress drop variation. Adams et al. (2017) pointed out that the difficulty to constrain the upper bound of the energy-based average stress drop |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| is due to the difficulty of constraining the shortest wavelength of the local stress drop variation with surface data. Thus, this result with static data is an expected artefact.
Figs 6(c) and (d) show the inverted slip distributions using eqs (7) and (8) as the objective function, with a |${\tau _{\mathrm{ target}}}$| of 7 and 20 MPa, respectively. We refer to them as ‘Mls-7MPa’ and ‘Mls-20MPa’, respectively. It can be seen that the constraint to the local stress drop has negligible impact to the general shape of the slip asperities, which are similar to the first set of models (Figs 6a and b). The peak slips for these models are similar to the first two models at 8.8 and 10.4 m, respectively. Also, both Mls-7 MPa and Mls-20 MPa exhibit fewer patches of high slip, that is, considering only the ruptured area, the asperity itself has a more homogeneous slip distribution, than Meb-7 MPa and Meb-20 MPa, respectively. Overall, these four models do have similar slip extensions and the asperity boundary in all cases considered here are well constrained.
Regardless of the promise to constrain the upper bound of |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| with additional information of the maximum local stress drop |$\Delta {\tau _{\mathrm{ max}}}$|, readers should be aware that |$\Delta {\tau _{\mathrm{ max}}}$| mentioned here is different than what has been found in the laboratory. The |$\Delta {\tau _{\mathrm{ max}}}$| is the maximum stress drop across a fault length proportional to the subfault size, that is, on the order of kilometres. In contrast, the rock units used in laboratory experiments are generally on the order of metres or less (e.g. Mclaskey et al.2014). The |$\Delta {\tau _{\mathrm{ max}}}$| should then be less than the maximum local stress drop found in laboratory experiments.
6 DISCUSSION
Our previous study of the 2014 MW 7.9 Rat Islands earthquake (Adams et al.2017) revealed that using teleseismic waveforms, only the lower bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| could be constrained. We interpreted it as a result of the limited resolution to the fine scale roughness of the fault slip. In this study, we reach the same conclusion for the 2015 MW 7.8 Gorkha, Nepal earthquake with its excellent near-source geodetic observations. We further illustrate that the upper bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| can be accessed only when the highest local stress drop is known. Also, we must caution readers about the potential for a pseudo upper bound caused by using too large a subfault size.
6.1 Preferred slip model of the 2015 Gorkha event and its uncertainty
All models with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}} > {\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| match the data indistinguishably. Considering the large differences amongst these models (e.g. Fig. 7a), it is not meaningful to discuss the slip distributions without additional constraints to |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$|. As mentioned above, |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| is related to the ‘apparent’ available strain energy (|${\rm{\Delta }}{W_0})$|, that is, |${\rm{\Delta }}{W_0} = \frac{1}{2}\ {\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}P,$| where P denotes the seismic potency (Kanamori & Rivera 2006). |${\rm{\Delta }}{W_0}$| is the portion of the total strain energy responsible for creating a new fault plane and radiating seismic waves. Since the seismic potency is often well constrained, the solution with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}} = {\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| is then the model with the minimum available strain energy |$( {{\rm{\Delta }}{W_0}^{\mathrm{ min}}} )$| that also matches the data.
We have conducted 10 additional inversions with different initial models but with the same target stress drop of |$\ \ {\tau _{\mathrm{ target}}} = {\rm{\ }}{\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| (7 MPa). While what is shown in Figs 6(a) and 8(a) can be viewed as one representative model, here Figs 9(a) and (b) show the mean slip distribution of all of these slip models and the corresponding standard deviations for individual subfaults, respectively. We find that besides the similar seismic moment and |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$|, the weighted average slip, |$\bar{D}$|, defined above is also incredibly stable, changing negligibly from 4.30 to 4.41 m among these models. The variation of slip on individual subfaults is still notable but much smaller than what is shown in Fig. 7(a). The average standard deviation of slip is only 0.42 m (5 per cent of the peak slip), though on individual subfaults the standard deviation of slip can be up to 1.5 m (18 per cent of the peak slip). However, amongst these models the peak slip varies from 8.5 to 11.7 m, suggesting that peak slip is not a well-constrained variable.
(a) Average slip distribution for 10 different models, varying only the choice of the initial model, with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}} = \ 7\ \mathrm{ M}\mathrm{ Pa}$|, and (b) distribution of the standard deviation of the 10 models. Contours outline the area where greater than 10 per cent of the peak slip occurs. (c) Average stress distribution of the 10 models and (d) the distribution of the standard deviation of the 10 different models used to create the average model. Contours outlining the slip distribution in (a) are superimposed on the stress distributions in (c) and (d).
Given the small standard deviation, it is unsurprising that the mean of these slip models (Fig. 9a) is visually identical to what is shown in Fig. 6(a) and has a cumulative seismic moment of 7.55 × 1020 Nm (Mw 7.8). The peak slip of this mean slip distribution is 8.36 m, and also has a compact slip distribution as only 389 (44 per cent) subfaults have significant slip, that is, slip larger than 10 per cent of the peak slip. Most of the subfaults with significant slip are within the aforementioned large slip patch with an inverse triangle shape (Fig. 7a). It extends 152 km along-strike and 56 km downdip and the corresponding depth extension is limited to 9.3–15.2 km.
Figs 9(c) and (d) further show the average stress distribution of the models with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| and its respective standard deviation. The average stress drop varies from −13.3 to 32.5 MPa over the fault (Fig. 9c) but has large uncertainty. The standard deviations on individual subfaults are up to 12.6 MPa with a mean of 3.7 MPa. Considering such a large uncertainty, it is not surprising that in Fig. 8(c) even when we set the maximum local stress drop to 25 MPa (Fig. 8c), the inverted model still can explain the data equally well. In other words, in the case where a local maximum stress drop is enforced, the distribution of the average stress drop across the total fault plane does not vary drastically, but there are more individual patches of subfaults with high stress drop, as |${\overline {{\rm{\Delta \tau }}} _{\rm{E}}}$| is still being forced to τtarget. In Fig. 9(c), the circles denote the epicentres of 379 M > 4.5 aftershocks that occurred between 2015 April 25 and 2015 May 12 and were relocated using the local network data (Adhikari et al.2015). Although some of them may have occurred outside the fault plane that hosted the mainshock, it can be seen that in map view that most of them are located in regions of negative shear stress drop, as expected (Mendoza & Hartzell 1988).
The group of slip models with |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| are our preferred models. These solutions are generally consistent with previous analyses using seismic and geodetic data (e.g. Avouac et al.2015; Galetzka et al.2015; Wang & Fialko 2015; Yue et al.2016) in terms of the location of the dominant asperity and its along-strike and downdip extension. However, we admit that the peak slips of our models are considerably higher than previous analyses, which are less than 6.5 m (Avouac et al.2015; Galetzka et al.2015; Wang & Fialko 2015; Yue et al.2016). The cause of this discrepancy is not unique. Differences in fault plane depth, subfault size, the weight of the smoothing constraint and the choice of data used are potential reasons. As we mentioned above, the value of the peak slip is not well constrained. Nevertheless, in comparison with the previous Laplacian smoothing constraints, it appears that our new approach tolerates more small-scale variations.
6.2 |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| and the energy budget
Our estimate of the seismic potency due to the 2015 MW 7.8 Gorkha earthquake is 2.28 × 1010 m3 (Meb-7 MPa). This value is close to the value inferred from GCMT, 2.54 × 1010 m3, as well as the potency from the USGS W-Phase, 2.01 × 1010 m3 (using the same rigidity of 3.3 × 104 MPa). Just considering |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| = 7 MPa, the minimum ‘apparent’ available strain energy, |$\Delta W_0^{\mathrm{ min}}$|, during this earthquake is 7.99 × 1016 J. Such an estimate is also data dependent. Using their slip model constrained by teleseismic body waves, Lay et al. (2016) estimated the |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| of the Gorkha earthquake to be 3 MPa, a factor of two smaller than our |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| value. We suspect that this discrepancy simply reflects the fact that the near-fault static observations better resolve the small-scale spatial variations of the fault slip than the teleseismic body waves, that is, |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| estimated using teleseismic data might be systematically lower. Consequently, an underestimation of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| will lead to an underestimation of ΔW0.
The radiated seismic energy of the Gorkha earthquake has been estimated using teleseismic data by multiple researchers (see a review in Lay et al.2016), varying within a factor of 2 from 7.1 × 1015 to 12.5 × 1015 J (e.g. Lay et al.2016). The corresponding ratio of ER and M0 is 0.8–1.5 × 10−5 and the apparent stress drop is 0.3–0.5 MPa. Both are compatible with the global average values of large thrust earthquakes (ER/M0, 1.06–1.81 × 10−5; Convers & Newman 2011; Ye et al.2016). Using the |$\Delta W_0^{\mathrm{ min}}$| mentioned above, the inferred seismic radiation efficiency (ηR = ER/ΔW0) is only 0.09–0.15. This estimate of ηR is close to its upper bound, that is, |$\eta _R^{\mathrm{ max}}$|, unless ER was significantly underestimated during the previous studies (e.g. Lay et al.2016). It is considerably lower compared with estimates of other large earthquakes. Ye et al. (2016) systematically analysed the finite fault solutions of 114 Mw > 7 earthquakes constrained by teleseismic body waves. The mean of their estimates of ηR is 0.34 (logarithmic) or 0.49 (linear) (Ye et al.2016). However, as mentioned above, ηR estimates in Ye et al. (2016) might be systematically higher because of the underestimated ΔW0 using teleseismic data.
Using |$\frac{{{V_\mathrm{ R}}}}{{{V_\mathrm{ s}}}}$| = 0.8–0.93 and approximating |$\varepsilon$| as the L/W ratio mentioned above, the predicted |${\eta _\mathrm{ R}}$| using eq. (10) is 0.09–0.12, consistent with our estimation. However, readers should be aware that in the model of Kikuchi & Fukao (1988), the slip within the main asperity varies smoothly. Then, the heterogeneous rupture process during the Gorkha earthquake suggested by the inverted slip (Figs 6 and 9) and stress drop distributions (Fig. 7) will radiate more seismic energy than a quasi-static elliptic crack rupture. Nevertheless, our result suggests that the earthquake rupture could propagate quickly and still dissipate most of the available strain energy in fracturing the fault plane. Hence, investigating the upper bound of |${\eta _\mathrm{ R}}$| provides an additional but interesting angle to study the earthquake rupture process.
6.3 Effect of subfault size
The choice of the subfault size is one of the most important a priori conditions affecting the source inversion (e.g. Hartzell & Helmberger 1982; Beresnev 2003). While small subfaults can provide a more detailed description of the earthquake rupture process, limitations in the station distribution and bandwidth of the data that can be modelled forces modellers to use a large subfault size and to sometimes even ignore slip in one dimension to decrease the number of unknowns (e.g. Bouchon & Vallee 2003). For the inversion using seismic data, spatial changes less than the smallest wavelength within a certain frequency range of the seismic waves for the analysis cannot be constrained by the data (e.g. Sekiguchi et al.2000). Thus, in some instances, this has been considered as a rough guide for the choice of the subfault size (e.g. Shao & Ji 2012). For geodetic inversions, efforts have been made to use irregular subfault sizes according to data resolution (e.g. Fialko et al.2001; Lohman & Simons 2005; Page et al.2009). Alternatively Yabuki & Matsu'ura (1992) proposed to represent the spatial distribution of each slip component as the linear combination of a finite number of cubic B-spline functions. However, the characteristic length of the basis functions needs to be pre-assigned. This then affects the finest spatial resolution, similar to the subfault size.
When we discuss the uncertainty of inverted parameters, the effect of the subfault size cannot be ignored. As stated in the methods section, the stress change is estimated at the centre of each subfault. Thus it can only represent a low-pass filtered version of the target stress distribution. The maximum wavenumber is |$\frac{1}{{2\Delta x}}$|, where |$\Delta x$| is the subfault size. In other words, solutions in the model space cannot represent the spatial variation of stress on a scale smaller than two times the subfault size. Such a limitation will inevitably affect the stress estimation of the inverted result (e.g. Dreger et al.2007), and as it is an a priori condition its impact can only be explored by comparing the final results using different subfault sizes. The comparison shown in Fig. 4(b) suggests that the low bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| defined here is stable in terms of the subfault size but the upper bound of |${\overline {{\rm{\Delta }}\tau } _E}$| is not. Readers should also be aware that such a conclusion is for the relatively ‘small’ subfault sizes of 4 to 8 km. If we instead use a much larger subfault size, even the low bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| could be affected. For a given earthquake, the ‘small’ or ‘large’ subfault sizes mentioned here are not absolute values. Instead it is correlated with not only the data used but also the earthquake itself. A stable |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| can, therefore, be used as one of the conditions during the search for an appropriate subfault size for the given observations.
The uncertainty of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| has received considerable attention in the source inversion community. Duputel et al. (2015) estimated the uncertainty of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| during the 2014 Mw 8.1 Iquique earthquake by analysing the ensemble of accepted solutions obtained during their Bayesian inversion. They found that the posterior probability density of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| has a peak around 10 MPa and is negligible for |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| < 5 MPa and |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| >25 MPa. In principle all accepted solutions during Bayesian inversions should have good misfit values, and to some extent, the probability of the solutions are correlated with the misfit values. The low posterior probability density for models with smaller |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| is then consistent with our result, and |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| is close to 5 MPa. At first glance, the low posterior probability density for the models with large |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| is inconsistent with our result. However, as we mentioned above (Fig. 4a), such an apparent |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ max}}$| of about 25 MPa can simply be an artificial result of using large subfaults. Besides, it is noteworthy that the inferred posterior probability density of one model also depends on the topology of the model space near this optimal solution. For example, even the two arbitrary optimal solutions have the exact same misfit value, the one associated with a broad and flat ‘basin’ type minimum in terms of the misfit value has a higher posterior probability density than the other with a narrow ‘valley’ type minimum. Therefore, it is also possible that the latter scenario is associated with the optimal models with large |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$|. Nevertheless, the low posterior probability density for the models with large |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| found by Duputel et al. (2015) cannot be used as unique evidence against our statement that the |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ max}}$|of an earthquake cannot be solely constrained by the surface geophysical observations. Instead, it might provide supplementary information on this interesting source property.
6.4 Using |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| as a new inversion constraint
We end our discussion by outlining the finite fault inversion strategy using the ‘minimum average stress drop constraint’. After the routine data and Green's functions preparation, this new inversion procedure starts with a grid search for |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$|. The subfault size will then gradually be reduced until |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| becomes stable and the flat pattern shown in Fig. 4(b) can be observed in the trade-off curve between the misfit and |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$|. A method such as the L-curve algorithm mentioned above is then applied to more rigorously determine the point of maximum curvature. The corresponding stress drop thus becomes |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| of this earthquake, which is data dependent. To explore the epistemic uncertainty in the inverted slip and stress drop distributions, multiple inversions with |${\tau _{\mathrm{ target}}} = {\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}^{\mathrm{ min}}$| are then conducted, with the only difference in this part being the choice of the initial model. Other types of perturbations, such as bootstraping the data (Custodio et al.2009) or allowing a perturbation in the alignment of data and synthetics (Hartzell et al.2007) when using seismic data, could be included in a future study. It is of interest to note that the new procedure has no parameter that needs to be assigned empirically. However, as hundreds of inversions need to be conducted, this procedure is time consuming. Nevertheless, all individual inversions are independent and identical in a computational perspective. The entire procedure is then easily made parallel and can take full advantage of desktops or clusters with a modern multiple core CPU or GPU.
7 CONCLUSION
Resolving the state of stress on the fault is a key component in understanding the physics of earthquakes, and it has important implications to the earthquake energy budget. As a continuation of our companion work (Adams et al.2017), we found that only the lower bound of the energy-based average stress drop |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| can be constrained using teleseismic data, and here we demonstrate that this argument holds even with excellent near-source geodetic data for the 2015 MW 7.8 Gorkha earthquake. The lower bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| or |${\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}} ^{\mathrm{ min}}}$| during this earthquake is ∼7 MPa. We also show numerically that the upper bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| may also be constrained if an extra constraint to the maximum local stress drop is available (e.g. via laboratory experiments). At the same time, readers should be cautious about a potential pseudo-upper bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| resulting from the use of a large subfault size to represent the source.
The lower bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| leads to the lower bound of the apparent available strain energy ΔW0 (7.99 × 1016 J) and an upper bound of the seismic radiation efficiency (ηR), though the latter is also affected by the large uncertainty of the observed seismic radiated energy. Both of these are important parameters in earthquake energy budget analyses. The estimate of |$\eta _\mathrm{ R}^{\mathrm{ max}}$| is 0.09–0.15, despite the high rupture velocity (VR = 83–94 per cent of the shear wave speed). We point out that the low ηR and high VR can be reconciled by considering the aspect ratio of the dominant slip patch. The earthquake rupture with large aspect ratio could propagate quickly but still dissipate most of the available strain energy in fracturing the fault plane. Finally, we have introduced a new, more physical constraint to smooth the inverted fault slip during finite fault inversions, the ‘minimum average stress drop constraint’, which is based on finding the slip and stress distribution with the minimum available strain energy of the earthquake.
ACKNOWLEDGEMENTS
We thank Prof. M. Matsu'ura and an anonymous reviewer for their constructive comments; and Ralph Archuleta for many fruitful discussions and suggestions. The unwrapped InSAR line-of-sight deformation data was processed and distributed by the Satellite Geodesy research group of UCSD (http://topex.ucsd.edu/nepal/) using the data of ALOS-2 satellite (operated by the Japan Aerospace Exploration Agency). The coseismic displacements at near-fault GPS stations deployed by the Caltech Tectonics Observatory (e.g. Avouac et al.2015) were processed by Advanced Rapid Imaging and Analysis Center for Natural Hazards at Jet Propulsion Laboratory (JPL). This study is supported by grant EAR-1215769 and the Southern California Earthquake Center (SCEC; contribution 7062), which is funded by NSF Cooperative Agreement EAR-0109624 and USGS Cooperative Agreement 02HQAG0008. Finally, Jinlai Hao was supported by the National Natural Science Foundation of China (grant 41541033).
