Seismogenic Potential of the Main Himalayan Thrust Constrained by Coupling Segmentation and Earthquake Scaling.

Recent studies have shown that the Himalayan region is under the threat of earthquakes of magnitude nine or larger. These estimates are based on comparisons of the geodetically inferred moment deficit rate with the seismicity of the region. However, these studies did not account for the physics of fault slip, specifically the influence of frictional barriers on earthquake rupture dynamics, which controls the extent and therefore the magnitude of large earthquakes. Here we combine an improved probabilistic estimate of moment deficit rate with results from dynamic models of the earthquake cycle to more fully assess the seismogenic potential of the Main Himalayan Thrust (MHT). We propose a straightforward and efficient methodology for incorporating outcomes of physics-based earthquake cycle models into hazard estimates. We show that, accounting for uncertainties on the moment deficit rate, seismicity and earthquake physics, the MHT is prone to rupturing in M w 8.7 earthquakes every T > 200 years.

Introduction

Historical records show that very large earthquakes have struck the Himalayan arc (e.g., Bilham, 2019; Bollinger et al., 2014; Wesnousky, 2020). The 2015 M7.8 Gorkha (e.g., Avouac et al., 2015; Elliott et al., 2016; Galetzka et al., 2015), 1950 ∼M8.7 Assam (Chen & Molnar, 1977) and 1934 ∼M8.4 Nepal‐Bihar (e.g., Chen & Molnar, 1977, 1983; Sapkota et al., 2013) earthquakes are the most recent examples of such hazard. These large events occurred on the Indian‐Eurasian plate interface called the Main Himalayan Thrust (MHT). Paleoseismological studies of this megathrust have documented multiple >15 m slip events suggesting multiple M > 8.0 earthquakes over the last millennium (e.g., Kumar et al., 2010; Lave et al., 2005; Wesnousky et al., 2017). Constraining this hazard and evaluating the magnitude of the largest possible earthquake, hereafter referred to as maximum magnitude, is crucial in such a populous region.

In general, plate tectonics dictates that a region accommodates a certain amount of deformation in a given time. Most evaluations of the seismic potential of a region estimate how this deformation rate is apportioned among earthquakes of different magnitudes and recurrence intervals, using historical and instrumental seismic catalogs as a guide. Significant efforts have been put forward in this direction in California (Field et al., 2014), in Europe (Woessner et al., 2015) and in more local studies (e.g., Chartier et al., 2017). Some models incorporate aspects of earthquake physics such as multi‐fault ruptures or slip profiles within ruptures (e.g., Field et al., 2014), and a few attempts have been made to include other aspects (e.g., quasi‐static elastic interactions; Shaw et al., 2018), but fault physics are not usually directly incorporated. This is largely because physics‐based models would significantly increase the computational time and power required to derive large scale earthquake models such as the Himalayan arc.

One particular missing ingredient in such estimates is the influence of the spatial distribution of locked sections of faults. The locking degree of a fault can be kinematically inferred from geodetic data, and is then referred to as fault coupling, (e.g., Chlieh et al., 2011; Hashimoto et al., 2009; Savage, 1983). At each point on a fault, coupling ranges between 0, corresponding to a fault that slips aseismically at its long‐term deformation rate (accruing no strain), and 1, corresponding to a locked fault (Avouac, 2015). In this framework, locked regions of the fault accumulate moment deficit, an indicator of how much elastic strain accumulates within the crust. Following the elastic rebound hypothesis (Reid, 1910), the moment deficit accumulated during the interseismic period should be balanced on the long‐term average by the moment released by earthquakes and aseismic slip (Avouac, 2015; Molnar, 1979). Evaluating this “moment budget” can provide first‐order estimates of a regions seismic potential (e.g., Michel et al., 2018; Rollins & Avouac, 2019; Rong et al., 2014; Shen et al., 2007; Stevens & Avouac, 2016, 2017).

The spatial variability of the slip budget along faults is not enough to constrain seismic hazard. On top of modulating the moment budget, regions slipping aseismically influence the propagation and arrest of dynamic ruptures, controlling the size and thus the magnitude of earthquakes (Avouac, 2015; Béjar‐Pizarro et al., 2010; Kaneko et al., 2010; Pritchard et al., 2007). The physics of fault slip, controlled by factors such as rheological properties and the history of stress along the fault, influences the timing and size of earthquakes, and so must be accounted for in hazard estimates. For instance, it has been recently proposed that evaluating the ability of a fault to rupture using estimates of fracture energy would be a significant step towards incorporating the physics of fault slip in hazard estimates (Weng & Ampuero, 2020).

On the MHT, moment budget evaluation based on previous coupling models and modern seismicity catalogs concludes that there is potential for millenary earthquakes (Stevens & Avouac, 2016). However, the corresponding probability density function (PDF) of the maximum magnitude earthquake does not exclude events. Such events have not been observed or inferred on any fault on Earth and, if they were to follow the empirical scaling between earthquake moment and rupture area , (e.g., Kanamori & Brodsky, 2004), would require rupture areas larger than the entire MHT. In addition, a recent Bayesian evaluation of coupling along the MHT (which relaxes the strong assumptions required by previous studies) suggests a more heterogeneous coupling distribution than previously inferred, including in particular an along‐strike segmentation of the MHT with three segments that may act as aseismic barriers (Figure 1, Dal Zilio et al., 2020). These low‐coupling sections of the megathrust should in principle influence the size of earthquakes in the Himalayas, a hypothesis supported by historical earthquakes (Bilham, 2019; Dal Zilio et al., 2021).

Figure 1

Regional setting of the Main Himalayan Thrust (MHT). (a) Interseismic coupling map along the MHT (Dal Zilio et al., 2020). Black dashed lines indicate the position of the potential barriers to earthquake propagation. The dots represent the microseismicity since 1995 from the ANSS catalog. Dark green dots are earthquakes selected for the seismic potential analysis (<200 km from the MHT trench). Thick green lines indicate the spatial extent of historical large earthquakes ( since 1500; Bilham, 2019). The thin solid blue line delimits the bottom extent of the coupled zone, based on a coupling threshold of 0.3. The bottom left inset shows the probability density function of the moment deficit rate. (b) Gray lines indicate the along‐strike distribution of the mean coupling of each of the 360 000 coupling models (Supporting Information Text S1). The black line represents the mean value of the 360,000 coupling models. The red dashed line indicates the coupling threshold of 0.3 used in this study. (c) Width of coupled zone based on a coupling threshold of 0.3. Green shading indicates the width of the aseismic barriers for this particular case as an example.

In this study, we develop a computationally efficient method to include the influence of fault physics on seismic hazard estimates in the Himalayas and by extension in any setting. Our results suggest that earthquakes of magnitude comparable to that of the 1950 8.7 earthquake are the most likely maximum magnitude earthquake in the Himalayas.

Overview of Method

Our approach combines geodetic, seismological and physical constraints to probabilistically evaluate the seismic potential of the MHT. We convert the coupling model of Dal Zilio et al. (2020) into a probability density function (PDF) of total moment buildup rate (Section 3). We assume that on the long‐term average, (a) earthquakes and aseismic slip release moment at this same total rate (i.e., the moment budget is balanced) and (b) earthquakes (excluding aftershocks) obey a power‐law magnitude‐frequency distribution (MFD) up to a maximum magnitude (e.g., Michel et al. 2018). We build a suite of long‐term earthquake rate models based on these criteria and, for each model, evaluate its probability in light of the seismic catalog. We then multiply these probabilities by two physics‐ and data‐driven priors (Section 4). First, we compute the probability that the earthquakes in a given model would fit on the kinematically coupled part of the MHT, using the empirical relation between seismic moment and fault area (e.g., Kanamori & Brodsky, 2004; Ye et al., 2016). Second, we include the probability that the earthquakes in a given model would have to cross an aseismic barrier (in order to attain their full extent) and would be able to do so. Such statistics have been evaluated in quasi‐dynamic models of fault slip (Kaneko et al. 2010) and can be directly used in a probabilistic sense.

We define the probability of each seismicity model as , with the probability of the models total moment rate (from the moment buildup rate PDF), the probability of the model considering the seismic catalog, the probability of the model considering moment‐area scaling laws and the probability of the model considering the effect of aseismic barriers (Supporting Information Text S4; Table S1). Note that any constraint can be added to this approach as long as it can be formulated as a probability density function.

Constraints From Coupling and Seismicity

We use the coupling model of Dal Zilio et al. (2020) (Figure 1) to build a PDF of moment buildup rate on the MHT, estimate the available fault area, and estimate the size of potential aseismic barriers, effectively propagating uncertainties in the original geodetic data to all of these constraints. Summed over the region, moment deficit rate is normally distributed with a mean, , and standard deviation, , of 9.88 0.29 1019 N.m.yr−1 (Figure 1; Supporting Information Text S1).

Each long‐term earthquake rate model is characterized by a MFD, which we assume characterizes earthquakes excluding aftershocks. These MFDs depend on four parameters:, which controls the relative rates of small and large events; , the maximum magnitude; , the ratio of moment released by the earthquakes in the MFD to that released in total including aftershocks and aseismic deformation; and the total moment rate of all earthquakes in each model. Assuming that the moment budget is balanced and the MFD has a power‐law form up to a maximum magnitude, the rates of earthquakes at all magnitudes can be expressed in closed form and tied to the geodetic moment buildup rate (Avouac, 2015; Molnar, 1979). In light of the debate on the shape of MFDs for large magnitude earthquakes (Frohlich, 2007; Kagan, 1993; Wesnousky, 1995), we test two families of seismicity models. The first, which we call a truncated model, is a cumulative power‐law MFD (the Gutenberg‐Richter law) truncated at (e.g., Avouac, 2015; Michel et al., 2018; Stevens & Avouac, 2016). The second, which we call a tapered model, is formed from truncating an incremental power‐law MFD at , and has a tapered form in cumulative magnitude‐log‐frequency space (Rollins and Avouac, 2019, Supporting Information Text S2, Figure S2). Note that it is not the same as the tapered model used by (Kagan & Jackson, 2000) and subsequent studies. To compute the long‐term MFDs, we sample the space of and , assuming that and have uniform priors and that is normally distributed. The total moment release rate of each long‐term MFD is a sample from the PDF of the moment buildup rate (multiplied by the MFDs ), so that the starting distribution of MFDs intrinsically incorporates the distribution of .

For each MFD, at each magnitude , we compute the probability that a Poisson process with the event rate given by the MFD would, within the timespan of the seismic catalog, generate the exact number of earthquakes of magnitude that are in the seismic catalog. The catalog‐based probability for each model is then the product of these probabilities over the magnitude bins. In this way, we probabilistically estimate the long‐term earthquake rates at each magnitude based on the moment budget and seismic catalog (Michel et al., 2018; Rollins & Avouac, 2019).

To estimate , we need a declustered seismicity catalog (Supporting Information Text S3). We combine two earthquake catalogs (Supporting Information Text S6 and Figure S9), one covering the instrumental period (1995–2020, ANSS catalog, https://earthquake.usgs.gov/earthquakes/search/), the other including historical earthquakes (Bilham, 2019; Table S2), in a MFD with magnitude binning 0.5, corresponding roughly to the uncertainty on historical earthquake magnitudes. We decluster the instrumental catalog using the probabilistic approach of Marsan & Helmstetter (2017) and retain events that are located within 200 km northward from the main frontal thrust, the surface expression of the MHT (Figure 1, Supporting Information Text S6). We assume the historical catalog is declustered by nature and complete after 1500 (Bilham, 2019) above . As the 1897 Shillong earthquake was not on the MHT, we do not include this event.

Constraints From Earthquake Physics

We separately evaluate the product , the probability that an earthquake follows the moment‐area scaling and behaves according to the physics of fault slip governed by rate‐and‐state equations. is obtained by considering a moment‐area scaling of (Figure S4, Supporting Information Text S5 and S7). We assume the coupling contour 0.3 defines the limit between the seismogenic and creeping areas of the fault. is the probability of an event to pass through an aseismic barrier, which depends on the non‐dimensional parameter characterizing barrier efficiency, (Figure S7; Kaneko et al. 2010). This efficiency is defined as the ratio between the barrier resistance to a perturbation and the effective size of the perturbation, where is the barrier normal stress, and the constitutive parameters of the rate‐and‐state law, the slip rate as the seismic event propagates through the barrier, the interseismic slip rate, and the width of the barrier (Figure S6). and are the earthquakes mean stress drop and length, respectively. The subscripts VS and VW stand for Velocity Strengthening and Weakening, respectively describing the frictional behavior of the barrier and the seismic asperity in terms of rate‐and‐state friction (Dieterich, 1978; Ruina, 1983). is the proportion of stress redistributed by the earthquake that perturbs the barrier and depends on the geometry of the system, assuming stress is redistributed uniformly by the event (Kaneko et al. 2010). For the Himalaya case, all parameters determining , and their uncertainties, can either be estimated from the coupling models (Figures 1b and 1c) or from other published studies. We chose a priori PDFs for the parameters involved within reasonable, but very large, bounds (Supporting Information Text S5). From this, we sample the space of magnitudes, compute the corresponding spatial extent of ruptures and evaluate their probability considering the scaling and their probability to cross barriers. For instance, an earthquake of magnitude 9 would have to potentially cross 2 or 3 barriers and the probability of crossing these barriers can be stochastically tested. All magnitudes and locations can be tested and lumped into the PDF of earthquakes for this region. The complete description of the test for a given earthquake is given in Supporting Information Text S7.

We combine the PDF of seismic models that satisfy the moment budget and the MFD of the current and historical catalogs with that of the physical consistency of earthquakes, and evaluate the marginal PDF of the maximum magnitude earthquake and of its recurrence time. Finally, from , thus knowing the earthquake rate of each magnitude for each seismicity model and their related probability to be plausible, we calculate the probability to have an earthquake over a certain magnitude for a given period of time, (Michel et al., 2018; Rollins & Avouac, 2019). Such estimates are important inputs for seismic hazard mapping.

The Seismogenic Potential of the MHT

We first evaluate , the probabilistic long‐term earthquake model given just the PDF of moment buildup rate and the seismic catalog, for both the truncated and tapered forms (Figure 2a). For the truncated model, we define as the recurrence interval of an earthquake of magnitude . The marginal PDF of for the truncated model in the ‐ space is shown by the gray shaded distribution in Figure 2a. The marginal PDF of indicates plausible values between 8.15 and (Figure 2a), with two modes (i.e., peaks of the PDF) at  = 8.65 and 8.95. Those are associated with modes of at  = 400 and 1000 years, respectively. The probability of having an event of magnitude at a certain rate, , taking all events of all seismicity models into account, can also be estimated (Supporting Information Text S4). For  = 8.65, has two peaks (Figure 2a), one at  = 400 years arising from the actual recurrence time of (same recurrence time as for ) and one at 4000 years arising from any events with present in seismicity models with larger than . The marginal probability of (solid blue line in Figure 2a inset) peaks at 0.92.

Figure 2

(a) Seismic potential analysis using only the moment budget and the seismicity catalog as constraints. The blue and green elements in the figure are associated to the truncated and tapered magnitude‐frequency distribution models, respectively. The dots indicate the mean of the marginal probability density function (PDF) of the seismicity models. The dashed lines indicate the spread of the 1% best seismicity models. Red and orange crosses and error bars indicate the rate of occurrence of historical earthquakes and earthquakes occurring over the instrumental period, respectively. The gray shaded distribution indicates the marginal probability of the maximum magnitude earthquake and its rate. The solid lines on the axis indicate the marginal probability of , . The blue solid line on the earthquake frequency axis indicate the probability of frequency, assuming that is equal to the mode (i.e., peak of the PDF) of , , for the truncated model: . The dashed lines on the earthquake frequency axis indicate the probability of having an event of magnitude at a certain rate, : . considers all magnitudes in the seismicity models, while considers only the recurrence rate of . The top‐right inset shows the marginal probability of the b‐value. Note that the seismicity MFDs in the figure are not cumulative. (b) PDF (expressed in counts) of earthquakes considering the moment‐area scaling law and the Main Himalayan Thrust segmentation, from 408,500 events uniformly sampled between and 10. The pink histogram shows the PDF of earthquakes considering only the moment‐area scaling law while the yellow histogram considers both the scaling law and the coupling segmentation into account. The black line is the ratio between the yellow and pink histograms, representative the effect of the segmentation, alone.

The marginal PDF for the tapered model (Figure 2a) rises at 8.25, peaks at  = 8.65 and approaches zero again at 9.65. It is associated with a peak of at  = 1250 years. The marginal probability peaks at  = 0.75. As with the truncated form, earthquakes of magnitude larger than 10 cannot be ruled out when scaling laws and segmentation of coupling are not accounted for.

The distribution of plausible events considering the moment‐area scaling law is shown as pink histograms in Figures 2b and S12. Events with magnitude larger than are less probable because of the scaling law, whereas events with magnitude over are impossible (Figure 2b) because they simply cannot fit on the coupled part of the MHT.

The distribution of plausible events considering the combined effect of the moment‐area scaling law and the MHT segmentation is shown as yellow histograms in Figure 2b and S12. The ratio between the yellow and the pink histogram allows us to visualize the effect of the MHT segmentation by itself (black line in Figure 2b and S12). From this ratio, the influence of segmentation becomes obvious for events with , as these events must propagate through one or multiple barriers. Multiple parameters highlight the role of barriers, including event magnitude, moment‐area scaling and stress drop. For instance, we note that events with necessarily need to pass through all 3 barriers. However, smaller events with low moment‐area scaling (i.e., less than 1.5) have an anomalous length with respect to the mean scaling, which involves these events must break through multiple barriers, hence a decreased probability. Stress drop, , has also a strong influence, with events with MPa being less likely to power through barriers than events with MPa, hence there lower final probability. Other parameters do not seem to have a striking influence although they theoretically should, either due to the range of values tested or because their influence is hidden by that of other parameters. Note that all parameter distributions are controlled by the initial range of tested.

The combined effect of moment budget, seismicity catalog, scaling law and fault segmentation on the marginal PDF of the maximum magnitude earthquake, , and of the corresponding recurrence time is shown in Figure 3 and S13. For the truncated seismicity model, the scaling law and segmentation effect does not change the mode of , and hence of (solid dark blue line in Figure 3), but the marginal PDFs shrink significantly with no probable earthquake above The corresponding recurrence time associated with reduces significantly from 3981 years to 2512 years. The second peak at 400 years does not move as it is arises from the recurrence rate of the seismicity models (for ) which hasnt changed from the scaling law and segmentation effects. For the tapered seismicity model, the scaling law and segmentation effect changes the mode of slightly from to . were already initially improbable and the marginal PDF of the maximum magnitude earthquake shrinks only slightly. Consequently, the corresponding recurrence time associated with just slightly reduces from 1259 to 1000 years.

Figure 3

(a), and (c) show the evolution of the marginal probability density function (PDF) of the maximum magnitude for an earthquake (a) and of the marginal PDF of the corresponding recurrence times (c) when progressively adding constraints, for the case of truncated MFDs. Dotted lines indicate the marginal PDFs for the moment budget analysis, the dashed lines indicate marginal PDFs when adding the earthquake scaling constraint and the continuous line indicate the marginal PDFs with all constraints. (e) Shows the probability of occurrence of earthquakes of magnitude larger than Mw over a period of X yrs for the case of truncated MFDs. Dashed lines indicate the PDF considering only the moment budget analysis while continuous lines show the PDF including all observational and physics‐based constraints. We show the probability of occurrence of such events for 4 time periods, 30, 100, 1e3, and 1e4 yrs. (b), (d), and (f) are the same as (a), (c), and (e) for the case of tapered MFDs. (g) Marginal probability of rupture extent for a nine earthquake. This PDF directly derives from our PDF of seismic models accounting for all observational and physics‐based constraints. K and D stand for Kathmandu and Delhi. (d) Cumulative PDF of a source distance from Kathmandu and Delhi.

The probability to have at least an earthquake with magnitude over 9 for a time period of 1000 years, , is equal to ∼12% and 3% for the truncated and tapered seismicity models, respectively (Figure 3e and 3f). The effect of scaling law and segmentation only arises for large time periods, on the order of the recurrence time of events with (∼1–4 103 yrs). Past this limit, the different physics‐based constraints play a more important role as larger event are expected to happen during this time period. For example, for a 10,000 year period using the truncated model, are less probable using the scaling law and segmentation constraints than if it was not taken into account. Whereas for shorter time periods (e.g., 100 years), events will anyways be improbable. Finally, we can evaluate the probability of an event of a given magnitude for a given position along the megathrust considering the segmentation of the megathrust and the moment‐area scaling law. Figure 3g shows an example of the probability of a point of the MHT to be part of a . Consequently, the distance between a chosen site and the seismic source can be calculated (Figure 3h), such estimate being one important element of the seismic hazard evaluation.

Discussion and Conclusion

Our estimate of the seismic potential of the Himalaya accounts for observations and for our current knowledge of the physics of earthquakes. In addition, we consider uncertainties on the MHT coupling distribution (Dal Zilio et al., 2020), based on uncertainties in geodetic measurements and rheology, and the uncertainties on the MFD, based on the uncertainties on the magnitudes of recorded and historical earthquakes. We consider uncertainties on the frictional constitutive parameters of the MHT (e.g., barriers size, position and their efficiency) as well as uncertainties on the global statistics of earthquakes (i.e., moment‐area scaling law). The combination of all these uncertainties suggest earthquakes of magnitude comparable to that of the 1950 earthquake (∼8.7, Bilham, 2019) are the most probable candidates for being the largest possible earthquakes in the Himalayas. The associated most probable recurrence time for such event exceeds 1000 years. However, given the model of magnitude frequency distribution considered, the probability of occurrence of an earthquake of magnitude larger than 8 over a 100 years ranges from ∼60% to 80%, a number high enough for large hence damaging earthquakes to be considered seriously (Figure 3).

Our analysis depends on several assumptions, including assumptions impacting the moment deficit rate and the seismic catalogs. First, the historical catalog includes the 1950 Assam earthquake although this rupture only partially overlaps with the coupling map. We obtain similar results with or without including this event in the catalog (Figure S13 and S14). Second, the width of the seismogenic zone and thus the associated barriers size are based on a coupling threshold of 0.3. Tests with coupling thresholds of 0.2 and 0.4 (Figure S12) shifts the PDF of maximum magnitude of only ∼0.1. Third, the current coupling model may not be representative of the long‐term coupling model (e.g., Hardebeck & Loveless, 2018). Fourth, we assume a Poissonian distribution of events. However, whether an event passes through barriers or not must influence the magnitude and the timing of the next large earthquake. More fundamentally, we do not account for time dependent processes while the timing of previous large earthquakes must be of influence. Fifth, and this is impossible to test yet, we assume the current or the historical earthquake catalog is representative of the long term seismicity rate. Large earthquakes ( might not follow the same MFD as small ones. In such case, neither of the seismicity models tested (i.e., truncated and tapered) would be representative of the seismicity for the MHT.

In conclusion, our study shows how the seismic potential of a fault can be evaluated taking into account first‐order physics of fault slip and all sources of uncertainties. We introduce constraints by considering global scaling laws and the frictional segmentation of the MHT. We show that in the case of the MHT, are implausible considering the moment‐area scaling law, and that get less and less probable due to the segmentation. Such estimation of the seismic potential can then be used to construct probabilistic seismic hazard assessments.

Supporting Information S1

Click here for additional data file.