Machine learning-based aftershock seismicity of the 2015 Gorkha earthquake controlled by flat-ramp geometry and a tear fault

1.   Introduction

Nepal is situated at the confluence of the Indian Plate and the Eurasian Plate, where continental collision has formed one of the most geologically active regions in the world, characterized by intense crustal deformation and earthquake hazards (Ding L et al., 2017; Bilham, 2019). Seismic activity in Nepal is predominantly governed by the Main Himalayan Thrust (MHT), a north-dipping low-angle megathrust that emerges from the convergence of multiple faults (Ni J and Barazangi, 1984; Pandey et al., 1995; Ader et al., 2012). Among the prominent faults contributing to this complex structure are the Main Frontal Thrust (MFT), Main Boundary Thrust (MBT), and Main Central Thrust (MCT) (Figure 1). The intricate geological framework was formed via millions of years of subduction and collision, accompanied by numerous seismic events (Dal Zilio et al., 2021).

Figure  1.  Map of station distribution. The inset shows the location of the study region. The two red stars mark the 2015 ${M}_{{\mathrm{W}}}$7.8 main shock (left) and ${M}_{{\mathrm{W}}}$7.3 aftershock. Red and blue triangles indicate selected seismic stations of the Himalayan-Xizang Continental Lithosphere during Mountain Building (Hi-CLIMB; Nábělek et al., 2009) array and the Nepal Array Measuring Aftershock Seismicity Trailing Earthquake (NAMASTE; Karplus et al., 2020), respectively. MFT: Main Frontal Thrust; MBT: Main Boundary Thrust; MCT: Main Central Thrust.

DownLoad: Full-Size Img PowerPoint

The 7.8 Gorkha earthquake, which occurred on April 25, 2015, serves as compelling evidence of the region's intense seismic activity, with a devastating toll of more than 8,700 deaths, 22,300 injuries, and 760,000 damaged buildings (Gautam et al., 2016). The mainshock rupture propagated eastwards along the MHT and triggered a significant ${M}_{{\mathrm{W}}}$7.3 aftershock at the eastern end of the rupture 17 days later (Avouac et al., 2015; Galetzka et al., 2015; Grandin et al., 2015; Lindsey et al., 2018; Bai L et al., 2016). Ultimately, the earthquake sequence affected an area measuring 150 × 80 ${{\mathrm{km}}}^{2}$ (Adhikari et al., 2015; Bai L et al., 2016), wherein aftershock activity was distributed along the tectonic strike and coincided with the co-seismic rupture zone of the mainshock (Avouac et al., 2015; Elliott et al., 2016; Hubbard et al., 2016).

Previous studies have suggested that the aftershocks were controlled by particular geological structures of the MHT (Hubbard et al., 2016; Baillard et al., 2017; Hoste-Colomer et al., 2017; Wang X et al., 2017; Mendoza et al., 2019; Yamada et al., 2020; Adhikari et al., 2023), and various interpretations have been proposed, such as a flat-ramp geometry (e.g., Hubbard et al., 2016), a duplex structure (e.g., Mendoza et al., 2019) along dip direction, or a tear fault (e.g., Hoste-Colomer et al., 2017) along strike. In the evaluation of the geological structures along dip direction, receiver function images (Hetényi and Bus, 2007; Nábělek et al., 2009; Caldwell et al., 2013) and seismic reflection images (Lemonnier et al., 1999; Berger et al., 2004; Robert et al., 2011) provide a general correlation with certain structures, but a comprehensive understanding and subsurface configuration of the MHT fault system is lacking. The proposed tear fault has been shown to influence the distribution of aftershock clusters along strike (Hoste-Colomer et al., 2017) by producing a "right-stepping" along the mid-crustal seismicity pattern (Adhikari et al., 2023). However, the precise characteristics of the tear fault and its regional implications are yet to be definitively determined. Considering that the region has remained a high seismicity rate for years (Adhikari et al., 2023), a meticulous examination of the aftershock catalog, specifically focusing on correlations with seismogenic structures, could reveal aspects of the intricate subsurface configuration of the MHT system.

To enhance the delineation of these distinct geological structures, a high-precision aftershock catalog is required, which led us to the utilization of machine-learning-based earthquake detection and location (Zhao M et al., 2023). This method, which mainly enhances the precision of phase picking (e.g., Perol et al., 2018; Zhou YJ et al., 2019; Zhu WQ et al., 2019; Zhu WQ and Beroza, 2019), has gained widespread adoption and demonstrated exceptional efficacy in more accurately locating fault structures (e.g., Tan YJ et al., 2021; Zhou LQ et al., 2022; Duan LF et al., 2023; Ding HY et al., 2023), advancing our comprehension of seismogenesis and fault systems (Feng T et al., 2022; Cai GY et al., 2023; Xu W et al., 2023).

In this study, we applied an end-to-end machine-learning-based earthquake location workflow called LOC-FLOW (Zhang M et al., 2022) to analyze two sets of data covering this region. Through this process, we aimed to adjust the parameterization of LOC-FLOW and generate a high-precision earthquake catalog for the region, particularly a high-precision aftershock catalog of the 2015 Gorkha earthquake. We investigated the spatial heterogeneities in the distribution of seismicity along strike by dividing the region into distinct boxes. We then analyzed the potential relationship between each box and the geological structures, especially the flat-ramp geometry and tear fault. This article summarizes our findings and highlights their implications for future seismic risk assessment in this region.

2.   Data and methods

Our study primarily followed the workflow of LOC-FLOW developed by Zhang M et al. (2022), which consists of machine-learning phase picking using PhaseNet (Zhu WQ and Beroza, 2019), phase association using Rapid Earthquake Association and Location (REAL, Zhang M et al., 2019), and determining the absolute location with HYPOINVERSE (Klein, 2002) and the double-difference relative location with HypoDD (Waldhauser and Ellsworth, 2000). We applied this method to two seismic arrays deployed at different times in the central part of the MHT, within 83°E–87.5°E and 26.5°N–29°N (Figure 1). The processing steps illustrated in Figure 2 were conducted to determine how best to use LOC-FLOW and obtain a high-precision earthquake catalog for this particular region.

Figure  2.  Workflow in this study. Two datasets were used to determine earthquake location, based on the LOC-FLOW method (Zhang M et al., 2022).

DownLoad: Full-Size Img PowerPoint

The data were obtained from two different sources which were the only two temporal seismic arrays covering this region that were available through the Incorporated Research Institutions for Seismology (IRIS). The first dataset was recorded from January to March, 2003, by the Himalayan-Xizang Continental Lithosphere during Mountain Building (Hi-CLIMB, Nábělek et al., 2009) seismic array, serving as a reference for a relatively quiet period. Continuous waveform data were downloaded from 66 stations in the south at a sampling rate of 50 Hz. The other dataset was recorded from July 2015 to April 2016 by the Nepal Array Measuring Aftershock Seismicity Trailing Earthquake (NAMASTE, Karplus et al., 2020), which was deployed after the Gorkha earthquake. We downloaded waveform data from 31 stations in this array, at a sampling rate of 100 Hz or 200 Hz. The distribution of these stations is shown in Figure 1, and the details of these two datasets can be found in Data availability. Of note, not all station data can be downloaded consistently from IRIS. The exact number of stations that one can download data for per day is listed in Table S1 (per month, plotted in Figure S1).

2.1   Data preprocessing

Data preprocessing was essential to ensure the accuracy of the subsequent workflow. This preprocessing involved converting the downloaded miniseed files to sac files, detrending, broad-brand filtering to 0.001–30 Hz, and removing instrument responses from the seismic data. Additionally, some of the NAMASTE data had to be downsampled to achieve a consistent sampling rate of 100 Hz because all the training data of PhaseNet were sampled to 100 Hz (Zhu WQ and Beroza, 2019). Once these preprocessing steps were completed, the datasets were ready for phase picking.

2.2   Earthquake detection and initial location

PhaseNet, an automatic phase-picking technique based on U-shaped convolutional neural networks, was applied as the first step of phase picking. To mitigate the risk of missing phases owing to model uncertainty, a P/S picking probability of 0.3 was set. We obtained 131,925 P-wave picks and 110,315 S-wave picks from 3 months of Hi-CLIMB data and 1,346,501 P-wave picks and 1,077,626 S-wave picks from the 10 months of NAMASTE data.

Next, we associated the picks belonging to the same seismic event. Here, we adopted REAL (Zhang M et al., 2019) to associate phases and obtain an initial location in two steps: grid-based location estimation and simulated annealing refinement. A small grid size was used to enhance the accuracy of this method. An eligible event must have at least three P picks, two S picks, and a total of eight P and S picks. The grid-search method was applied, with the search area centered around the station recording the earliest P phase and covering a horizontal range of 1° and a depth range of 50 km with search intervals of 0.04° and 2.5 km, respectively. The averaged 1D velocity model used for the calculation of travel time in REAL, calculated from the S-velocity model by Xu et al. (2013), is shown in Figure S2. Ultimately, 754 events from the Hi-CLIMB data and 13,656 events from the NAMASTE data were associated. Further relocation was then conducted to mitigate the effects of long-range, low-precision seismic phases in the catalog.

2.3   High-resolution relocation and magnitude calculation

To improve the accuracy of the initial catalog, we determined the location in two steps: the absolute location was obtained using HYPOINVERSE (Klein, 2002) and the relative location was obtained using HypoDD (Waldhauser and Ellsworth, 2000). For the absolute location, we set the initial calculated depth in HYPOINVERSE to 10 km, considering that both the main shock and largest aftershock occurred at mid-crustal depths (10–15 km). Additionally, we assigned weights to the stations based on their epicentral distance, setting a weight of 1 for picks contributed by stations within 50 km and a weight of 0 for picks contributed by stations above 100 km, with a linear decrease between 50 and 100 km, reducing the weighting of remote stations to enhance the accuracy of the absolute location. This process resulted in 14,214 events being relocated by HYPOINVERSE (561 events from Hi-CLIMB and 13,653 events from NAMASTE). We further filtered them and retained 13,062 high-quality events (450 from Hi-CLIMB and 12,612 from NAMASTE) with horizontal and vertical errors within 5 km, station gaps within 270°, and travel-time residuals within 0.5 s, as the final HYPOINVERSE catalog.

We then applied HypoDD for further relocation, leveraging the cluster distribution of earthquakes to refine the locations of these events. The relative location was obtained using the catalog differential travel time and 16 iterations (with NSET=4, damp=70, and the other parameters listed in Table S2). Considering the possible error of the absolute location, we used all absolute location events (13,653) instead of the filtered events (12,612). Finally, we obtained 13,057 high-precision events (388 from Hi-CLIMB and 12,669 from NAMASTE), providing a comprehensive representation of the seismicity during each period. Although the location process included a magnitude calculation (Figure S3), further calculations were required to obtain a more accurate value. The local magnitude was calculated using the following formulas (e.g., Feng T et al., 2022):

Here, A is the maximum amplitude, ${A}_{{\mathrm{N}}}$ and ${A}_{{\mathrm{E}}}$ are the maximum amplitudes of the two horizontal components of the signal after removing the instrument response and then convolving the obtained signal with the Wood-Anderson response. $R(\varDelta )$ is the distance correction factor, which was calculated based on a typical amplitude-magnitude relationship (Bakun and Joyner, 1984).

We calculated the travel times of the HypoDD picks and the magnitudes of the HypoDD catalog (Figure 3). The time-distance curves (Figure 3a) for P-waves (red) and S-waves (blue) show that the wave speed is in good agreement with the average speed we selected for the region (Figure S1). The magnitudes of the HypoDD catalog (Figure 3b) ranged from ${M}_{{\mathrm{L}}}$ 0.14 to 4.88. The magnitude of completeness (${M}_{{\mathrm{C}}}$) is estimated as approximately 1.5, smaller than that of previous studies in this region, such as those by Yamada et al. (2020; ${M}_{{\mathrm{C}}}\approx 2$) and Adhikari et al. (2023; ${M}_{{\mathrm{C}}}\approx 2.5$).

Figure  3.  Results from HypoDD. (a) Time-distance curves and (b) Gutenberg-Richter plot of the cumulative number of events from the HypoDD catalog.

DownLoad: Full-Size Img PowerPoint

3.   Results

3.1   Earthquake distribution

Catalogs obtained using the REAL, HYPOINVERSE, and HypoDD methods are shown in Figure 4. The REAL catalogs (Figures 4a and 4b) show some clustering of the distribution and contain the most events. Although REAL (Zhang M et al., 2019) can be used to refine locations using a simulated annealing method, some events in the REAL catalogs show a latticed distribution. The HYPOINVERSE catalogs (Figures 4c and 4d) display a more centralized distribution than the REAL catalogs. The number of events located far from the clusters was significantly reduced and was primarily influenced by the distance weighting applied to the stations (Klein, 2002). This trend was further emphasized in the HypoDD catalogs (Figures 4e and 4f). Finally, we successfully obtained high-precision catalogs specific to the target region.

Figure  4.  Earthquake distribution of the two data sets with locations obtained via different methods: (a, b) REAL catalogs for Hi-CLIMB data and NAMASTE data, respectively; (c, d) HYPOINVERSE catalogs for Hi-CLIMB and NAMASTE; (e, f) HypoDD catalogs for Hi-CLIMB and NAMASTE. The two red stars in (b, d, f) mark the epicenters of the ${M}_{{\mathrm{W}}}$>7 main shock and aftershock.

DownLoad: Full-Size Img PowerPoint

In the Hi-CLIMB catalogs (Figures 4a, 4c and 4e), which were recorded during a relatively quiet period, earthquakes are relatively infrequent. However, the distribution of certain earthquake clusters is clear. A major vein line runs along the strike of the MHT, parallel to the surface trace of the MCT. At 85°E along this vein line, the events are significantly denser. Additionally, a shallower cluster was observed at the same longitude situated within the southern surface lines of the MBT and MFT. Moreover, the more events occurred in the western region than in the eastern region. These two patterns can be mostly attributed to the distribution of Hi-CLIMB stations.

The spatial distribution of seismic events in the Hi-CLIMB catalogs implies that seismicity during the quiet period was distributed along the MHT’s mid-crustal ramp (Brown et al., 1996; Hauck et al., 1998; Hubbard et al., 2016), which is thought to have played a vital role in the formation of the Gorkha-Pokhara anticlinorium (GPA) and can be traced by the axis of the GPA (Hubbard et al., 2016). Rather than spreading extensively along the dip direction of the MHT, most events are situated along the strike of the MHT. Notably, the two most concentrated clusters observed at 85°E exhibit depths ranging from 10 to 15 km. These lateral clusters also indicate a noticeable discontinuity along the dip direction. Of note, the limited availability of data (overall sparse coverage and an acquisition period of only three months) restricts the comprehensive characterization of other scattered events.

For the NAMASTE catalogs (Figures 4b, 4d and 4f), it is evident that almost all aftershocks were distributed within the rupture zone of the 2015 Gorkha Earthquake, spreading from the main shock to the largest aftershock, as in other studies of aftershocks (e.g., Adhikari et al., 2015; Baillard et al., 2017). Similar to the distribution observed in the Hi-CLIMB catalogs, the clusters spread primarily along the strike of the MHT and were controlled by the MHT flat-ramp geometry (e.g., Hubbard et al., 2016). Specifically, a continuous seismic belt and four additional clusters are discernible. Among the four clusters, the northernmost cluster, which extends further into the Qinghai-Xizang Plateau, exhibits a normal-faulting focal mechanism (Li L et al., 2017), which is beyond the scope of our discussion.

In our catalogs, the along-strike patterns are similar to those from previous studies (e.g., Adhikari et al., 2015; Baillard et al., 2017; Bai L et al., 2019; Yamada et al., 2020) but have clearer details owing to the machine learning-based phase-picking method and the well-designed workflow (Zhang M et al., 2022). Such improvements provide a solid foundation for analyzing the underlying causes of these seismic events. To facilitate a more comprehensive investigation of these heterogeneities, we partitioned the aftershock catalogs into six distinct areas (Figure 5, Table S3). The delineation of these partitions was based on careful consideration of the spatial and depth distribution of each cluster, as well as the seismicity segmentation proposed by other studies (Tiwari et al., 2022; Adhikari et al., 2023). In the subsequent section, we discuss the detailed seismicity within each delineated area.

Figure  5.  Map of seismicity (circles colors correspond to depth) based on NAMASTE data, located using HypoDD. Red rectangles delineate the areas named one to six, which were examined in the subsequent aftershock analysis.

DownLoad: Full-Size Img PowerPoint

3.2   Seismicity of each area

The seismicity of these six areas is shown in Figure 6. Area One (Figures 6a and S4) contains 1,905 events with an average depth of ~10.7 km. This area was primarily affected by the mainshock. The events in this area can be further divided into three smaller clusters: (a) one distributed along the strike of the MHT, with its west side close to the main shock's epicenter; (b) one distributed along the up-dip direction of the MHT, with its north side near the main shock's epicenter; (c) one distributed along the strike of the MHT, with its west side connected to the south side of cluster (b), but with a lower event density than the other two clusters and a different spreading angle (angle of the general trend of seismicity, with respect to the North) to cluster (a). Notably, shallower events are distributed along the junctions of the three clusters. Together, these three clusters potentially show the rupture boundary and propagation path of the mainshock.

Figure  6.  Maps of seismicity for each area. The stars in (a) and (d) indicate the ${M}_{{\mathrm{W}}}$7.8 Gorkha earthquake and the ${M}_{{\mathrm{W}}}$7.3 largest aftershock. More detailed information about these areas is provided in Table S3.

DownLoad: Full-Size Img PowerPoint

Area Two (Figure 6b) contains 1,008 events with an average depth of ~14.5 km. This area was also primarily affected by the main shock and the events form only one seismic cluster. This cluster, connected to cluster (a) in Area One, was also distributed along the strike of the MHT. However, it is characterized by a greater focal depth and is distributed along a different angle, closer to an E-W direction.

Area Three (Figures 6c and S5) contains 2,495 events with an average depth of ~14.2 km. This area was affected by both the mainshock and the largest aftershock. The clusters here can be divided into major and minor branches: (a) The major branch spreads along the strike of the MHT, but is not continuous with the cluster in Area Two. Towards the western side of this cluster, an anomalous densification of seismicity was observed, aside from very little seismicity to the west, indicating the possible presence of a slanting barrier. (b) The minor branch appeared to originate southwest of the major cluster, indicating that the barrier did not hinder the propagation of this minor cluster. The distinctive branching structure observed in this area indicates the likelihood of barrier(s) impeding the westward expansion of the major branch.

Area Four (Figures 6d and S6) contains 5,948 events, with an average depth of ~13.4 km. This area was mainly affected by the largest aftershock and exhibited the highest level of seismic activity among all the areas. Clusters within this area display a complex pattern, leading us to categorize them into two linear and three peripheral clusters. (a) The first linear cluster is distributed along strike of the MHT although at a more SE orientation and its west side is connected to the eastern extremity of the major branch of Area Three. (b) The southern side of the second linear cluster connects to the southeastern part of the preceding cluster and the seismic events were distributed along the dip direction of the MHT. These two continuous linear clusters delineate the rupture boundary of the largest aftershock, and many shallower events occurred at their junction, suggesting an earthquake mechanism similar to that of the clusters observed in Area One, potentially resulting from the two ${M}_{{\mathrm{W}}}$>7 earthquakes. The three peripheral clusters (at ~85°48'E, ~86°00'E, and ~86°10'E, respectively) are situated beyond the linear clusters in the northern part of the region, and comprise deeper events. Unlike those in the linear clusters, these events occurred along the downdip side of the MHT and were not connected to each other.

Area Five (Figure 6e) contains 854 events with an average depth of ~10.5 km. This area is located in northwest Kathmandu. The clusters can be divided into two distinct but connected branches: one distributed almost N-S, and the other along strike of the MHT. Both of them exhibit shallower earthquake distributions than those in Area Two along the same longitude.

Area Six (Figure 6f) contains 221 events with an average depth of ~13.5 km, which is ~3 km deeper than that of Area Five (although they are at a similar distance from the MBT), suggesting a different earthquake mechanism. This cluster exhibits a paucity of seismic events, and most of them are related to the reactivation of some local structures (Baillard et al., 2017), with its shape fitting well with the surface trace of the MCT next to it.

4.   Discussion

4.1   Comparison with previous studies

Several aftershock catalogs of the 2015 Gorkha Earthquake have been created using multiple seismic datasets and various methods (e.g., Adhikari et al., 2015, 2023; Bai L et al., 2016, 2019; Yamada et al., 2020). These catalogs show consistent aftershock distributions. This consistency has endured over the subsequent 5 years, as confirmed by the later seismicity of this region (Adhikari et al., 2023). Our catalogs show good agreement with the large-scale distribution of these previous catalogs and reveal additional details on a small scale, especially in the western part of the catalogs, which not only ensures the reliability of our catalogs but also reveals more geological information related to the Gorkha earthquake sequence.

The use of the machine learning phase-picking method PhaseNet (Zhu WQ and Beroza, 2019) is the biggest advantage of our catalogs over other catalogs (Figure S7), as this method effectively improves the spatial resolution by increasing the number of small-magnitude events. This is reflected in the smaller ${M}_{{\mathrm{C}}}$ of our catalogs and the denser distribution of western clusters than those in previous studies, allowing for a more detailed partitioning discussion and better revelation of certain geological structures (e.g., the tear fault).

4.2   Aftershocks controlled by specific tectonic structures

4.2.1   Step distribution along dip direction controlled by MHT flat-ramp geometry

While previous studies have analyzed the data recorded in this study area using receiver function (Hetényi and Bus, 2007; Nábělek et al., 2009; Caldwell et al., 2013; Duputel et al., 2016) and other seismic imaging methods (Lemonnier et al., 1999; Berger et al., 2004; Robert et al., 2011; Kurashimo et al., 2019), substantial uncertainties persist about the nature of the MHT structures. These include a mid-crustal ramp (e.g., Elliott et al., 2016; Baillard et al., 2017) or two ramps (Hubbard et al., 2016; Wang X et al., 2017), flat-ramp geometry or a duplex structure (Mendoza et al., 2019), and on the specific depth profile of the MHT decollement (which varies among research groups). Additionally, variations in the velocity models and location algorithms employed in different studies have contributed to discrepancies in the depth calculated for earthquake locations. Given these challenges, this study focused on identifying potential subsurface structures based on the distribution patterns of the observed clusters.

To better describe the flat-ramp structure, we created multilocation profiles of this region along the dip direction of the MHT, as shown in Figure 7. Each profile represents a distinct area characterized by different clusters along the dip direction. The subsequent analyses were conducted individually for each region.

Figure  7.  Map and profile views of the aftershock sequence locations obtained from NAMASTE data. (a) Map of seismicity of the region. AA', BB', CC', and DD' indicate the projection lines of profiles shown in (b), (c), (d), and (e), respectively. Earthquakes within ~20 km have been projected onto the profiles.

DownLoad: Full-Size Img PowerPoint

1) The clusters within Area One (Figure 7b) exhibit pairwise connections, wherein the two clusters distributed along the strike are linked by a north-south cluster. The two clusters do not exhibit significant depth variations. The northern along-strike cluster is uniformly distributed between depths of 5 and 15 km, primarily influenced by the main shock, and is situated on the hanging wall of the MHT.

2) The two clusters in Areas Two and Five are shown in Figure 7c. The earthquake cluster of Area Five, in the south, is situated at a depth of approximately 5 km shallower than that in Area Two, displaying a distinct step distribution. Considering the different physical properties of flat and ramp geometry (e.g., Ramsay, 1992; Pedrera et al., 2012), we infer that these two clusters are located in regions where the slope of the MHT undergoes changes, most likely a ramp-flat or flat-ramp transition. We suggest that both clusters are located near ramps and that the elongated part between them corresponds to a flat decollement of the MHT flat-ramp geometry in this particular area. An enlarged schematic diagram of flat-ramp geometry is shown in Figure S8. The observed distribution pattern, characterized by two clusters at different depths resembling stairs, supports the presence of double ramps (Hubbard et al., 2016; Wang X et al., 2017) in this region.

3) The two clusters in Area Six and the junction of Areas Three and Four are presented in Figure 7d. In contrast to that in Figure 7c, there was no systematic increase in the depth or step distribution between the two clusters. Considering the previous description of Area Six, the structure of the MHT cannot be used to explain the localized depth distribution in the southern part of this profile.

4) The distribution of clusters in Area Four also exhibits a prominent stepwise pattern in the profile (Figure 7e). The linear cluster in the southern part shows a uniform depth distribution concentrated between approximately 10 and 15 km. It broadens and densifies towards the center, with a major depth range of ~5–17 km. In contrast, the distribution of the three discrete clusters in the northern region is primarily concentrated below ~15 km on the profile, lacking significant development towards shallower depths. However, the steps in this region are more closely interconnected, unlike those in Areas Five and Two (Figure 7c), where a clear flat decollement exists between the two clusters. The linear clusters of Areas Two and Four are both distributed near the GPA, suggesting that the linear cluster in Area Four may also be influenced by a ramp, similar to that in Area Two. The three deeper northern clusters (enlarged in Figure 8) in Area Four might also be attributed to three localized, smaller ramps further north with less separation from the main ramp compared to Figure 7c. Notably, there were some shallower events (~10 km or above) in the north, which may be related to the occurrence of the ${M}_{{\mathrm{W}}}$7.3 aftershock in this area.

Figure  8.  Map (a) and profile views of Area Four’s aftershock sequence locations, obtained from NAMASTE data. The three black dashed rectangles A, B, and C indicate the projection areas of profiles shown in (b), (c), and (d), respectively.

DownLoad: Full-Size Img PowerPoint

To examine the three peripheral clusters of Area Four, we created three smaller-scale profiles through Area Four along the dip direction of the MHT, as shown in Figure 8. These profiles (Figures 8b–8d) show a step pattern similar to that shown in Figure 7e, although the clusters appear more separated. Without higher-spatial-resolution models or images, it is difficult to uniquely determine the specific form of the MHT. For our catalogs, we plotted a profile of the eastern end of the aftershocks, as done by Mendoza et al. (2019), and obtained a similar distribution with three possible steep structures inside the large cluster in Figure 8d (an enlarged picture can be seen in Figure S9), which may represent a series of steeply dipping imbricate faults within a duplex structure. However, the other profiles (Figures 8b, 8c as well as supplementary figures in Mendoza et al. (2019)) do not show such a noticeable pattern, suggesting different mechanisms along the strike in Area Four. We infer that the flat-ramp geometry still plays a major role here, whereas seismicity in the eastern part may also be affected by a duplex structure (Mendoza et al., 2019).

4.2.2   Segmentation along strike caused by a tear fault

In addition to the impact of the MHT flat-ramp geometry and duplex structure, a tear fault also influences the distribution of aftershocks (Baillard et al., 2017; Hoste-Colomer et al., 2017; Adhikari et al., 2023). By analyzing the seismicity of each area, a distinct segmentation emerges along the strike between Areas Two and Three (around 85.4°E). We therefore created a profile view along strike through the study area to ascertain the relationship between each segment (Figure 9). Events that occurred in different areas are denoted by different colors.

Figure  9.  Map (a) and profile (b) views of aftershock sequences in Areas One, Two, Three, and Four. Black dotted line indicates the projection line. Red dotted box outlines the enlarged area shown in Figure 10. Red, yellow, blue, and green circles indicate the earthquakes occurring in Areas One, Two, Three, and Four, respectively.

DownLoad: Full-Size Img PowerPoint

The spreading angles of the main clusters along the strike differs between Areas One and Two as well as between Areas Three and Four, even though these areas are connected. However, the spreading angle of Area One closely resembles that of Area Four. Although the clusters in Area Two and the major branch in Area Three are not connected, they exhibit similar spreading angles. Furthermore, in Areas One and Four, the clusters are distributed over a relatively constant depth range, aligning predominantly along a horizontal line. Conversely, in Areas Two and Three, the depth range of the clusters display evident changes, illustrating a sinking trend from the two ends towards the central point.

Enlarged maps and profiles (Figure 10) were created for the clusters in Areas Two and Three in two directions: along the strike of the MHT and along the possible tilt direction of the barrier. The clusters in Areas Two and Three are clearly distributed at a similar angle but also over similar depths, indicating their potential affiliation with the lower edge of the locked MHT (e.g., Avouac et al., 2015; Grandin et al., 2015; Hoste-Colomer et al., 2016). This observation supports the notion that the slanting barrier was created by a strike-slip fault, owing to the similarities between these two clusters (depth, distribution angle, and both as the lower edge of the locked MHT). This strike-slip fault is most likely to be further identified as a tear fault, given its location within a thrust fault system (Dahlstrom, 1970) and is also supported by kinematic analysis (Kumahara et al., 2016; Hoste-Colomer et al., 2017). The tear fault contributes to the segmentation along strike in the distribution of the aftershock catalog (e.g., Baillard et al., 2017; Hoste-Colomer et al., 2017; Adhikari et al., 2023).

Figure  10.  Map (a) and profile (b) views of aftershock sequence in Areas Two and Three obtained from NAMASTE data. AA' and BB' indicate the projection lines of (b) and (c), respectively. Yellow and blue circles indicate locations from clusters in Areas Two and Three, respectively.

DownLoad: Full-Size Img PowerPoint

We propose that the entire seismic belt, which is composed of along-strike-distributed clusters in these four areas, may have originally constituted a single geological structure representing the lower edge of the locked MHT (e.g., Avouac et al., 2015; Grandin et al., 2015; Hoste-Colomer et al., 2017), prior to the emergence of this tear fault. This geological structure developed at a stable angle along the strike of the MHT (Figures 11a and 11c). However, possibly owing to evolution of the MHT, the structure exhibited non-uniform stress, resulting in greater stress accumulation in the middle section than on the ends. This may also have led to the appearance of along-strike heterogeneities such as a central decollement, which separates the central ramp from the deep ramp (Hubbard et al., 2016). Consequently, the structure underwent compression in the middle, gradually bending downwards, and eventually rupturing to form a tear fault. This process changed the form of this geological structure, which is reflected in the seismicity distribution. Specifically, the earthquake clusters are deeper at the junction of Areas Two and Three (Figure 11d), and their spreading angles have undergone significant changes.

Figure  11.  A simple schematic diagram of the tear fault. (a, b) shows the influence of tear fault formation on this geological structure horizontally. (c, d) shows the possible vertical impact on this geological structure.

DownLoad: Full-Size Img PowerPoint

By analyzing the overall cluster distribution, we discerned the dislocation pattern associated with this tear fault. Close to the tear fault, the spreading angle in Area Two tilts towards the north in comparison to Area One, whereas the spreading angle in Area Three tilts towards the south relative to Area Four. This suggests that the tear fault was likely a right-lateral strike-slip fault, which altered the shape of the geological structure. As the middle section is torn, the two sides of the structure were separated, eventually causing the west side to bend northward and the east side to bend southward (Figure 11b).

4.3   Implications for future seismic hazard assessment

From our analysis of the aftershock catalog, it is clear that the aftershock seismicity of the 2015 Gorkha Earthquake was mainly controlled by the flat-ramp geometry along the dip direction and a tear fault along strike, which could be used for future hazard assessment in this region. The catalog obtained here can be compared with those of previous or future studies (e.g., Hubbard et al., 2016; Baillard et al., 2017; McNamara et al., 2017; Arora et al., 2017; Bai L et al., 2016, 2019; Adhikari et al., 2023), improving our ability to estimate the location of future earthquakes in different regions and their possible rupture characteristics (Dal Zilio et al., 2021).

In previous studies (e.g., Tan YJ et al., 2021; Zhou LQ et al., 2022; Feng T et al., 2022), machine learning-based methods have proven useful in retrieving a greater number of smaller events and revealing more detailed structures. In our study, these smaller events enriched the details of the earthquake clusters, in particular, the western clusters in Areas One and Two provided clearer evidence of the tear fault, as discussed above. Considering that previous studies on tear faults in this region have focused on individual focal mechanisms (Hoste-Colomer et al., 2017) or discontinuities in seismic density maps (Baillard et al., 2017; Adhikari et al., 2023), our analysis contributes to our understanding of this tear fault. In addition, our interpretation of the seismic events in the catalogs and their relationship to deep geological structures accounts for both the flat-ramp geometry and the duplex structure in this area. The enriched catalog details indicate that these two structures may control the aftershock seismicity at different locations along the MHT, with the flat-ramp structure playing the major role. Improved identification and understanding of the structures that comprise fault systems are helpful for future seismic hazard assessments and long-term earthquake predictions (e.g., Stevens et al., 2018; Chamoli et al., 2021; Mittal et al., 2019; Li ZF, 2021).

5.   Conclusions

In this study, we applied LOC-FLOW (Zhang M et al., 2022), a machine learning-based earthquake location workflow, to two sets of seismic data collected by the Hi-CLIMB and NAMASTE arrays and obtained 13,057 high-precision earthquakes. Based on these results, the following conclusions were drawn:

1) The 388 events that occurred between January and March 2003 indicate a period of relatively low seismicity and were mostly distributed along strike of the MHT, demonstrating a strong relationship with the mid-crustal ramp.

2) The 12,669 events between July 2015 and April 2016 form an aftershock catalog, with similar distribution patterns to catalogs created in previous studies (e.g., Baillard et al., 2017; Adhikari et al., 2023). However, the events in this catalog were more strongly correlated with certain geological structures of the MHT. Through zoning analysis, we further demonstrated the differences in seismic activity among different areas.

3) Seismogenic structures of the MHT: based on a detailed analysis of seismicity maps and profiles, we suggest that the flat-ramp geometry controls the overall pattern and depth distribution of seismicity in this region, and that a tear fault is the primary cause of the sharp change along strike at approximately 85.4°E. Moreover, a duplex structure may have affected the seismicity in the easternmost part of the study region.