The mechanism of partial rupture of a locked megathrust: The role of fault morphology

Earthquake magnitudes are directly related to the spatial extent of the fault patch that ruptures and the amount of slip on that patch. The level of slip available to generate an earthquake is determined by the loading rate, the level of interseismic coupling, and the history of past ruptures. Seismic hazard estimates, therefore, rely on assessments of the potential spatial extent and magnitude of rupture patches. Geodetic data can be used to estimate coupling maps, which give an estimate of the spatial patterns of fault locking at present (Chlieh et al., 2008; Ader et al., 2012; Stevens and Avouac, 2015). Analyses of seismological data and numerical modeling have been used to assess the role of thermal and lithological boundaries in controlling the extent of possible rupture, as well as the mode of failure likely at different depths (Oleskevich et al., 1999; Kaneko et al., 2010; Lay et al., 2012).

This simple picture, however, does not address the fact that we regularly see only partial rupture of fault patches that were fully locked and at a late stage in the seismic cycle with no documented evidence for large-scale structural discontinuities. Examples include the 2014 M_w 8.2 Iquique, Chile (e.g., Ruiz et al., 2014), and the 2007 M_w 8.4 Bengkulu, Sumatra, earthquakes (Konca et al., 2008). How do we explain the occurrence of moderate to large earthquakes in the middle of patches that have, by this simple model, been forecast to produce great earthquakes?

The 2015 M_w 7.8 Gorkha (Nepal) earthquake was one such conundrum. The magnitude of this event was considerably smaller than expected from forecasts of a great earthquake that were based on the fact that this area of the Main Himalayan thrust (MHT) is fully locked and at a seismic gap (Ader et al., 2012; Bilham, 2015; Stevens and Avouac, 2015). The rupture did not reach the surface, despite considerable evidence that neighboring ruptures have (Lavé et al., 2005; Sapkota et al., 2013; Bollinger et al., 2016). Why did this event not cascade into complete failure of the seismic gap, all the way to the surface, generating a great earthquake?

The answer lies with the morphology of the fault at the local scale (Grandin et al., 2015; Hubbard et al., 2016). In this study we use the three-dimensional (3-D) fault morphology of Hubbard et al. (2016) for inversion of geodetic data to estimate the extent of rupture, showing that it was bounded by ramps on all sides, and for dynamic models of slip evolution that explain the role of this morphology in generating a complex earthquake history with single and multisegment ruptures.

The model for fault morphology was determined by studying the surface geology, topography, and seismicity of the area (Fig. 1; Hubbard et al., 2016). While broadly compatible with other geophysical data (e.g., Duputel et al., 2016), the geometry differs by including two blind ramps updip and downdip of the shallow-dipping décollement that generated the earthquake. The deeper ramp, at depths of ∼16–30 km, is required to produce the large Gorkha-Pokhara anticlinorium, which extends across central Nepal. The shallow ramp, which extends from depths of ∼8–10 km, explains the surface stratigraphy at this location. In some places, the two ramps join together, resulting in a single large ramp; this is the case at both ends of the Kathmandu klippe (Fig. 1).

We use the data from 18 GPS stations from the network run by the California Institute of Technology (Pasadena, California, USA) and the Department of Mines and Geology, Nepal (Avouac et al., 2015; Galetzka et al., 2015). Details on our data processing procedure were given in Feng et al. (2015). Three GPS stations located close to the rupture region show a maximum southward movement reaching ∼1.8 m and uplift of ∼1.3 m. The remaining stations, located further from the source, show negligible displacements.

We used all available interferometric synthetic aperture radar (InSAR) data for this event from the Sentinel-1A (European Space Agency; an ascending path), RADARSAT-2 (Canadian Space Agency; a descending path), and ALOS-2 (Advanced Land Observing Satellite; Japan Aerospace Exploration Agency; ScanSAR descending paths 048 and 047, and ascending path 157) satellites (Fig. DR1 and Table DR1 in the GSA Data Repository¹). We processed the Sentinel-1A and RADARSAT-2 data using GAMMA (GAMMA Remote Sensing Research and Consulting AG, Gümligen, Switzerland; Scheiber and Moreira, 2000; Wegmuller et al., 2015). For the ALOS-2 data, we used GMTSAR (http://topex.ucsd.edu/gmtsar/; Sandwell et al., 2011) to generate the unwrapped interferogram (Lindsey et al., 2015).

To avoid redundancy of the InSAR data due to the high degree of spatial correlation, and to achieve an efficient inversion, we used a quadtree downsampling scheme (Simons et al., 2002) to extract the deformation pattern with a limited number of data points. Our inversion is therefore based on ∼9000 data points concentrated mostly over the rupture patch. We used triangular elastic half-space dislocation Green’s functions (Nikkhoo and Walter, 2015) to capture the fault geometry, and a linear least-squares scheme to estimate fault slip with positivity constraints and smoothing regularization. The parameter controlling the level of spatial smoothing was determined using Akaike’s Bayesian information criterion (Fig. DR3; Yabuki and Matsuura, 1992). The model and corresponding residuals for InSAR are shown in Figure DR2. The relative weight between InSAR and GPS data was chosen by trial and error so that fit to both data sets is satisfactory (Fig. DR4).

We find that the slip extended along strike on the shallow décollement for ∼130 km (Fig. 2A). The rupture area was bounded at either end by pinch points formed by the two-ramp system merging into a single ramp at either end of the Kathmandu klippe (Fig. 2A). Similarly, it was bounded updip and downdip by two ramps. The hypocenter was actually located on the deeper ramp; the intervening décollement appears to have provided a prime region for rupture expansion. Because the rupture was confined to the décollement, previously published slip models (e.g., Avouac et al., 2015) that were based on a single planar fault appear relatively similar to ours. We note that published patterns of coupling (Ader et al., 2012; Stevens and Avouac, 2015) also seem to correspond to the fault morphology (Fig. 2A). The rupture increased the Coulomb stress on the updip section of the Himalayan megathrust by a few megapascals, bringing this section closer to failure (Fig. DR5).

An important question is the extent to which the MHT at this location sometimes ruptures in a piecemeal fashion with ruptures bounded by morphological barriers, or sometimes cascades into complete rupture to the surface, generating a great earthquake. In the absence of a comprehensive historical earthquake catalogue, we look for insight from fault mechanics, and investigate the seismic potential of the Kathmandu section of the megathrust using dynamic models of slip evolution in the framework of rate and state friction (Lapusta et al., 2000). We discretize the geometry of the MHT using triangular elements, and the associated stress kernels are obtained using the analytic formulation of Nikkhoo and Walter (2015). We consider a simple scenario where the seismogenic zone extends from the surface to ∼20 km depth with uniform dynamic friction properties. The only complexity of the model is in the fault morphology (see the Data Repository). Our modeling shows that some large (M_w > 7) ruptures can be confined between two ramps, or on the topmost ramp. However, these relatively smaller events can also be preceded or followed by rupture of the entire seismogenic zone in a single, great (M_w > 8) earthquake (Fig. 2B). Our models show areas of high stress accumulation at ramp boundaries, which often initiate ruptures. These ruptures may or may not cascade into larger events, depending on the slip history and current stress status. The rupture sequences are unpredictable in terms of the location, time, and amount of slip. However, the sum of events will close the slip budget on the fault. The scarcity of throughgoing ruptures in our simulations may explain why megathrust ruptures do not often reach the surface, even in cases where the shallow fault is velocity weakening and capable of sustaining seismic ruptures. Overall, the nonplanar geometry introduces complexity in the dynamics of rupture propagation and heterogeneity in the slip distribution of earthquakes that make single events inadequate representatives of the full range of events.

Based on our results, we propose that both single-segment and multisegment ruptures are possible near Kathmandu (Fig. 3). We note that the A.D. 1833 M_w 7.6–7.7 earthquake had striking similarities to the 2015 earthquake; it also generated shaking in Kathmandu, and did not rupture to the surface (Bollinger et al., 2016). It seems plausible that both events ruptured the same confined décollement. However, the occurrence of earthquakes like the 1833 M_w 7.6 and the 2015 M_w 7.8 do not preclude the possibility of a great earthquake on this segment.

Other parts of the MHT have historically behaved differently from this section, with evidence for great M_w > 8 events rupturing all the way to the surface (Fig. 2A; Sapkota et al., 2013). The Gorkha section may be relatively unique in having such clear morphological constraints. However, we note other published inferences of geometrical controls on historical earthquake ruptures; for example, lateral ramps may have played a role in segmentation of the 1905 M_w 7.8 rupture (Middlemiss, 1910; Molnar, 1987).

The effect of fault morphology on earthquake size is well known in strike-slip systems (Wesnousky, 1988, 2006; Klinger, 2010), where the bends are considerably easier to map than on megathrusts. This idea is also grounded in fault mechanics (Rice, 1980; King, 1986; Kame et al., 2003). For strike-slip systems, there is evidence that fault morphology can influence not only earthquake size, but also rupture behavior. For example, particularly straight segments may be more prone to supershear rupture velocities (e.g., Bouchon et al., 2010). Local megathrust fault morphology may therefore be an additional clue to understanding rupture behavior within the broad frictional regimes outlined by Lay et al. (2012). In Figure 3 we illustrate the idea that fault morphology may result in asperities (within the broad seismogenic zone) that have a different behavior to their neighbors and a tendency to create along-strike and downdip segmentation.

In our quest for more predictive earthquake models, these results suggest that gaining detailed images of fault morphology may be the key to better understanding fault segmentation, the variability of earthquake size, and potentially the rupture style. Collaboration between structural geologists and groups modeling rupture propagation may be essential to the next generation of earthquake forecasting.

This research was supported by the National Research Foundation (NRF) Singapore under its Singapore NRF Fellowship scheme (National Research Fellow Awards NRF-NRFF2010–064, NRF-NRFF2013–04, and NRF-NRFF2013–06) and by the Earth Observatory of Singapore (EOS) and the National Research Foundation Singapore and the Singapore Ministry of Education under the Research Centres of Excellence initiative. Advanced Land Observing Satellite data were provided by the Japan Aerospace Exploration Agency. Sentinel-1A data were copyrighted by the European Space Agency. We thank the Canadian Space Agency for providing RADARSAT-2 data. We are grateful to Jean-Philippe Avouac for useful discussions and for making the global positioning system data immediately available. This is EOS Paper 108. The simulation code used in this study is available at http://bitbucket.org/sbarbot/.