Geodetic plate coupling and seismic potential on the main Himalayan thrust in Nepal

Abstract

We model interseismic plate coupling distribution on the Main Himalayan Thrust (MHT) in Nepal through the inversion of secular GNSS velocity data. To complement previously published data, we compile velocity data from ten additional continuous stations to improve spatial resolution in central and mid-western Nepal. A regional non-planar structural model is adopted to reproduce the MHT fault plane. In general, the coupling pattern seems nearly binary, indicating that transition from full coupling to decoupling is occurring sharply in very narrow zones. In eastern Nepal (> 84.5ºE), plate coupling is very strong from the surface to intermediate depths, including the source region of the 2015 Mw 7.8 Gorkha earthquake. Geodetically estimated slip deficit rates are consistent with the rupture history of great earthquakes in the east revealed by geomorphological observations. In the west (< 83.0ºE), a weakly coupled zone extends laterally at intermediate depths, whereas the coupling in the shallower part remains very strong. Although the slip deficit rate in the west is significantly smaller than that in the east, seismic moment accumulated, since the last complete rupture in 1505 may be capable of generating a future great event. In central Nepal, estimated slip deficit rates are comparable with those in the east, and no great event has been documented over the past several centuries. Thus, the seismic risk may be most urgent in central Nepal. The development of local seismicity and crustal deformation should be carefully monitored.

Graphical Abstract

Introduction

The Mw 7.8 Gorkha-Nepal earthquake on April 25, 2015 has reminded us again of the high seismic hazard in the Himalayan countries, in the convergence zone between the Indian and Eurasian plates. Interseismic tectonic shortening across the Himalayas has been released by slip along a major basal thrust fault, the Main Himalayan Thrust (MHT), which divides India and Eurasia (Bilham et al. 2001; Avouac 2007). Major destructive earthquakes have been documented in the past several centuries along the MHT (Ambraseys and Douglas 2004), and the 2015 Gorkha earthquake was the largest event in Nepal since the 1934 Mw 8.2 Bihar earthquake (Sapkota et al. 2013) (Fig. 1). In central and western Nepal, just west of the 2015 rupture zone, is the 400-km-long “Central Seismic Gap”, where no destructive earthquake has occurred over the last several centuries since a major seismic event in 1505. Recent surface displacement data from the Global Navigation Satellite System (GNSS) network across the Himalayan arc indicate that about half of the India–Eurasia convergence has been accommodated in this arc (Ader et al. 2012; Stevens and Avouac 2015). Plate coupling distribution and its temporal evolution are of great concern for understanding elastic strain build-up on the plate boundary zone and estimating ultimate strain release by future earthquakes.

Fig. 1

Great earthquakes in Nepal since 1500 AD. Gray-shaded boxes show approximate extents of source region of the historical events (Ambraseys and Douglas 2004; Bollinger et al. 2016). Yellow arrows show Indian plate motion relative to the Eurasian plate predicted from the ITRF2014 plate motion model (Altamimi et al. 2017). Red line delineates surface trace of the Main Himalayan Thrust, the plate boundary between India and Eurasia (Hubbard et al. 2016). Circles show locations of continuous GNSS stations used in this study; purple: NERDiM stations constructed in 2016–2018 and white: pre-existing stations whose data are available through GAGE or IGS

Interferometric Synthetic Aperture Radar (InSAR) measurements delineated detailed surface deformation patterns caused by the 2015 Gorkha earthquake (Lindsey et al. 2015; Kobayashi et al. 2015; Wang and Fialco 2018). This may be the best example showing that InSAR is particularly effective in monitoring this region where continuous GNSS network deployment is limited. However, InSAR measurement can collect only line-of-sight (LOS) displacements, and failed to detect the most dominant coseismic signal at the time of the 2015 Gorkha earthquake; that is, horizontal elastic rebounds exceeding 1.8 m nearly parallel to the descending path of the SAR satellite (Zhang et al. 2015). Thus, densification of the GNSS network, if it can be realized soon, may contribute to provide a more complete dataset of surface deformation in conjunction with InSAR measurements.

Soon after the 2015 event, the Nepal–Japan joint project “Integrated Research on Great Earthquakes and Disaster Mitigation in Nepal Himalaya” (hereinafter called NERDiM: Nepal earthquake research and disaster mitigation) was started (Japan Science and Technology Agency 2023). Under this project, the authors organized a geodetic group and constructed continuous GNSS stations in central and mid-western Nepal to enhance the spatial resolution of the pre-existing domestic network. The main purpose was to estimate long-term earthquake generation potential in the Himalayan Frontal Zone using the space–time variation of crustal deformation. In this paper, we present a plate coupling distribution in Nepal using GNSS data and a more realistic plate boundary fault model. We then discuss seismic potential along the MHT in Nepal from a geodetic point of view.

Data

Installation of continuous Global Positioning System (GPS) stations was initiated in the late 1990s in Nepal, but station coverage remained sparse until 2004. Soon after the 2015 Gorkha earthquake, many GNSS stations were renewed or added, especially in and around the source region, including Kathmandu Valley. The total number of the continuous stations in Nepal exceeds 50, but the number of stations where stable observations have been maintained is much lower. We collected raw RINEX (Receiver Independent Exchange Format) data obtained at these stations since 2004 through the Data Center of GAGE (the Geodetic Facility for the Advancement of Geoscience, formerly UNAVCO: University NAVSTAR Consortium).

Additionally, the geodetic group of NERDiM installed ten GNSS stations in 2016–2018 to reinforce the pre-existing domestic network and particularly improve spatial resolution near the plate boundary. Altitudes of the NERDiM stations vary from 10 m near the Indian border to 2500 m in the highland. The distribution of GNSS stations is shown in Fig. 1.

We process all RINEX data up to 2023.0 based on the International Terrestrial Reference Frame (ITRF) 2014 (Altamimi et al. 2016) along with data from International GNSS Service (IGS) reference stations in South Asia, using the IGS precise orbits and the Bernese GNSS Software v5.2 (Dach et al. 2015). Assuming a linear trend and annual and semi-annual seasonal variations, secular site velocities (displacement rates) are determined from time series of daily positions. However, the periods and lengths covered by the time series vary widely from site to site. In addition, the occurrence of the 2015 Gorkha earthquake may have severely disturbed the local deformation field in and around the source region. For example, coordinate time series of eight selected sites are shown in Additional file 1: Fig. S1. For sites with enough observations both before and after the earthquake, we calculate velocities separately from the pre-earthquake data, from the post-earthquake data, and from all data covering the whole observation period after correcting the coseismic jumps. When the post-earthquake rate differs significantly from the pre-earthquake rate, likely because of coseismic stress change and postseismic transients, the pre-earthquake rate is adopted as the secular, more stationary velocity. For sites with insufficient or unstable observations before the earthquake, we carefully check the location and omit those sites from further analysis if they are included in the postseismic deformation area (Tian et al. 2020). Sites where the observation period is shorter than 2 years are not taken into consideration. Seven sites are omitted under these criteria.

Uncertainties in the GNSS time series caused by unmodeled sources other than the linear trend, seasonal variations, and coseismic jumps are difficult to quantify. We follow the empirical formula by Ader et al. (2012) to evaluate these uncertainties. Typical velocity uncertainties are 0.60, 0.57, and 1.87 mm/yr for the east, north, and up components, respectively. Estimates of the secular velocities and corresponding uncertainties are summarized in Additional file 1: Table S1.

We compile additional velocity data from previous studies to cover regions outside of Nepal. However, in addition to the different observation periods, the version of ITRF that the GNSS data processing is based on has not been the same throughout the studies; versions used have included ITRF2005 (Ader et al. 2012), ITRF2008 (Wang and Shen 2020), and ITRF2014 (this study). First, we check apparent velocity changes produced by different versions of ITRF using transformation parameters published by the International Earth Rotation and Reference System Service (2023). The apparent velocity changes are smaller than 0.3 mm/yr in three components in central Nepal. Therefore, we no longer take into consideration the different ITRF versions. Second, we check the velocity values at the sites overlapping between this study and that of Ader et al. (2012). The root-mean-squares of the velocity discrepancies at 20 overlapping sites are 1.2 mm/yr and 1.7 mm/yr in the east and north components, respectively. These are slightly larger than 1-σ uncertainties of the velocity components (Additional file 1: Table S1) and may be attributed to different observation periods and lengths between the studies and the different data processing strategies employed. We do not combine the three sets of velocity data (Ader et al. 2012; Wang and Shen 2020; this study) into one before the inversion analysis. Adjustments of velocity data are done during the inversion to better fit the reference frame.

Finally, we transform the reference frame of the velocity field into the Indian plate-fixed frame using the ITRF2014 plate motion model (Altamimi, et al. 2017). This is done, because the Indian subcontinent can be treated as a nearly rigid plate with negligible internal deformation, whereas Nepal and the Tibetan Plateau are included in an active intraplate deformation zone (e.g., Kreemer et al. 2014). Therefore, the Indian plate is more suitable as the reference in investigating Indo-Eurasia plate convergence. The horizontal velocity field is shown in Fig. 2.

Fig. 2

Horizontal velocities in the Indian plate-fixed reference frame. Green, red, and blue arrows represent velocities derived in this study, those from Ader et al. (2012), and those from Wang and Shen (2020), respectively. Transformation of the reference frame is performed using the ITRF2014 plate motion model (Altamimi et al. 2017). Yellow allows show velocities predicted from the ITRF2014 plate motion model when assuming that the Tibetan Plateau belongs to the main part of the Eurasian plate

In addition, we further derive the strain rate field from the velocity data by spatial interpolation and filtering (Shen et al. 1996), with the distance-decaying constant taken as 50 km (Fig. 3). Because the strain field is free from the block motion, it is more suitable than the displacement field as a direct overview of the internal deformation.

Fig. 3

Horizontal strain rate field. a Areal dilatation, b maximum shear, and c principal strain. Observed velocities in Fig. 2 are differentiated by spatial interpolation and filtering with a distance-decaying constant of 50 km (Shen et al. 1996)

Methods

Previous studies of Indo-Eurasia plate convergence adopted simple geometries for the plate boundary, consisting of laterally connected planar sub-fault segments (e.g., Ader et al. 2012; Stevens and Avouac 2015; Dal Zilio et al. 2020) where the dip angle was fixed to 10° throughout. However, based on numerical experiments, Dal Zilio et al. (2019) demonstrated that the seismic cycle in the Himalayan region is controlled by the frictional properties and non-planar geometry of the MHT. The geometry of the plate interface should be carefully considered when plate coupling strength in the interseismic stage is estimated.

In this study, we adopt a regional non-planar structural model to reproduce the MHT fault plane in Nepal (Hubbard et al. 2016). This model, hereafter called the MHT model, was derived from geological investigation and is consistent with seismic observations. It is characterized by down-dip variation of the dip-angle, with a steep ramp at intermediate depths separating upper and lower décollements (Additional file 1: Fig. S2a). A depth contour of 30 km in this model seems to well define a lower limit of regional seismicity (Hubbard et al. 2016). In parallel, we test another nearly planar fault plane to check how the fault geometry affects plate coupling estimates. This second model is based on a global subduction zone geometry model (Hayes 2018). The fault plane has a gentle slope of about 10–11° from the surface to a depth of 30 km (Slab-2 model hereafter, Additional file 1: Fig. S2b). In both models, we set nodes with a lateral spacing of 30–40 km along depth contours every 5 km.

The horizontal GNSS velocities in Fig. 2 are inverted to simultaneously estimate the slip deficit rate distribution on the plate boundary and rigid rotation of the overriding block using TDEFNODE (McCaffrey 2009, 2023). Incorporating the latter is necessary, because the horizontal motion in the Nepal–Tibet region is notably different from that predicted from the global plate motion model (Fig. 2). The software can solve both the long-term steady motion of the sites and short-term transients, but we examine only the steady motion in this study. The surface response to a unit slip at each node on the fault is calculated using an elastic half-space dislocation model (Okada 1992). The inversion is performed in TDEFNODE through a combination of simulated annealing and grid search. This process is iterated several times.

Decoupling, or fault unlocking, is imposed on the lower margin of the model fault at a depth of 30 km, taking into account the correlation between seismicity, temperature, and interseismic coupling in central Nepal. According to Ader et al. (2012), seismicity has a peak at a mid-crustal depth, where the temperature is 300–350ºC, and the coupling then starts to drop abruptly; coupling and seismicity are expected to fade out at a depth of around 30 km with a higher temperature. On the other hand, we do not set any boundary conditions at the surface or at the eastern and western ends. Upper limit and direction of the slip deficit on the plate boundary is implicitly constrained by rigid rotation of the overriding block. In TDEFNODE, the plate coupling ratio at the interseismic stage, defined as the ratio of the estimated slip deficit rate divided by the block velocity on each fault patch, is regulated to fall between 0 and 1. This ensures that the slip deficit is neither negative nor faster than the relative block motion (McCaffrey 2009).

The strain rate field suggests that a single block assumption for the overriding block may not be adequate to explain the areal dilatation and east–west extension observed in the Tibetan Plateau (Figs. 3a and c). However, we do not divide it into multiple segments. This is chiefly because the major concern in this study is regional deformation within Nepal and its surrounding area, and the velocity field there seems consistent with a single block treatment. For the Tibetan Plateau, Wang and Shen (2020) reported that the deformation field was predominantly continuous, and that the deformation gradient only existed perpendicular to the Indo-Eurasia relative plate motion. These results do not support a simple segmentation of the overriding block into east and west components. Therefore, we keep a single block treatment for the overriding block, but avoid incorporating far-field GNSS velocity data more widely distributed in the active Tibetan Plateau.

Prior to the inversion of the real data, we carry out a spatial resolution test, also called a checkerboard test, to confirm the performance of the applied method. The simulated pattern of the slip deficit rate for the modeled MHT fault is well reproduced through the inversion (Additional file 1: Fig. S3).

Results

The distribution of the interseismic slip deficit rate and plate coupling ratio along the modeled MHT fault, along with the observed and calculated site velocities, are shown in Fig. 4, where the depth of the fault’s lower margin is set to 30 km. The block velocity is predicted from the estimated rotation pole and angular velocity of the overriding block relative to India. They are estimated, respectively, as (53.85 ± 6.84ºE, 28.17 ± 7.96ºN) and −0.370 ± 0.031º/Ma in this inversion, predicting that the velocities of the overriding block relative to India would change from 16.0 mm/yr in rate and 185º in azimuth at 80ºE to 21.2 mm/yr and 191º at 89ºE along the plate boundary. Because the predictions from the global ITRF2014 plate motion model (Altamimi et al. 2017) at the same locations are 35.4 mm/yr and 192º and 39.2 mm/yr and 196º, respectively, 45%–54% of the relative motion of Indo-Eurasia is taken up as shortening in a direction deflected by 5–7º counter-clockwise in the plate boundary zone in Nepal and its surroundings. Block velocities at GNSS sites are shown in Additional file 1: Fig. S4, along with elastic velocities produced by plate coupling.

Fig. 4

Interseismic plate coupling in Nepal. The MHT model (Additional file 1: Fig. S2a) is used for the plate interface. a Distribution of plate coupling ratio along with the observed (black arrows) and model-predicted (yellow arrows) site velocities relative to the Indian plate. b Distribution of the slip deficit rate vector on the plate interface

Figures 4 shows that nearly full coupling is realized in a shallower (< 10 km) part of the MHT fault, and that the slip deficit rate itself becomes larger toward the east because of the rigid rotation of the overriding block. In contrast, weakly coupled patches are identified in places at intermediate depths (10–20 km). In western Nepal, the decoupled zone extends from the intermediate to lowermost depths. Correlation between the plate coupling strength and fault dip-angle is not clear. Looking more closely at the source region of the 2015 Gorkha earthquake at about 84.5–86.5ºE and 10–20 km in depth, nearly full coupling is realized there. Because we made efforts to exclude the effect of the earthquake from the dataset, the results likely reflect robust strain accumulation prior to the earthquake. Very weak coupling at greater depths is natural because of the constraint of unlocking at the lower margin of the fault. In general, the coupling pattern seems nearly binary (1 or 0), and transition from full coupling to decoupling is occurring sharply in very narrow zones.

Calculated site velocities well reproduce the observed data in central and eastern Nepal (Fig. 4a), where most residuals are within a few mm/yr in root-mean-square. Conversely, the fitting becomes worse toward the outer margin of the study area. As expected based on the previous section, the assumption of a single overriding block is inadequate to reproduce the observed east–west extension in the Tibetan Plateau. This will be further investigated in a subsequent study. Because the upper limit of the slip deficit rate is constrained in the inversion by rigid rotation of the overriding block, there may be some correlation between these estimates. To check this empirically, we repeat the inversions for slip deficit rates, fixing the rigid rotation to several pre-defined values. Based on the inversion results described above, we fix the pole position of the rigid rotation but vary its angular velocity by ± 5%, ± 10%, ± 15%, ± 20%, and ± 30% from the original estimate. We use the root-mean-squares of the residuals between observed and calculated site velocities to evaluate the robustness of the slip deficit rate estimates. The results are shown in Additional file 1: Fig. S5. Because the relative block motion in this region is directed nearly north–south, the root-mean-square of the north component of velocity is indicative, whereas that of the east component is less informative. Although the minimum root-mean-square is obtained when the angular velocity increment is −5%, very similar values are also seen at 0% and −10%. The spatial pattern of the plate coupling ratio on the plate boundary does not change significantly. Summarizing the results, the slip deficit rate estimate in this inversion may contain uncertainty of several to 10% of its amplitude.

The strain rate fields derived from the calculated site velocities are shown in Additional file 1: Fig. S6. Comparing Figs. 3 and S6 reveals that the strain rate fields predicted from the model well reproduce the deformation observed in central and eastern Nepal, whereas the fitting is slightly worse in the west. In general, crustal shortening of about 10^–7 /yr is identified in a narrow plate boundary zone with a width of about 200 km. Interseismic strain accumulation decays sharply beyond the Himalayas.

Discussion

Fault geometry and plate coupling

How the geometry of the fault plane affects plate coupling estimates is the crucial question in this study. In addition to the non-planar, structural MHT model, we test the nearly planar Slab-2 model, keeping the depth of the lower margin of the fault to 30 km, and the constraint of decoupling at the same depth. In the results of the Slab-2 model (Fig. 5), lateral variation of the plate coupling at intermediate depths (10–20 km) seems more emphasized. The weakly coupled zone in western Nepal is wider than in the MHT model, and more small patches of weak coupling are identified in central Nepal. Furthermore, weak coupling is estimated in places near the surface, in contrast to the nearly full coupling throughout in the MHT model. Because the residuals between observed and calculated site velocities are slightly larger than in the MHT model along the plate boundary, and because surface creep motion has not been detected in Nepal (e.g., Bilham et al. 2017), the weak coupling near the surface estimated in the Slab-2 model may be less reliable.

Fig. 5

Interseismic plate coupling for the Slab-2 model (Additional file 1: Fig. S2b). a Distribution of plate coupling ratio along with the observed (black arrows) and model-predicted (yellow arrows) site velocities relative to the Indian plate. b Distribution of the slip deficit rate vector on the plate interface

We repeat the same inversion analysis for the MHT model with a minor modification of the depth of the lower margin of the fault to 35 km. This modification results in little change of the plate coupling distribution from the original model except near the lowermost depth (Additional file 1: Fig. S7). Based on the results shown in Figs. 4, 5 and Additional file 1: Fig. S7, we suggest that two characteristics of the plate coupling pattern in Nepal, the eastward increase of the slip deficit rate and the existence of a wide weakly coupled zone in western Nepal, are robust regardless the fault geometry, and these characteristics may be mainly driven by data and modeling.

Recently, Dal Zilio et al. (2020) estimated interseismic geodetic coupling along the MHT between roughly 73ºE and 96ºE, based on four planar sub-fault segments and a probabilistic Bayesian approach. They detected lateral variability of the coupling characterized by large (300–450 km) and discrete highly coupled segments separated by a smaller (~ 100 km) weak zone. Conversely, previous studies using similar planar sub-fault segments (Ader et al. 2012; Stevens and Avouac 2015) estimated a nearly homogenous distribution of plate coupling along the plate boundary with no significant lateral variation. Dal Zilio et al. (2020) concluded that the homogenous coupling in the previous studies were mainly attributed to a strong constraint on the slip distribution along the fault.

According to Dal Zilio et al. (2020), there are two distinctive zones at around 83ºE and 87ºE in the area of this study (79–91ºE), where very weak coupling was estimated probabilistically from near the surface to the deeper part. For the former zone (around 83ºE), we installed four NERDiM GNSS stations, two each on both sides of the plate boundary (Fig. 1), to overcome the observation gap. As a result, we estimate that the coupling in the shallower and intermediate depths is still strong (Fig. 4). Similarly, we found no weak coupling in the latter zone (87ºE) without additional GNSS stations. These discrepancies of plate coupling estimates may not be attributed to the differences of the fault geometry, because different results were derived from previous studies using similar planar sub-fault segments (Ader et al. 2012 and Stevens and Avouac 2015 vs. Dal Zilio et al. 2020), and we confirmed that the MHT and Slab-2 model yield similar plate coupling distributions. Plate coupling estimates are more likely to depend on the spatial distribution of data and the modeling approach.

We find systematic residuals between the observed and calculated site velocities in Figs. 4, 5, and Additional file 1: Fig. S7; that is, eastward near the eastern margin and westward in the western. These discrepancies may result from the single block treatment of the active Tibetan Plateau. In addition, they strongly suggest that a more realistic Earth model should be adopted in further modeling. To date, realistic topography of the Earth, including its sphericity, has not been considered in the modeling of plate coupling in the Himalayan region. Strictly speaking, however, such topography should be incorporated into modeling, because the study area includes the highest mountains in the world, and the spatial extent exceeds 1000 km. Langer et al. (2019) calculated coseismic surface displacements produced by the 2015 Gorkha earthquake for a realistic topographic Earth model with a height difference of 8.33 km. They found a maximum difference of about 10% from the displacements on a flat Earth model. On the other hand, Banerjee et al. (2005) discussed how Earth’s sphericity affects the horizontal displacement field produced by thrust faulting on a long fault. Divergence of the displacement amplitude between homogeneous spherical Earth and half-space increases quickly beyond about 10 geocentric degrees from the source, whereas agreement is found within a few percent out to about 5 geocentric degrees. Although it is difficult to evaluate this precisely, we should keep in mind that the slip deficit rates estimated using an elastic half-space dislocation model and ignoring realistic topography of the Earth may be partially biased. Hori et al. (2021) provided an answer to this problem, calculating Green’s functions based on a numerical model that incorporated realistic topography and subsurface structure of the subduction zones in Japan, along with the ellipsoidal shape of the Earth. Although performing such calculations is generally challenging, the approach of Hori et al. (2021) is also promising for modeling the Himalayan region.

Slip deficit rate and seismic potential

Figure 6 shows the current along-strike variation of the slip deficit rate and its vector azimuth estimated in this study and the rupture history of great earthquakes in Nepal derived from Bollinger et al. (2016). Here, the slip deficit rates at nodes are averaged across the strike from the surface to the depth of 25 km. The longitudes along the horizontal axis indicate the locations of surface nodes. In addition, the moving averages of the slip deficit rates with 3-nodes-widths are shown. We evaluate long-term seismic potential from a geodetic point of view, dividing the country into west (roughly < 83.0ºE), east (> 84.5ºE), and central regions (83.0–84.5ºE).

Fig. 6

Lateral variation of the slip deficit rate and its azimuth along the MHT. Slip deficit rates estimated at nodes on the MHT model fault are averaged across the strike from the surface to the depth of 25 km. Orange broken line shows the moving averages of the rate with 3-nodes-widths. Rupture history of the great earthquakes in Nepal is illustrated in the bottom of the figure (Bollinger et al. 2016), where red and blue lines indicate approximate lateral extent of the complete rupture propagated to the surface and the blind slip remained in the down-dip part, respectively

In the east, two types of rupture patterns have been evidenced: blind events (Mw ~ 7.8) remaining in the down-dip part of the seismogenic zone, and great events (Mw 8 +) propagating to the surface (e.g., Sapkota et al. 2013; Bollinger et al. 2016; Dal Zilio et al. 2019). The 1934 Bihar earthquake (Mw 8.2), the successor of the 1255 great event after 679 years, ruptured the surface at least 150 km in length between 85.8ºE and 87.3ºE (Sapkota et al. 2013). Assuming a source length of 150 km, a depth of the seismogenic zone of 25 km, and a rigidity of the crust of 30 GPa, and adopting 18.5 mm/yr as the long-term, depth-averaged slip deficit rate from Fig. 6, we estimate the total seismic moment accumulated in 679 years as 8.14 × 10²¹ Nm, which is equivalent to a Mw 8.54 event. This estimated magnitude may be the maximum value because the accumulated moment would be partially released by intermittent blind ruptures such as the 1833 event or taken up by potential transient movement caused by the inelastic property of the crust. Therefore, the geodetic slip deficit rate in the eastern part of Nepal is not far beyond the long-term slip rate expected based on geomorphological investigation.

On the other hand, no great earthquake has occurred in the western part of Nepal since the 1505 great earthquake. In the source region of this event, roughly between 80.0ºE and 83.0ºE, the depth-averaged slip deficit rate is about 7.3–16.3 mm/yr, significantly smaller than the rate in the eastern part (Fig. 6). This smaller slip deficit rate results in a lower rate of seismic moment accumulation, which would lead to a longer recurrence time to the next great event. Thus, the geodetic slip deficit rate along the MHT seems consistent with the long period of historical quiescence in western Nepal. We estimate that current seismic potential in the region of 80.0–83.0ºE, adopting 10.4 mm/yr as the depth-averaged slip deficit rate. Assuming a source length of 300 km, a depth of the seismogenic zone of 25 km, and a rigidity of the crust of 30 GPa, the total seismic moment since the 1505 event amounts to 6.98 × 10²¹ Nm. Thus, the source region of the 1505 event is presently capable of generating a Mw 8.50 event. However, it should be noted that the weakly coupled zone extending at intermediate depths in the west may contribute to the differences in rupture history compared with the east.

In the central part, roughly between 83.0ºE and 84.5ºE, paleseismological study clarified that a complete rupture to the surface occurred in 1344 or 1408, but no comparable event has been identified since then (Bollinger et al. 2016). Taking into consideration the longer quiescence time and the larger slip deficit rate than in the western part (< 83.0ºE), seismic risk in the central part may be the most urgent in Nepal. Full attention should be paid to the temporal progress of ongoing local seismicity and crustal deformation.

Conclusions

Based on secular GNSS velocities compiled from previous studies and our original observations, we estimate the interseismic plate coupling distribution on the MHT fault plane in Nepal, along with rigid rotation of the overriding block relative to India. We adopt a regional non-planar structural model to reproduce a fault plane characterized by a steep ramp at intermediate depths separating upper and lower décollements. The block velocity is estimated to increase toward the east along the MHT, but is about half of the relative motion of India and Eurasia predicted from the global plate model. In general, the results indicate crustal shortening of about 10^–7 /yr occurring in a narrow plate boundary zone with a width of about 200 km, decaying sharply beyond the Himalayas. The coupling pattern estimated for the MHT is nearly binary, indicating a sharp transition from full coupling to decoupling in very narrow zones. At shallower depths (< 10 km), the coupling is very strong from east to west along the MHT, whereas the slip deficit rate itself is shown to become larger toward the east because of the rigid rotation of the overriding block. Conversely, the coupling pattern shows lateral variation at intermediate depths (10–20 km). In the east, the coupling is very strong from the surface to intermediate depths, including the source region of the 2015 Mw 7.8 Gorkha earthquake. In the west, however, a weakly coupled zone extends laterally at intermediate depths. In general, correlation between plate coupling strength and fault dip-angle is not clear.

We evaluate seismic potential along the MHT in Nepal using geodetically estimated slip deficit rates and the rupture history of great earthquakes derived from geomorphological investigations. In the east, the geodetic slip deficit rate is large enough to compensate for seismic moment released by a sequence of great earthquakes over the past centuries. In the west, the seismic moment accumulated in the source region of the 1505 rupture event is estimated to be equivalent to a Mw 8.50 event, although the depth-dependent coupling pattern differed from in the east. Seismic risk may be most urgent in the central Nepal, taking into consideration the longer quiescence time of the earthquakes and the larger slip deficit rate than in the western part of the country.