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

Assessment of seismic hazard relies on estimates of how large an area of a tectonic fault could potentially rupture in a single earthquake. Vital information for these forecasts includes which areas of a fault are locked and how the fault is segmented. Much research has focused on exploring downdip limits to fault rupture from chemical and thermal boundaries, and along-strike barriers from subducted structural features, yet we regularly see only partial rupture of fully locked fault patches that could have ruptured as a whole in a larger earthquake. Here we draw insight into this conundrum from the 25 April 2015 Mw 7.8 Gorkha (Nepal) earthquake. We invert geodetic data with a structural model of the Main Himalayan thrust in the region of Kathmandu, Nepal, showing that this event was generated by rupture of a décollement bounded on all sides by more steeply dipping ramps. The morphological bounds explain why the event ruptured only a small piece of a large fully locked seismic gap. We then use dynamic earthquake cycle modeling on the same fault geometry to reveal that such events are predicted by the physics. Depending on the earthquake history and the details of rupture dynamics, however, great earthquakes that rupture the entire seismogenic zone are also possible. These insights from Nepal should be applicable to understanding bounds on earthquake size on megathrusts worldwide. INTRODUCTION 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 largescale structural discontinuities. Examples include the 2014 Mw 8.2 Iquique, Chile (e.g., Ruiz et al., 2014), and the 2007 Mw 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 Mw 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).


Introduction
This supplementary information contains details about the dynamic earthquakecycle modeling procedure (Text S1); dynamic model assumptions and limitations (Text S2); the combined multiple tracks of InSAR data used in the study (Table S1); and forward model prediction and residuals between the data and the forward model prediction based on our geodetic inversion modeling ( Figure S1 to S4); and the Coulomb stress changes on the megathrust from our preferred slip model ( Figure S5).

Text S1.
We simulate the slip evolution on an isolated segment of the Main Himalayan Thrust (MHT) using a rate-and-state friction framework. We employ the boundary integral method with radiation damping (Rice et al., 2001;Kato, 2003). As with the inversion, we discretize the geometry of the MHT using triangular elements. We obtain the stress kernels using the analytic formulation of Nikkhoo and Walter (2015). We use uniform friction properties with the static friction coefficient µ 0 = 0.6 at steady state with the reference velocity V 0 = 10 -6 m/s, a= dµ/d(ln V)= 10 -2 , the effective confining pressure σ= 300 MPa, and the characteristic weakening distance L= 7.5 cm (µ is the frictional resistance and V is the slip velocity).
We define the region that could potentially rupture as an area larger than the rupture patch of the 2015 M w 7.8 Gorkha earthquake, extending from the surface to 20 km depth and about the length of the Kathmandu Klippe ( Figure 2B). This region is modeled with velocity-weakening properties, with (a-b)= dµ ss /d(ln V)= -4x10 -3 , where µ ss is the frictional resistance at steady state. Outside this region, the fault is modeled with the velocity-strengthening properties of (a-b)= 4x10 -3 . We load the fault at a uniform and constant rate of V pl = 0.02 m/yr.
In the simulation, we obtain a series of earthquakes of varying magnitudes and with ruptures covering sometimes a fraction of the velocity-weakening patch and some other times the entirety of the patch. The complexity of the slip history emerges spontaneously from the dynamics of the system. In the context of our simulation, the complexity of the slip history comes from the morphology of the megathrust, not from heterogeneous initial conditions or spatial distribution of friction properties. (If the fault were perfectly planar, such a distribution of frictional properties would lead to a single type of earthquakes with a characteristic size and a single recurrence time and both recurrence times and rupture slip would be predictable.) The simulated fault history contains events with characteristics that resemble those of the 2015 M w 7.8 Gorkha earthquake. The Gorkha-like ruptures extend through the deep décollement but are confined by the two ramps. Other single-segment ruptures occupy the upper décollement or the topmost ramp, resulting in earthquakes of M w <8. Occasionally, ruptures spread from the surface to the bottom of the seismogenic zone in a through-going, multi-segment event, resulting in a great earthquake with M w >=8. The simulation is only intended to explain the role of geometry on the complexity of earthquake dynamics in the context of a long sequence of earthquakes. In reality, strong heterogeneities in composition, rheology and other physical properties make the earthquake history even more complex. As we present a simulation with a simple distribution of dynamic frictional properties, some important differences with the Gorkha setting arise. For example, the simulated afterslip that follows the M w 7.7 event, which is similar to the Gorkha earthquake, has hundreds cm of afterslip downdip of the coseismic rupture, presumably more than expected for this event.
In this study, we highlight one particular combination of effective normal stress and dynamic friction properties to show that fault morphological barriers can control the extent of earthquake ruptures and introduce complexity in earthquake cycles. Admittedly, not all distributions of dynamic friction properties produce similar results. Our study merely provides a proof of existence of a combination of frictional parameters and geometrical properties spontaneously giving rise to segmentation. We do not aim to reproduce every aspect of the 2015 Gorkha event. Theoretical insight about the propagation of elasto-dynamic rupture across fault bends can be found in the work of Poliakov et al. (2002) and Rice et al. (2005).

Text S2.
We model earthquake cycles using the rate-and-state framework and the boundary integral method under the radiation-damping approximation. The method allows us to evaluate earthquake cycles on geometrically complex, non-planar, faults. This is accomplished by damping of the seismic waves, ignoring the stress concentration at the rupture tip induced by the wave-mediated stress transfer. The numerical simulation might produce biased estimates of rupture size, recurrence times, and other rupture characteristics (Thomas et al., 2014).  Figure S1. Footprint of InSAR data considered in this study, including from the Sentinel-1A, RADARSAT-2, and ALOS-2 satellites. The beach balls represent the moment tensors of the mainshock and largest aftershock. The coloured circles correspond to aftershocks (same as in Figure 1). Second panel   Figure S4. Examples of estimated slip distributions obtained using different relative weights for data input to the slip inversion. The GPS data prefer a little more shallow slip, but in general the distributions remain similar regardless of data weight.  Figure S5: Coulomb stress changes on the megathrust from our preferred slip model, assuming a coefficient of friction of µ 0 =0.6. The dash lines show the location of the shallow (8-10 km) and deep ramp (16 km), respectively. The coloured dots correspond to aftershocks (same as Figure 1). The beach balls represent the moment tensor of the mainshock and largest aftershock.