The first results from the data analysis reveal the peculiar behavior of the Himalchuli swarm, both persistent over weeks and with a possible strong seasonal modulation. With 6793 M_L>1.8 events, this is the largest swarm ever recorded in the Himalayas. The onset of the seismicity in 2017 is sharp and not following at all the typical mainshock-aftershocks pattern.

Indeed, after a first M_L 2.8 earthquake on August 7, 2017 at 00:29:21 UTC, the swarm expanded quickly over a surface similar to the final extension of its overall seismicity. The beginning of this seismic sequence do not coincide with the occurrence of a large earthquake at regional or teleseismic distance. We found therefore no candidate for a dynamic triggering by long-period surface waves, a phenomenon described further west in Kumaon-Garhwal Himalaya, with local microearthquake activity triggered by the surface waves of the 2007 Mw 8.5 earthquake (Mendoza et al., 2016). The number of events per day increased from 1 to 294 in 7 days, contributing then to a dense seismic activity near the barycenter of the swarm (Figure 2 and Supplementary Figures 2, 3). The rate of earthquakes decreased more progressively than it rose until early September 2017. A few additional earthquakes occurred from 16 September to November 2017 (Figure 2 and Supplementary Figure 2). The seismic activity rate reached a low level with 1 event every 10 days until March 2018 and the beginning of a new period of activity (Figure 2 and Supplementary Figure 2). The rate of earthquakes increased more progressively from the May 17, 2018, reaching a peak of activity on July 14, 2018, with 15 earthquakes that day. This episode was followed by a decrease of the seismicity that ended on September 12, 2018. Then, almost no earthquakes were detected until the summer 2019, when sporadic local earthquakes occurred (Figure 3A). The higher seismicity rates of the swarm during the summer for these three successive years cannot be attributed to a temporal decrease of the seismic noise. Indeed, the high frequency seismic noise at the NSC stations is higher in summer than in winter due to river and landslide activity (Burtin et al., 2009), leading to less event detection (Bollinger et al., 2007). The largest event, with M_L 3.8, occurred on November 5, 2017.

The seismic moment released by the crisis was evaluated from the M_w converted from the M_L of the earthquakes obtained at NSC, using the following relation (Ader et al., 2012): M_w=0.84M_L+0.21. We obtained that the 5992 events recorded in 2017 progressively released 9.92 × 10¹⁶ Nm, the 759 events in 2018 released 1.19 × 10¹⁶ Nm, while 2019 events released only 5.13 × 10¹³ Nm. The total seismic moment released reached then 1.11 × 10¹⁷ Nm, a value equivalent to a single earthquake of M_w 5.3.

The earthquake catalogues also show that the seismicity spreads within a single-continuous patch (Figures 1, 2, 3A). However, the genuine size and orientation of the spread are not sufficiently well constrained in the original NSC catalogue. Indeed, the addition of one or the other temporary stations (PRK, between the 21 and 28/08/2017 or JL, between the 02 and 18/09/2017), which efficiently closes the azimuthal detection gap, reshapes the cluster. The synthetic tests performed for every network geometry available during the seismic crisis help to quantify the hypocentral location biases (Supplementary Material “Synthetic tests” and Supplementary Figure 6). They demonstrate the robustness of the catalogue generated when the network was complemented by these temporary stations closing the azimuthal gap. This relocated catalogue, which contains only the 20% of best located earthquakes, is illustrated in Figure 3B.

The structure of the seismicity deduced from the relocated catalogues is significantly different from the original one (Figure 3A). The seismicity swarm does not spread over nor is elongated along a NE-SW axis, but is organized within a narrow (2–4 km-thick) steeply northward dipping, 7 × 15 km² region located at the foot of the Himalchuli summit (Figure 3B). The difference observed in comparison with the original catalogue is expected, since the shape and width of the original seismicity spread were largely controlled by the network geometry (a result illustrated by synthetic catalogue testing described in Supplementary Material and Supplementary Figure 6). The difference is also not a bias related to the short period and limited number of events covered by the temporary stations. Indeed, the density map of events considered for the relocation shows similarities with the density map of all earthquakes happening in 2017 or 2018, the location and shape of the area of highest density remaining similar (Supplementary Figure 3). This suggests that the geometry described by a small number of well-relocated events may be representative of the overall distribution of the swarm. The events in 2018 and 2019 could also line up with the same inferred fault, no migration of the densest area of the seismicity or change of the shape of this region being detectable (Supplementary Figure 3). The depth of the relocated events ranges between 0 and 28 km, when considering all events, most of the best located events (95%) falling between 0 and 15 km. An average seismic slip of 2.85 cm on such a 7 × 15 km² plane would have been enough to release the total seismic moment of the swarm (1.11×10¹⁷ Nm). However, the swarm geometry defined by the microseismicity remains associated with large uncertainties, while the kinematics on the fault remains uncharacterized, limiting the possibility for a convincing interpretation of this seismicity taken out of its geological or seismotectonic context. We therefore attempted at determining the focal mechanisms of the largest earthquakes of the seismic swarm in order to benefit from a fault plane orientation and slip vector. The attempts at waveform inversion at regional distances failed because all earthquakes are too small (below M_w 3) for being studied at KKN, the nearest permanent broadband station installed 90 km from the cluster. In parallel, we attempted at determining focal mechanisms with first motion polarities of the direct P-wave. For most of the largest earthquakes selected, the signals are not very impulsive and the polarities are difficult to ascertain. Furthermore, for all of them, the distribution of the polarities on the focal sphere is far from optimal, due to the large azimuthal gap, the range of hypocentral distances covered by the stations and the shallow depths of the earthquakes. The events of M_Lv 3.5 and 3.1 (resp. M_w 2.8 and 2.5) that occurred on August 26 and 29, 2017, respectively (and were recorded at PROK), are among the best recorded for determining a focal mechanism (see Supplementary Material). However, the result we obtained for the two are not consistent, precluding the validation or invalidation of some possible geometries based on the observed polarities (Supplementary Figure 7).

The updip-end, shallow trace, of the swarm defined by the best located earthquakes, recorded during the few days when the local station of PROK was available, falls precisely at the toe of Himalchuli Southeastern face, along a glacial valley. They did occur on a steep plane. The geological survey of the region by Colchen et al. (1986) revealed no active fault trace but a large slab of High Himalayan Crystalline (HHC) units, slightly dipping to the North. The slab is associated with almost East-West structural fabrics. The bedding of the local quartzite/gneisses, complemented by the bedding of the Lesser Himalayan unit southward, helped tying a local balanced cross-section (Figure 3A and Figure 4A). Most, if not all, of the seismic swarm develops through the Formations I and II of the HHC, above a change of dip of the units in the hanging-walls of the Main Central Thrust (MCT) and MHT fault systems (Figure 4A). A structure activated with the same orientation as the main structural fabrics is likely, but the activation of a structure across the main fabrics remains possible. To our knowledge, the valley along the south-eastern slope of the Himalchuli was not the focus of any extensive tectonic survey, and its morpho-glacial activity probably precludes the preservation of the neotectonic landforms of active faults. Without an opportunity to survey the area, we visually inspected the high resolution GoogleEarth™ and ArcGIS™ Wayback imagery of the world images available in the vicinity of the seismic swarm before and after the Gorkha earthquake and the first swarm sequence (Supplementary section Imagery and Supplementary Figure 8). We checked in particular possible ∼NS and NE-SW structures consistent with the orientation of the closest active normal faults, mapped 60 km to the North-East in Gyirong graben (Armijo et al., 1986). This examination confirmed the main East-West orientation of the rocks bedding, schistosity and/or structural fabrics (Supplementary Figure 8), but was not relevant to identify any active fault trace in the area, along or across fabrics within that area, highly covered by snow and ice.

This shallow seismicity swarm of more than 6500 earthquakes with M_L>1.8 located at the foot of the High Himalayan range is unusual, falling in a region where swarms had never been recorded since 1994. The clustered seismicity, however, developed 30 km North of the 2015 M_w 7.9 Gorkha earthquake epicenter in a region prone, since the main shock, to large variations of strain and stresses. We found no evidence of dynamic triggering, founding no large earthquake possibly associated with the onset of the swarm.

The earthquake ruptured a 140 × 50 km² locked fault patch to the South-East which was responsible for the subsidence of the high Himalayan range (e.g., Avouac et al., 2015; Grandin et al., 2015; Figure 4A). This subsidence locally exceeded 1 m. Despite lower values in the Himalchuli area, the coseismic strain and stress changes applied to the geological medium were likely significant (Figure 4A). Indeed, the peak of coseismic slip on the MHT was observed 80 km to the southeast of the study area, exceeding 5 m (e.g., Avouac et al., 2015; Grandin et al., 2015). A displacement of 1–2 m at depth is also estimated on the fault segments 10–20 km to the South of the Himalchuli. The discrepancy between the source models and the relatively low resolution of the slip accommodated at depth on the midcrustal ramp of the MHT, precludes any precise estimate of the stress buildup at short distances. Nevertheless, the whole hanging wall of the MHT was necessarily considerably affected by the main shock in this region of high slip gradients and localized extension (Figure 4A). The earthquake was immediately followed by deep afterslip, poroelastic rebound and viscous relaxation (Figure 4B) (e.g., Gualandi et al., 2017; Zhao et al., 2017; Wang and Fialko, 2018; Jiang et al., 2019; Jouanne et al., 2019; Tian et al., 2020; Liu-Zeng et al., 2020). The deep afterslip, within the 3 months following the mainshock, exceeded 1 cm and reached locally as much as 5 cm after a year, a value which corresponds to 2.5 years of interseismic loading (Gualandi et al., 2017; Zhao et al., 2017; Ingleby et al., 2020; Liu-Zeng et al., 2020; Tian et al., 2020).

Knowing the large postseismic relaxation and afterslip of the Gorkha earthquake, in the following, we discuss different types of forcing that may modulate, induce or trigger the seismicity in the Himalchuli swarm and the mountain collapse including surface hydrological loading and fluid migration.

The strain and stress rates that resulted from the early postseismic behavior of the fault were significantly larger than those induced by surface hydrological loading of the Earth. However, these postseismic effects dampen with time as the tectonic system relaxes, becoming locally less influential on the strain and stress within the upper crust than the hydrological surface loading. Indeed, in the Himalayas, surface mass redistribution of continental hydrology induces measurable transient and seasonal deformation. Seasonal amplitudes of regional horizontal and vertical displacements, as recorded by cGPS stations located in the high range, reach a few mm to a few cm (Bettinelli et al., 2008; Chanard et al., 2014; Larochelle et al., 2018) (Figure 5). Seasonal hydrology not only induces seasonal strain but also seasonal stress variations at depth on the fault systems. These stress variations are suspected to be responsible for the seasonal modulation of the midcrustal seismicity along the MHT (Bollinger et al., 2007; Bettinelli et al., 2008; Ader and Avouac, 2013).

Here, we quantitatively estimated the strain and stress induced by regional and local surface loading (see Supplementary Methods 1, 2). We first take advantage of satellite gravity measurements from the Gravity and Recovery Climate Experiment (GRACE) to quantify spatial and temporal variations in regional surface loading. We then calculate the full displacement field characterizing the solid-Earth’s response to surface variations in hydrological loading inferred from GRACE, using a spherical layered elastic Earth model (Chanard et al., 2018). Our calculation is validated by comparison with observed displacements at cGPS stations (Supplementary Figure 9), and then used to evaluate stress perturbations on steeply dipping fault planes N080-N100E, dipping 70N-90N, a fault geometry consistent with the pattern of the swarm (Figure 3) and the main structural fabrics of the area. We compute change in Coulomb stress (CFS) as an indicator of the fault susceptibility to failure under annual stresses derived from GRACE. This shows a maximum of CFS in spring and a maximum rate of CFS at the onset of the monsoon, with amplitudes reaching up to 800 Pa (Figure 5, Supplementary Figures 10–12).

The spatial resolution of the GRACE data is only of a few hundred of km (∼300 km). This leads to averaging loads over large areas, resolving the large wavelengths of the surface forcing on the faults. This may lead in this area to a significant mis-determination of the static stress change applied on shallow faults, which are also sensitive to localized loads, such as short wavelength snow accumulation on top of faults. Indeed, during summer monsoon, the freezing-thawing altitude reaches nearly 6000 m. The precipitations under this altitude favor rapid snow melt. Water is evacuated into rivers, while a thick snowpack is formed above 6000 m. According to the precipitations recorded at Gharedunga and Chhekampar, two meteorological stations at close distance from the Himalchuli (Figures 1, 5), surface load of the summer monsoon due to snow accumulation could reach 1 m -water equivalent- above 6000 m in the hanging wall of the suspected northward dipping fault, while no snow accumulates in the footwall which remains below the zero isotherm. This leads to a potential large stress increase on the fault plane. This would favor the normal fault failure during the summer monsoon. To test this hypothesis, we compute stress fields induced by 1 m snow loads localized on the surface of the hanging wall, over the Himalchuli or Manaslu mountain ranges using the Boussinesq solution (Boussinesq, 1878; see Supplementary Methods 2). Supplementary Figure 13 shows the two settings of snow surface load at 6 km on the Himalchuli or Manaslu, as well as the resulting depth profiles of the largest stress component σ z z .   We then evaluate the CFS at a depth of 10 km for an East-West steeply dipping (80°N) fault plane (see Figure 5, Supplementary Methods 2 and Supplementary Figure 13). Supplementary Figure 13 shows results of CFS as a function of the distance to the center of snow load, with a maximum amplitude reaching up to 600 Pa, 10 km BSL, below the load for the combination of both the Manaslu and Himalchuli snow loads. A larger load in term of water equivalent, or a depth shallower than 10 km BSL will lead to values that can overcome the 2 kPa (see Supplementary Figures S13a,b) which correspond to the CFS variations suspected to modulate the secular stress buildup of a few kPa/yr and its associated seismicity along the MHT at midcrustal depth (Bettinelli et al., 2008).

This is not surprising, snow loads having been suspected to be the forcing which modulate seismicity elsewhere, inhibiting the seismicity in winter in Japan on shallow dipping thrust faults and strike slip faults parallel to the snow loads (Heki, 2003). The tectonic setting here is different from the Japanese case, because of the orientation of the geological structures and their kinematics, but the effect could be similar in term of the amplitude of the Coulomb stress modulation on the fault system. Note that if the receiver faults are normal faults oriented NS or NE-SW, steeply dipping toward the highest topography to the West, the results would have not been drastically different.

Note that this local effect adds to the seasonal modulation of the stress by regional forcing previously described with a time structure favoring earthquakes in April and enhancing the seismicity rate in summer, the seismicity rate being expected to vary in proportion to the stress rate (Figures 4B,C) (Bettinelli et al., 2008).

To decipher the processes behind the seismic swarm, another important information is that large carbon dioxide (CO₂) emissions and hydrothermal disturbances followed the Gorkha earthquake at the front of the High Himalayan range (Girault et al., 2018). Several plausible causative scenarios for the failure and seismogenesis of the area can then be considered.

We consider a scenario involving the progressive collapse of the mountain under its own weight along a normal fault system. Indeed, the basement of the Himalchuli-Manaslu, which is the highest mountain group affected by subsidence along the Gorkha rupture, was probably weakened due to instant stretching and damages resulting from the mainshock (Figure 4A). Although not mapped by previous geological surveys, possibly hidden under the local glacier, a northward dipping fault is highly probable and both consistent with the shape of relocated swarm of earthquakes (Figure 3B and Figure 4D) and the local stress field imposed by the high topography. The failure was probably delayed by the deep afterslip on the MHT, responsible for clamping the steep northward dipping normal fault (Figure 4B). During this period of afterslip at depth, the fault is also solicited by the seasonal variations of the surface loads (Figure 4B for summer 2015 and 2016 and 4C for fall and winter 2015 and 2016). In August 2017, the normal fault is activated at the peak of local seasonal stresses (Figure 5C). During the two years that follow, the seismicity is then modulated by the stresses related to the surface hydrological loads (Figure 5C). The transient slip and swarm are reactivated in two stages in 2018: 1) at the peak of the regional Coulomb failure stress applied on the fault, involving a first increase of the seismicity rate in spring, and 2) during the summer monsoon, (a) when the local load in the hanging wall of the normal fault contributes to facilitating the rupture and (b) when the stress rate is optimal. The time structure of the sporadic events recorded in 2019 suggest that this scenario repeated two years after the initiation of the swarm (Figure 5C).

These mechanisms could be possibly accompanied by fluid migration that affects fracture strength. This would explain the high b-value obtained (Figure 2), the seismicity exhibiting a deficiency of large earthquakes, a signature often associated with the presence of fluids (e.g., Murru et al., 1999; Bachmann et al., 2012). In addition to pore pressure changes which could be the driver of the failure, alternative scenarios can involve changes in poroelastic stresses contributing to the generation of induced seismicity. Instead of resulting from the sole fault clamping by afterslip, the delay between the Gorkha mainshock, and/or an initial fluid injection, and the beginning of the swarm development could be due to the time necessary to increase enough the fluid pressure for triggering seismicity. These delays between injection and earthquake realizations have been illustrated elsewhere (e.g., Zhai et al., 2019). The seasonal deformations that follow the first swarm realization in 2017 can lead to repeated episodes in which the rock is hydraulically fractured due to additional poro-elastic stresses.

Several examples of seismicity driven by fluids involved deep CO₂-rich fluids upwelling through fractured rocks (e.g., Cappa et al., 2009; Weinlich, 2014; Ingebritsen et al., 2015; Miller et al., 2004; Chiodini et al., 2020). In the Nepal Himalayas, CO₂ produced by metamorphic reactions is released mainly in the MCT zone in the dissolved form as hot spring waters and in gaseous form as tectonic fumaroles and diffuse degassing structures (e.g., Girault et al., 2014; Figure 1). The 2015 M_w 7.9 Gorkha earthquake strongly affected Central Nepal triggering significant perturbations of the whole CO₂ emitting region (Girault et al., 2018): increases and decreases of gaseous CO₂ emissions at diffuse degassing structures and of dissolved CO₂ emissions in springs, up to a complete cessation of flow in the Trisuli valley (Inset map Figure 1).

Other valleys, less than 30 km away from the swarm, such as the Budhi Gandaki and Marsyandi valleys, also release significant CO₂ and are at closer distances (Figure 1). Repeated CO₂ flux measurements at the Bahundanda hydrothermal system (Inset map Figure 1) during active and inactive periods of the swarm did not show any clear change. The available time-series of soil-gas radon-222 concentration at this site show variations associated with the monsoon and steady signals during winter. Carbon isotopic values of CO₂ (δ¹³C) in the Budhi Gandaki valley remain relatively similar: for gaseous CO₂ emissions, −5.4 ± 0.1 ‰ in January 2017 and −5.6 ± 0.1 ‰ in January 2018 at Khorlabesi, −3.0 ± 0.1 ‰ versus −3.1 ± 0.1 ‰ at Machhakhola; for dissolved CO₂ in hot spring waters, −0.5 ± 0.2 ‰ versus +1.1 ± 0.1 ‰ at Machhakhola. In the Marsyandi valley, at Bahundanda, the δ¹³C values of gaseous CO₂ emissions (−2.2 ± 0.1 ‰ in January 2017, −3.6 ± 0.1 ‰ in August 2018) and of dissolved CO₂ in hot spring waters (+0.2 ± 0.2 ‰ versus +0.9 ± 0.1 ‰, and −0.2 ± 0.3 ‰ versus +0.6 ± 0.1 ‰) show a slight difference possibly due to an increase in CO₂ degassing or to monsoonal rain effects (Girault et al., 2018; C. France-Lanord and P. Agrinier, private communications). At this stage, we cannot identify any clear excursion of carbon isotopic values of CO₂ in gas and hot water at these hydrothermal sites during the period covered by the swarm.

Thus, no hydrothermal changes at known hydrothermal sites nearby were noticed in measurements available during the specific period of activity of the Himalchuli swarm. However, the potential presence of CO₂ emission zones or CO₂-rich springs in the immediate vicinity of the swarm has not been investigated yet. Despite the very high b-value, the particular behavior of the seismicity, the presence of numerous fluid/gas vents in the region affected by the seismicity, there are no other compelling evidence for a fluid-driven swarm activity.

Whereas some of these plausible forcing generate unique seismic temporal structures, or particular seismic spatial patterns, their respective contributions to the triggering of the seismic events emphasized here are difficult to further resolve. Indeed, we were unable to resolve any migration of the seismicity front, nor to determine the volume and time of fluid injection possibly triggering the observed seismicity, as well as to overcome the absence of well-constrained focal mechanisms. Finally, while we lack time-series of near-field InSAR and cGPS data covering the whole crisis, this seismic cluster located below the highest Himalayan summits and found in the trace of the Gorkha earthquake, at the edge of the western Nepal seismic gap, offers unique and challenging observations.