# Annual variations in Nepalese seismicity: b-values and seismicity rates

Martin Colledge, Kristel Chanard, Clara Duverger, Alexandre Schubnel, Lok

Bijaya Adhikari, Laurent Bollinger

► To cite this version:

Martin Colledge, Kristel Chanard, Clara Duverger, Alexandre Schubnel, Lok Bijaya Adhikari, et al. Annual variations in Nepalese seismicity: b-values and seismicity rates. Geophysical Journal International, 2025, 242 (3), pp.ggaf259. 10.1093/gji/ggaf259. hal-05340942

HAL Id: hal-05340942

https://hal.science/hal-05340942v1

Submitted on 31 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Annual variations in Nepalese seismicity: *b*-values and seismicity rates

Martin Colledge<sup>1,2</sup> Kristel Chanard<sup>1,3</sup> Clara Duverger,<sup>2</sup> Alexandre Schubnel,<sup>1</sup> Lok Bijaya Adhikari<sup>4</sup> and Laurent Bollinger<sup>2</sup>

Accepted 2025 June 29. Received 2025 June 23; in original form 2025 February 7

## SUMMARY

Seasonal variations in Nepalese seismicity have been reported with varying degrees of confidence. We re-investigate these claims by analysing 20 yr of Nepalese seismicity before the 2015 Gorkha earthquake, as detected by the Nepalese national network, and focusing on earthquakes located along the eastern and central sections of the Nepalese Main Himalayan Thrust. Using several declustering techniques, we find no statistically robust evidence of seasonal seismicity in the studied record, regardless of magnitude threshold above completeness. This suggests that previously reported seasonality may be restricted to the western section of the Nepalese orogeny, may be an artefact, or may indicate that nucleation times of earthquakes are longer than the year. We also investigate potential annual variations in the Gutenberg–Richter *b*-value, given its recent observed modulation by transient stressing. Additionally, we use large-scale mass redistribution estimated from the monthly gravity field retrieved from the Gravity Recovery And Climate Experiment and Follow-On (GRACE/-FO) missions, to resolve stress variations at depth induced by transient surface loads. We find that the mean annual *b*-value peaks when seasonal Coulomb stress rates reach their minimum value at the height of the summer rainy season. When considering the combined effect of tectonic and seasonal loading, this corresponds to a recurring period of stress reversal, when Coulomb stress momentarily decreases. This suggests that periodic clamping of the Main Himalayan Thrust reduces the likelihood of earthquakes growing to larger magnitudes in accordance with hierarchical rupture models. The susceptibility of *b*-value to stress variations of roughly 0.1 points.kPa⁻<fnref n="1" /> is consistent with recent estimates of *b*-value sensitivity to transient loading, although it remains high when compared to the stress-dependence associated with both static differential stress, and with long-term evolution during the seismic cycle. This discrepancy points to the large impact of stress transients on the dynamics of seismic rupture.

**Key words:** Transient deformation; Asia; Statistical seismology; Fractures, faults, and high strain deformation zones.

## 1 INTRODUCTION

The effect of transient stress changes on the occurrence of earthquakes remains a debated question, despite earthquake periodicities being a recurring subject of study since the early days of seismology (e.g. Merian 1834; Perrey 1861; Schuster 1897; Davison 1938). The spatiotemporally clustered nature of seismicity complicates the problem, requiring a high-quality catalogue, a precise declustering strategy and either a tool to quantify periodicities, or stress estimates to compare periodicities to.

Recent improvements in seismicity catalogues resolutions and their extended durations have enabled studies over the last 20 yr to identify earthquake periodicities at annual frequencies (e.g. Heki 2003; Dutilleul et al. 2015; Zhan & Shearer 2015; Johnson et al. 2020; Hsu et al. 2021; Pantobe et al. 2024; Ueda et al. 2024), and tidal frequencies (e.g. Cochran et al. 2004; Métivier et al. 2009; Wilcock 2009; Wang et al. 2022). These studies suggest that stress variations on the order of the kilopascal can effectively modulate seismicity rates. Additionally, recent studies have also shown dependencies of the magnitude-frequency distribution on transient

<footnote n="1"><sup>1</sup>Laboratoire de Géologie, CNRS, ENS, 75005 Paris, France. E-mail: colledge@geologie.ens.fr</footnote>

<footnote n="2"><sup>2</sup>CEA, DAM, DIF, F-91297 Arpajon, France</footnote>

<footnote n="3"><sup>3</sup>Université Paris-Cité, IPGP, CNRS, IGN, F-75005 Paris, France</footnote>

<footnote n="4"><sup>4</sup>Department of Mines and Geology, National Seismological Center, Kathmandu 44600, Nepal</footnote>
stresses of earthquakes (Ide et al. 2016; Tan et al. 2019), laboratory acoustic emissions (Colledge et al. 2023) and in numerical quasi-static models (Pétrélis et al. 2021), with seemingly much higher susceptibility to transient stress changes than to well-established static stress dependencies (Scholz 1968, 2015).

In the Himalayas, large hydrological loads driven by the Monsoon cycle cause significant annual stress variations along the seismo-genic thrust structures. This annual cycle has led to the investigation of potential seasonal variations of seismicity. Bollinger et al. (2007) demonstrated that between 1995 and 1999, seismicity rates were higher during the winter than during the summer. This was further studied by Bettinelli et al. (2008) which showed that these seismicity rates variations correlated with hydrology-driven elastic Coulomb stress rate variations. However, Ader & Avouac (2013) analysed an extended catalogue, spanning 1995 to 2008, using the Schuster Spectrum and suggested that the previously observed annual periodicity could have been an artefact of incomplete declustering. This observation was nuanced by the dual result of annual signals when considering only large magnitude events, and when considering the International Seismological Centre's longer instrumental catalogue. These results did however come through the analysis of a relatively small number of earthquakes (16 earthquakes $M_L \ge 5.5$ and 210 earthquakes $M_w \ge 4$, respectively), and seemingly indicated a puzzling lack of observed periodicity at lower magnitudes where more earthquakes were available to study. Still, localized transient seismicity episodes are well documented in Nepal, with the occurrence of synchronous yet distant swarms across Nepal potentially indicating region-wide seismicity responsivity (Hoste-Colomer et al. 2017, 2018), whilst other swarms' activation is seemingly correlated with seasonal phenomena (Adhikari et al. 2021).

In this study, we re-explore the relationship between hydrological surface loads variations inferred from the Gravity Recovery and Climate Experiment (GRACE) and Gravity Recovery and Climate Experiment Follow-On (GRACE-FO) and the seismicity of Nepal prior to the $M_w$ 7.9 Gorkha earthquake over the 1995 to 2014 time period. We start by computing stress variations at depth induced by the hydrological cycle, using GRACE data and an elastic half-space model. We then calculate annual seismicity rate variations in eastern and central Nepal, and apply two types of declustering procedures to evaluate the influence of declustering on seismicity rates and their relation to magnitude thresholds. Finally, we investigate potential annual b-value variations, assess their significance and relate them to seasonal stress variations mainly driven by hydrological surface loads.

# 2 METHODS

## 2.1 Seismicity catalogue

The earthquake catalogue studied herein is a product of the operational efforts of the Nepalese Earthquake Monitoring and Research Center (NEMRC) and its ongoing collaboration with the Laboratoire de Détection et de Géophysique (DASE-CEA).

The Nepalese Seismological Network consists of 21 ZM500 short period (1s) vertical component seismometers. The bulk of these were deployed in 1994 for the eastern and central subnetwork, and between 1994 and 1998 for the western subnetwork (Pandey et al. 1999). The seismic network deployment occurred primarily in 1994, with no major modification thereafter and only local and brief station failures (Fig. S1). The specifics of the workflow used for the seismic bulleting generation can be found in Duverger et al. (2021).

The $M_w$ 7.9 2015 Gorkha earthquake resulted in heightened clustered seismic activity in central Nepal (Adhikari et al. 2023) with seismicity rates in 2020 remaining more than twice as high as in the preceding interseismic period. As such, we restrain the temporal extent of the catalogue from 1995 to 2014 to preserve the homogeneity of detection capabilities.

With the intent of building upon the works of Bollinger et al. (2007), Bettinelli et al. (2008) and Ader & Avouac (2013), we focus on the seismicity located along the Nepalese portion of the Main Himalayan Thrust (MHT; Fig. 1a). Indeed, most of the seismicity in Nepal occurs on the MHT, with focal mechanisms indicating shallow-dipping thrust faulting. Like previous studies, we select the seismicity based on the 3500 m altitude contour line which roughly correlates with the northern extension of Nepalese MHT seismicity. This is interpreted as a combined effect of vertical stress becoming the highest principal stress below the high chain, thus inhibiting underlying seismicity (Bollinger et al. 2004), and a rheological transition from a velocity-weakening to a velocity-strengthening regime around the 350 °C isotherm (Ader et al. 2012). Finally, we only retain earthquakes located above 40 km depth, a threshold much deeper than the MHT that corresponds to large catalogue depth uncertainties but ensures that selected earthquakes are almost entirely mid-crustal.

Due to the progress of the Nepalese Civil War, the western portion of the Nepalese seismic network became inoperable starting in 2004, significantly reducing the detection of western Nepalese seismicity. In addition to this, the Nepalese Seismic Network consists of two distinct subnetworks (western Nepal and eastern/central Nepal), requiring occasional merging and mixing of the bulletins created by each subnetwork in order to create a comprehensive national catalogue, an arduous and difficult process which can introduce duplicates or location biases. Accordingly, we focus only on central Nepal, between the 83°E and 87°E meridians, where the seismic record is most complete and use the catalogue built exclusively using the eastern/central Nepal subnetwork. To remove duplicate detections, we discard events relocated 5 s after another event, within a longitude distance of less than 0.1°, a latitude distance of less than 0.1° and at a magnitude difference of less than 0.5. This procedure removes 52 earthquake detections in our area of interest, most of which occur at times of large seismicity crises, as expected due to the increased waveform association complexity.

Whereas the subnetwork detected 34 299 earthquakes located within the borders of Nepal between 1995 and 2014, the final catalogue contains 6945 earthquakes of local magnitudes $M_L$ between 1 and 6, of which 1275 are earthquakes of magnitude greater or equal to 3.

To determine statistical variations within a seismicity catalogue, it is crucial to consider seismicity only above the completeness magnitude, $M_c$. This requires an accurate assessment of the completeness magnitude and its variations. For instance, high-frequency background seismic noise levels increase significantly during the summer monsoon due to the increase in surface hydrology, debris flows and their evacuation in rivers (Burtin et al. 2008), thereby resulting in a seasonally elevated $M_c$ (Fig. 2). Furthermore, the completeness magnitude changes over time with network degradations (Fig. S1) and workflow evolutions (Fig. 1b). To avoid bias from an underdetermined $M_c$, we adopt a conservative estimate of $M_c = 3$, based on the value reached by a b-value stability analysis (Woessner & Wiemer 2005) using a 30-d-long sliding window and considering only the day of the year of earthquake occurrence (Fig. 2).
<figure><img src="image_2.png" /></figure>

**Figure 1.** Spatiotemporal distribution of seismicity in Eastern and Central Nepal. (a) Spatial distribution of seismicity in relation to the main tectonic structures. Seismicity aligns itself along topography, the extent of MHT seismicity being roughly controlled to the north by the 3500 m contour line and to the south by the Main Boundary Thrust (MBT). Focal mechanisms are consistent with sub-horizontal thrust faulting. Seismometers of the NEMRC central network are represented by white triangles. Black circles depict surrounding seismicity ($M_L \ge 2.5$) relocated by the network between 1995 and 2014. Focal mechanisms are sourced from ISC bulletin up to 2023, Chen & Yang (2004), Burtin (2005) and Wang et al. (2017) and include the 2015 Gorkha earthquake and several of its aftershocks to underline the homogeneity of focal mechanisms in the area. (b) Temporal distribution of seismicity reveals several small seismic crises throughout the catalogue, with the most significant occurring in 1997, 2001 and 2005 (See year by year seismicity in fig. S2). Note that the lowest magnitudes events were no longer reported after early 2001, due to the smallest events measured only at short distances (less than 90 km) no longer being reported.

## 2.2 Seasonal rate variations

### 2.2.1 Declustering procedure

As alluded to by Schuster (1897) and demonstrated by Ader & Avouac (2013), temporal clustering of seismicity can lead to statistical artefacts that resemble periodicity. However, declustering is a non-trivial problem and several declustering approaches exist, each with their limitations. We decluster our catalogue using the growing cluster method of Reasenberg (1985) implemented in

ZMAP (Wiemer 2001) in order to approximate previous results obtained by Bollinger et al. (2007) and Ader & Avouac (2013). The parametrization details are provided in the supplementary information. We also use the Nearest Neighbour method described by Zaliapin et al. (2008) to separate clustered seismicity from 'background seismicity' with a constant nearest neighbour distance of $\eta = 10^{-4}$. This conservative threshold was chosen after inspecting the rescaled time/rescaled distance (Fig. S3). To limit biases due to large earthquakes outside our observation window, we apply these declustering methods before spatially cropping earthquakes,
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Completeness magnitude estimation and earthquake count as a function of day of year for the non-declustered catalogue. Completeness magnitude is evaluated through *b*-value stability analysis (Woessner & Wiemer 2005) with a 30-d-long sliding window to resolve seasonal variations and consider a representative number of earthquakes. The magnitude of completeness $M_c$ (solid white line) changes throughout the year due to background noise level changes, with a maximum of $M_c = 3$ during the summer monsoon. The computed maximum completeness magnitude is higher than the maximum curvature estimation of the completeness magnitude (solid black line).</figcaption></figure>

and include events from 1989—year of the first earthquakes considered from the previous local network—and 1994—deployment year of the current seismic network. Furthermore, to minimize biases linked to catalogue incompleteness on our declustering procedure, we systematically apply a magnitude threshold truncation before declustering, using a threshold magnitude higher or equal to the completeness magnitude. After Reasenberg declustering, 1080 earthquakes above $M_L = 3$ remain, while the Nearest-Neighbour declustering leaves 988 earthquakes above $M_L = 3$.

### 2.2.2 Annual periodicity

For each declustered catalogue, we determine annual variations of seismicity by calculating the *p*-value of the Schuster test (Schuster 1897; Shlien 1972):

$$p = \exp\left(-\frac{R^2}{N}\right), \qquad (1)$$

where *N* is the number of earthquakes considered, and *R* is the distance travelled by the phasor walkout at the considered period *T* such that:

$$R = \sqrt{\left[\sum_{n=1}^{N} \cos(2\pi t_n / T)\right]^2 + \left[\sum_{n=1}^{N} \sin(2\pi t_n / T)\right]^2} \qquad (2)$$

with $t_n$ the time of occurrence of the $n^{th}$ earthquake. The Schuster *p*-value describes the likelihood that a succession of punctual events result from a random sampling with uniform probability density. As

such, it is crucial for the catalogue to be properly declustered before calculating the Schuster *p*-value in order for the assumption of independence of events to be respected. As seasonal phenomena exhibit natural variance, the recurrence interval of seasonal forcing is not usually exactly one year. However, this variance is relatively low compared to the length of the year, with variations of ±10 per cent. As such, a yearly Schuster test remains appropriate to probe annual periodicity, although it may slightly underestimate said periodicity. Keeping earthquakes occurring between 1995 January 1 and 2014 December 31 ensures that each day of the year is equally represented. For simplicity’s sake, earthquakes that occur on December 31 of leap years are ignored at this stage. For each declustered catalogue we also investigate the effect of lower magnitude limit selection by calculating the annual Schuster *p*-value of catalogues that have undergone truncation with different magnitude thresholds before declustering, as described in the previous section. We may therefore see the influence of each declustering procedure on the periodicity analysis at different magnitude thresholds.

## 2.3 Magnitude frequency distribution

Due to the self-similarity of earthquake processes, seismic moments are widely thought to follow a power law (e.g. Kagan 2010). Thus, earthquake magnitudes, which traditionally scale logarithmically with seismic moments (Hanks & Kanamori 1979), are described as distributed according to an exponential distribution of the GR law (Gutenberg & Richter 1944):

$$N = 10^{a-bM}, \qquad (3)$$

where *N* is the number of earthquakes of magnitude greater than or equal to *M*, and *a* and *b* are constants. *b*, often referred to as ‘*b*-value’, represents the ratio between the abundance of small and large earthquakes. The dependence of the *b*-value on static stress has long been recognized (Scholz 1968), as has its evolution during the seismic cycle (Main et al. 1989; Rivière et al. 2018). Recent studies have also shown the impact of transient stressing on *b*-value (Ide et al. 2016; Pétrélis et al. 2021; Colledge et al. 2023). In this context, we investigate the possibility that annual hydrological loading on the MHT fault interface causes *b*-value variations.

We use the binned maximum likelihood estimator (Guttorp & Hopkins 1986; Tinti & Mulargia 1986) to calculate *b*-values, following the best practices laid out by Marzocchi & Sandri (2003). *b*-value is estimated as such:

$$b_{\text{mle}} = \frac{1}{\log_e(10)\Delta M} \log_e \left( 1 + \frac{\Delta M}{\bar{M} - M_c} \right) \qquad (4)$$

with $\Delta M$ the magnitude bin size, $\bar{M}$ is the mean magnitude and $M_c$ is the completeness magnitude below which events are discarded. In addition, we apply the *b*-positive approach (van der Elst 2021) which estimates the *b*-value, for binned magnitude distributions as:

$$b_{\text{pos}} = \frac{1}{\log_e(10)\Delta M} \log_e \left( 1 + \frac{\Delta M}{\bar{M}' - M_c'} \right), \qquad (5)$$

where $M'$ are positive magnitude differences and $M_c'$ is the completeness magnitude difference. Typical declustering methods induce biases in *b*-value estimations (Mizrahi et al. 2021). In consequence, we estimate *b*-values of non-declustered catalogues solely.
<table><caption>**Table 1.** Main shocks leading to significant short-term aftershock incompleteness (STAI), as identified by our outlier *b*-value recognition procedure.</caption><thead><tr><td>Main shock date</td><td>M<sub>L</sub></td><td>Longitude</td><td>Latitude</td></tr></thead><tbody><tr><td>1997/01/31</td><td>5.8</td><td>85.34</td><td>28.07</td></tr><tr><td>2001/07/16</td><td>6.0</td><td>84.71</td><td>27.96</td></tr><tr><td>2003/11/22</td><td>5.1</td><td>83.91</td><td>28.44</td></tr><tr><td>2004/07/11</td><td>6.1</td><td>83.49</td><td>30.73</td></tr><tr><td>2007/10/05</td><td>4.8</td><td>84.45</td><td>28.20</td></tr><tr><td>2007/10/29</td><td>5.3</td><td>84.49</td><td>27.96</td></tr><tr><td>2011/08/15</td><td>5.0</td><td>86.27</td><td>27.44</td></tr></tbody></table>

### 2.3.1 Transient incompleteness

Two main sources of incompleteness can bias *b*-value evaluations: variations in seasonal completeness magnitude and short-term aftershock incompleteness (STAI). We assume the former is appropriately addressed as we consider a high completeness magnitude, and focus on addressing the latter. To identify and account for STAI, we apply a jackknife strategy, that is, removing years one by one and calculating a sliding *b*-value estimate as a function of day of year. STAI is magnitude-dependent and primarily affects smaller magnitude earthquakes. Therefore, if STAI impacts our selection, the mean magnitude of detected earthquakes will increase, resulting in an underestimated *b*-value, according to eq. (4). If the removal of a specific year leads to transient higher *b*-values, we infer that STAI is present in the removed year. Using this strategy, we identify STAI by calculating the interquartile range (IQR) of the *b*-values estimates obtained through the jackknife method for each Julian day, and consider segments that repeatedly exceed $Q3+1.5 \times \text{IQR}$ (Fig. S4). Although our *b*-value estimates are not independent in this procedure, considering sliding *b*-values estimates for each year independently is not possible due to the insufficient number of earthquakes. To ensure that we are identifying STAI, we cross-reference our catalogue for the largest magnitude events when our procedure indicates anomalous *b*-values. As such, we remove segments corresponding to the earthquakes in Table 1 and the 10 d following these events, based on the typical length of STAI within our catalogue.

Whether this procedure is magnitude-independent is ambiguous. While we specifically remove the largest magnitude earthquakes and many large aftershocks, removing all events within a time window also excludes many smaller aftershocks. A standard ETAS synthetic analysis would not be helpful due to the imposed fixed *b*-value which would lead to no change of *b*-value by removal of segments with the largest magnitudes. The bias introduced by our STAI strategy can be assessed by examining some of the largest earthquakes with little to no STAI. Indeed, their removal from the stack in the yearly jackknife does not lead to a significant *b*-value drop. This points to the validity of our procedure to remove large STAI, assuming that *b*-values do not naturally vary too significantly.

### 2.3.2 *b*-value meaningfulness

For each catalogue subdivision, we test whether the magnitude frequency distribution is compatible with an exponential distribution by computing the Lilliefors test statistic (Lilliefors 1969). However, exponential distributions which have undergone binning are, by definition, no longer exponential distributions. More specifically, binned magnitudes extracted from a population that conforms to the GR law follow a geometric distribution, not an exponential one. The exponential Lilliefors test *p*-value statistic is therefore systematically underestimated when considering binned magnitude

distributions with $\Delta M = 0.1$. To mitigate the effect of binning, we recreate a pseudo-continuous distribution of magnitudes by adding uniformly distributed random noise in $[-\Delta M/2, \Delta M/2]$ to the magnitudes. Although this approach cannot recreate perfectly an exponential distribution, it leads to conservative yet reliable Lilliefors test statistics. This approach is popular when applying the exponential Lilliefors test to binned magnitudes (e.g. Herrmann & Marzocchi 2020; Collettini et al. 2022). We consider Lilliefors *p* values below 0.1 to indicate a departure from the exponential distribution, in which case *b*-values cannot be considered meaningful. To evaluate the significance of changes in *b*-value, we apply the bootstrap procedure proposed in Amorèse et al. (2010), which has been demonstrated to be a more reliable indicator than the traditionally used Akaike Information Criterion (AIC). We consider a threshold *p*-value of the two-sided bootstrap test of 0.05 for significance.

## 2.4 Load-induced stress variations

In the Himalayas, seasonal stress changes in the crust are significantly influenced by redistribution of water-mass, driven by the monsoon cycle. To model these stress variations, we use estimates of surface hydrological loads inferred from the monthly space and time varying GRACE and GRACE-FO Earth's gravity field (e.g. Tapley et al. 2004; Landerer et al. 2020). Specifically, we use the GRACE/-FO Multichannel Singular Spectrum Analysis (M-SSA) solution from Gauer et al. (2023). This solution combines multiple GRACE/-FO data sets and employs M-SSA to fill observational gaps and reduce noise. Although the solution is corrected for contributions of Glacial Isostatic Adjustment using ICE6-G-D model from Peltier et al. (2015, 2018), secondary contributions of atmospheric and non-tidal loading effects (Atmosphere and Ocean Dealiasing Level-1B model (AOD1B; Dobslaw et al. 2020)) are restored to compute the full load-induced stress changes. Data is provided as monthly Equivalent Water Height (EWH) on 1° by 1° grids with a spatial resolution of mass distributions of roughly 300 km (Figs 3a and b). For brevity, we herein refer to this data set as GRACE data.

To compute the stresses induced by GRACE-derived surface loading, we use the analytical Boussinesq model (Boussinesq 1885), following prior studies (e.g. Bettinelli et al. 2008). This model represents the Earth as an elastic, isotropic, homogeneous 3-D half-space on top of which loads are applied. It allows for the computation of resulting displacement components and both symmetric strain and stress tensors in response to surface loading [see Steer et al. (2014) and equations in supplementary material]. GNSS data in the study area shows large annual variations of both horizontal and vertical displacements, consistent with elastic deformation due to GRACE-derived seasonal hydrological loading (Bettinelli et al. 2008; Fu et al. 2013; Chanard et al. 2014), despite potential GNSS biases and other seasonal physical processes (Chanard et al. 2020). We therefore assume that hydrological loading dominates seasonal stress variations, compared to GNSS biases, thermoelastic and poroelastic deformation (Johnson et al. 2017).

Loads are considered only within a 5° radius of the calculation point, as most of the stress generated at shallow seismogenic depths results from local loads (Fig. S11). Since GRACE data is provided in spherical coordinates while the Boussinesq model uses a Cartesian geometry, we project GRACE loads onto local tangent plane coordinates around the coordinates where stresses are calculated. This is a reasonable first-order approximation given the spatial extent of the loads considered.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Annual surface loading and fault stressing variations in and around central Nepal between 2003 and 2022. (a) EWH annual variation amplitude across the Ganges watershed and surrounding region. The rectangle delimits the area in which stress is calculated to obtain average stress variations in central and eastern Nepal. (b) EWH annual phase variations in the same region with cyclic colour bar. Peak loading across the Ganges watershed occurs in phase with the end of the summer rainy season, that is, around September. (c) Stress variations on the MHT caused by mass redistribution across eastern Nepal assuming a friction coefficient of 0.1. Light lines in panel (c) represent one pixel contained in the black box of panels (a) and (b) and darker lines represent average values of the former. (d) Stress variations over one year.</figcaption></figure>

Furthermore, to accurately evaluate stresses at depth, we re-distribute the GRACE surface loads onto 0.1° by 0.1° grids (Figs. S12 and S13). We validated this approach by comparing it to results from a spherical harmonic layered-Earth surface loads model, yielding similar results (Fig. S14).

We then project the stress tensor ($\sigma$) onto the MHT, modelled with a 280° trend and 10° N dip, decomposing the resulting traction vector $T = \sigma \cdot n$ into its normal ($\sigma_n$), downdip shear ($\tau_d$) and alongstrike shear ($\tau_a$) components:

$$ \sigma_n = n \cdot T \qquad (6) $$

$$ \tau_d = (T - \sigma_n) \cdot t_d \qquad (7) $$

$$ \tau_a = (T - \sigma_n) \cdot t_a \qquad (8) $$

with $n$, $t_d$ and $t_a$ the normal, downdip and alongstrike vectors respectively. Assuming a friction coefficient $\mu = 0.1$ (Dal Zilio et al. 2019) and considering that tectonic stress accumulation occurs in the direction of convergence, we compute Coulomb failure stress $S$ variations as:

$$ \Delta S = \tau_d - \mu \sigma_n. \qquad (9) $$

A decrease in downdip shear stress direction results in an unloading of the MHT.

# 3 RESULTS

## 3.1 Annual stress variations in Nepal

South-east Asia experiences large annual mass redistribution, with EWH variations reaching up to several tens of cm. In Nepal, these variations can be as large as 20 cm (Fig. 3a). Driven primarily by the monsoon cycle, these annual mass redistributions are in phase across most of the Indian subcontinent (Fig. 3b). Hydrological masses usually reach their peak around September, quasi-synchronously across Nepal, following an increase starting in early summer and before a gradual withdrawal in the following months. These seasonal mass redistributions induce correlated normal and shear stress oscillations on the MHT (Fig. 3c). Normal stress variations exhibit the largest amplitude, with peak-to-peak amplitude of roughly 4 kPa across eastern and central Nepal. Shear stress variations are conversely of much smaller amplitude due to the low (10°) dip angle of the MHT and the absence of large load gradients in our modelled surface loads, mainly due to the GRACE data coarse spatial resolution. This leads to Coulomb stress variations of the order of 400 Pa (Fig. 3d).

## 3.2 Seismicity rate variations

Declustering procedures yield similar seismicity rates estimates, with roughly one earthquake every 8 d above magnitude 3.0 (Fig. 4). Schuster $p$ values do not indicate statistical significance for an annual period, regardless of the declustering procedure or magnitude
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Yearly averaged seismicity rate of $M_L \ge 3.0$ with a 30-d sliding window. The nearest neighbour declustered catalogue uses a nearest neighbour distance threshold of $\eta = 10^{-4}$. The full catalogue contains 1275 earthquakes, the Reasenberg declustered 1080 earthquakes, and the Nearest Neighbour declustered 988 earthquakes.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Schuster *p*-value at annual period for (top) Reasenberg declustering and (bottom) nearest neighbour declustering.</figcaption></figure>

range (Fig. 5). However, despite the overall high similarity at low threshold magnitudes, the Reasenberg declustering method produces slightly lower *p* values at higher magnitudes compared to the nearest neighbour method. This difference in declustering approaches might explain the seemingly significant periodic variations observed by Bollinger *et al.* (2007) and Ader *et al.* (2012), and suggests that the choice of declustering method can have a significant influence on the detection of earthquake seasonality.

## 3.3 Seasonal *b*-value variations

*b*-values start the year at roughly 0.9–1 and gradually increase to around 1.2 in July. *b*-values then decrease to their minimum of 0.9 by October before remaining relatively stable until the end of the year (Fig. 6 a-top). These heightened *b*-values during the summer

months cannot be explained by an underestimation of the completeness magnitude. Indeed, incompleteness above the determined completeness levels would result in artificially lowered *b*-value estimates. Nor can they either be caused by STAI as STAI was effectively mitigated, as evidenced by the absence of outlying *b*-values in our yearly jackknife analysis. We therefore interpret these variations as true seasonal variations in *b*-value. These variations are statistically significant at the 95 per cent confidence level when comparing the summer peak with most of the rest of the year (Fig. 7). To assess the quality of our *b*-value evaluation, we also perform a *b*-positive estimation (van der Elst 2021) using two approaches: one with magnitudes above the completeness magnitude of 3.0 and another with all magnitude differences exceeding the suggested threshold of 0.2. In both cases, we considered the entire catalogue, including STAI. We find that the *b*-positive results using the first of these approaches
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Yearly variations of *b*-value and Lilliefors *p*-values. Each light line corresponds to annual variations calculated after removal of a year, while bold lines correspond to annual variations considering all years. (a) *b*<sub>mlc</sub> with a sliding window of 200 earthquakes after removal of the seismic crises discussed in subsection 2.3 and for *M*<sub>L</sub> ≥ 3.0, and *b*<sub>pos</sub> with a sliding window of 100 positive magnitude differences using *M*′ ≥ 0 and *M*<sub>L</sub> ≥ 3.0, and keeping STAI. (b) Lilliefors *p*-values for magnitude distribution and positive magnitude difference distribution at each point in the corresponding plots of panel (a). Both panels have sliding windows spanning between roughly 40 and 60 d. Versions with precise window length and inclusion of seismic crises are available in supplementary figs S6 to S10.</figcaption></figure>

(Fig. 6 a-bottom) are consistent with *b*<sub>mlc</sub> estimates (Fig. 6 a-top). Conversely, the second approach yields magnitude difference distributions that fail to pass the Lilliefors test (Fig. S8). This makes the *b*-values estimated through using the 0.2 magnitude difference

threshold method unreliable. Exploring the combination of magnitude threshold and positive magnitude difference threshold reveals that a much larger magnitude difference threshold than 0.2 should be used for *b*-value estimates to be stable (Fig. S5). When using
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Significance of *b*-value variations over two moving windows of length 50 d using $M_L \ge 3$. The *b*-values around days 200 are significantly different from *b*-values during most of the rest of the year at a 95 per cent threshold. Low *p*-values of the bootstrap test around day 300 are not relevant and marked with an X. Indeed, the magnitude frequency distribution around that time of year does not conform to an exponential distribution (Fig. 6). The bootstrap *p*-values are determined using 5000 draws for each *b*-value comparison.</figcaption></figure>

the magnitude threshold of 3.0, we find Lilliefors *p*-values consistent with an exponential distribution of magnitudes for most of the year (Fig. 6 b-top) except around October. Likewise, using the same magnitude threshold, the positive magnitude difference distributions are similarly compatible with an exponential distribution (Fig. 6 b-bottom) except for a short amount of time around day 190 and, as with the magnitude distribution, around the month of October. These departures from the exponential distribution do not appear to be linked to specific years, as the results of our yearly jackknife indicate.

# 4 DISCUSSION

## 4.1 Transient stressing susceptibility of b-value

When comparing *b*-value variations to stress variations, we observe that the largest *b*-values do not correlate directly with hydrology-induced stress variations, but rather with stressing rate (Figs 8a–b). To ensure that these results are not affected by potential incompleteness in the subcatalogues, we follow the approach suggested by Marzocchi *et al.* (2020) by increasing the completeness magnitude. Increasing the threshold magnitude by 0.2 points to $M_c = 3.2$

leads to similar results, although the shape and amplitude of the seasonal variations are slightly altered (Figs 8c–d).

This behaviour is not the one usually observed in the literature, where *b*-value sensitivity to stress transients usually correlates with stress rather than stress rate (Ide *et al.* 2016; Tan *et al.* 2019; Colledge *et al.* 2023), with only limited evidence of stress rate control (e.g. Heki 2003). However, considering both hydrology and tectonic loading, and using reasonable estimates for secular stress increase in the seismogenic zone of the order of 1 kPa (Stevens & Avouac 2015), the summer peak in *b*-value also correlates with a short period of potential stress reversal, that is, negative stressing rate and stress reduction. That is, *b*-value increases when absolute stress levels are reduced at the peak of the hydrological load increase. This could be interpreted similarly to previous studies that reported a clear correlation with stress. However, these studies often focused on areas where the tectonic loading rate was relatively low compared to the periodic loading rate considered. In such cases, periodic stressing was likely the primary factor influencing *b*-value, while the importance of tectonic loading could be mostly ignored.

The observed stress reduction that occurs during stress reversal seemingly does not significantly prevent nucleation from occurring, as we do not detect significant reductions of seismicity rates.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** b-value as a function of (a,c) Coulomb stress variations and (b,d) Coulomb stress rate. b-value assessment is done with earthquakes above or equal to (a,b) $M_L$ 3.0, and (c,d) $M_L$ 3.2, comprised in a 30-d-long window, making each point independent. Error bars correspond to 1 standard deviation.</figcaption></figure>

However, considering hierarchical ruptures controlled by energy balances (Fukao & Furumoto 1985; Ide & Aochi 2005), punctual stress reduction may hinder seismic rupture propagation by reducing available elastic energy in the medium, thus making the propagation of earthquakes more likely to be arrested and increasing b-value, as suggested by Ide et al. (2016). Another potential explanation, related to the hierarchical model, involves the effect of the nucleation size. According to numerical investigation of a hierarchical model by Dublanchet (2020), b-value increases with normal stress due to reduced critical nucleation size which increases the number of partial ruptures occurring on each patch, which aligns with our observations. However, the effect of loading rate on nucleation size has also been investigated numerically (Kaneko et al. 2016) and experimentally (Guérin-Marthe et al. 2019), indicating that higher loading rates decrease the rupture nucleation size, thus the same principles could cause a theoretical decrease in b-value at the peak of the summer monsoon. These variations in nucleation lengths would however cause seismicity rate modulation, which we do not observe. We therefore favour the idea that the main control on the seismicity is fracture energy, rather than nucleation size or loading rate variations. Given the relatively small magnitude of

these stress variations, a large earthquake could easily offset the stress decrease induced by hydrology through static stress changes. Consequently, aftershocks sequences would be expected to have season-independent b-values. This may help explain the anomalous b-values observed during summer seismic crises in this study.

## 4.2 Seismicity rate variations and state of stress

Using the framework of rate-and-state (Dieterich 1994; Heimisson & Segall 2018; Heimisson & Avouac 2020), assuming that the effective normal stress $\sigma_0$ is much larger than the amplitude of the effective normal stress variations, the theoretical seismicity rate reads:

$$R(t) \approx r \frac{\exp\left(\frac{S(t)}{A\sigma_0}\right)}{1 + \frac{1}{t_a} \int_0^t \exp\left(\frac{S(t')}{A\sigma_0}\right) dt'} \quad (10)$$

where A is the rate-and-state velocity constant, and $t_a$ is the earthquake nucleation time equal to $\frac{A\sigma_0}{\dot{t}}$, with $\dot{t}$ the background shearing rate.
Using the stressing history along the MHT, we can make predic-
tions as to the impact of seasonal stressing on seismicity rates. These
predictions are however closely tied to the value of the key param-
eter Aσ₀. The seemingly significant periodicity found by Bettinelli
et al. (2008) led the authors to conclude that Aσ₀ held a value of
roughly 5 kPa. However, as pointed out by the authors themselves,
post-seismic creep has repeatedly pointed to Aσ₀ values of the or-
der of the MPa, with recent modelling consistently agreeing with
such values (e.g. Wang & Fialko 2018; Ingleby et al. 2020). Such
estimates would be consistent with a lithostatic pressure gradient
and hydrostatic pore pressure gradient, and would correspond to
seismicity rate variations of a fraction of a percent. These Aσ₀ val-
ues are, however, usually determined after large earthquakes such
as the 2015 M<sub>w</sub> 7.9 Gorkha earthquake when increased fracture per-
meability is likely to cause any existing pore pressure to be greatly
dissipated (Sibson 1992), thus leading to increased effective normal
stress. The value of Aσ₀ is therefore likely to be lower than these
estimates at the time of nucleation of earthquakes in the interseismic
period.

Regardless of the considered Aσ₀ and using stress build-up rates of 1 to 3 kPa.yr⁻¹ from Stevens & Avouac (2015), we find that the nucleation timescale is of the order of several years, indicating that the response of seismicity should be that of the high-frequency response regime (Beeler & Lockner 2003) and correlate with stress and not stressing rate. This result is opposed to that obtained by Bettinelli et al. (2008), and potentially explains the absence of a clear response of seismicity to periodic stress in our study as the response to the high-frequency regime is more damped than the low-frequency Coulomb regime. When considering the lowest values of Aσ₀, we find that stress oscillations would lead to seismicity rate variations of up to ±15 per cent with the peak seismicity rate occurring in spring. The clear absence of any such signal in background seismicity rates indicates that the amplitude of the seismicity rate modulation is much weaker than this (Fig. 4). Indeed, with our catalogue containing 900 earthquakes over a period of 20 yr, a 15 per cent seismicity rate modulation amplitude should result in a Schuster test value of less than 0.05 (Ader & Avouac 2013), while our results clearly do not support this (Fig. 5). This supports our conclusion that Aσ₀ is actually much larger than 5 kPa, with seismicity rate modulation amplitude being much lower than 15 per cent.

Furthermore, the presence of seismicity transients associated
with fluid movement indicate that stress along the MHT could at
least in part be controlled by fault-valve mechanisms with transient
up-flow triggering seismicity (Sibson 1992). Whether these tran-
sients are controlled by seasonal phenomena is beyond the scope
of this study but is crucial for understanding the effect of periodic
loading on seismogenic zones.

5 CONCLUSIONS

Re-investigating the claim that earthquake productivity in Nepal is
modulated by seasonal phenomena, we analyse 20 yr of earthquake
catalogue in eastern and central Nepal prior to the 2015 M<sub>w</sub> 7.9
Gorkha earthquake. While we cannot formally exclude the possi-
bility of a modulation of background seismicity rates, we are unable
to find a statistically significant annual signal regardless of consid-
ered magnitude thresholds and declustering procedure. Periodicity
significance test results did however vary moderately with declus-
tering procedure at high magnitude thresholds. As the declustering
procedure favoured by previous studies leads to the least conser-
vative results, we believe that declustering procedure choice could

have played a part in the results obtained by previous studies. How-
ever, we demonstrate the existence of annual b-value variations
during the same time period, with higher b-values being observed
in the summer. We theorize that these variations are linked to the
hydrological loading and unloading of seismogenic structures, with
increased b-value occurring at times of Coulomb stress decrease
along the Main Himalayan Thrust when considering background
tectonic loading. This is consistent with a fracture energy control
of seismic rupture, with a reduction of available fracture energy
reducing the growth potential of each rupture. This also indicates
potential variations in seismic hazard with total moment release
being lower at the height of the summer monsoon.

ACKNOWLEDGMENTS

We thank Jean-Philippe Avouac, Satoshi Ide and Éric Beaucé for their comments, Louis-Marie Gauer for providing GRACE data and the Department of Mines and Geology, in particular the team at the National Earthquake Monitoring and Research Center (NEMRC) for their work in data collection and curation. This work was supported by the LRC Yves Rocard (Laboratoire de Recherche Conventionné CEA-ENS-CNRS).

SUPPORTING INFORMATION

Supplementary data are available at *GJIRAS* online.

**Table S1.** Values used for Reasenberg (1985) declustering with Zmap (Wiemer, 2001). Parameters correspond to values mentioned in Bollinger et al. (2007) and in part used in Ader & Avouac (2013). (1) In both studies, it is our interpretation that the so-called maximum distance ‘D’ (expressed in kilometres) is in fact the interaction zone size scaling factor ‘Q’ as described in Reasenberg (1985). As such it is a dimensionless number.

**Figure S1.** Availability of network seismic stations between 1995 January 1 and 2015 January 1. Orange lines indicate periods where data exists while white lines denote a lack of data. Stations are presented East to West and are restricted to east and central Nepal stations. The apparent lack of availability in 1998 is an artefact due to the loss of the corresponding data in the arrival time database. No particular problems were, however, occurring on the network at the time and earthquakes were adequately recorded.

**Figure S2.** Seismicity rate for M<sub>L</sub> ≥ 3 over 60-d windows for each year. Seismic crises are easily identifiable. The sliding window is centred around the considered day, as such early 1995 and late 2014 data should be interpreted with this in mind.

**Figure S3.** Nearest neighbour distance density [as defined in Zaliapin & Ben-Zion (2013)]. The high density of distances around a roughly constant value of the rescaled distance η and the low density of a second mode of constant rescaled space component are the signature of low clustering.

**Figure S4.** Transient incompleteness identification using interquar-
tile range of jackknife b-value at each day of the year. When tran-
sient incompleteness is removed from the stack used to calculate
b-values, the b-value can significantly increase. We use a thresh-
old magnitude of 3.2, and a 10-d-long sliding window to group
successive significant deviations from the median.

**Figure S5.** b-positive as a function of threshold magnitude and magnitude-difference threshold. A catalogue with well-defined b-value should have stable b-positive above completeness magnitude and above completeness magnitude difference. The suggested magnitude difference threshold value in van der Elst (2021) of 0.2 is not
sufficient here, due to the large portion of the catalogue below com-
pleteness. First, removing earthquakes below magnitude M$_L$ = 2.5
leads to more stable estimates.

**Figure S6.** Annual jackknife of bmle, window duration and Lilliefors *p*-value for *M*<sub>c</sub> = 3.2, with a window size of 50 earthquakes, before removal of seismic crises. Large deviations from the mean behaviour are observed when the segments are removed, indicating a large bias caused by their inclusion. This is the tell-tale sign of Short-Term Aftershock Incompleteness (STAI). Each grey and coloured line corresponds to annual variations calculated after removal of a year, while black lines correspond to annual variations considering all years.

**Figure S7.** Annual jackknife of bmle, window duration and Lilliefors *p*-value for *M*<sub>c</sub> = 3.2, with a window size of 200 earthquakes after removal of seismic crises. Increasing the threshold magnitude does not significantly alter the result, but does cause some smoothing of *b*-value variations. Each grey and coloured line corresponds to annual variations calculated after removal of a year, while black lines correspond to annual variations considering all years.

**Figure S8.** Annual jackknife of $b_{\text{pos}}$, window duration and Lilliefors $p$-value for $M' \ge 0.2 M_c = 0$, with a window size of 200 magnitude differences, before removal of seismic crises. The Lilliefors $p$-value of the positive-magnitude-difference frequency distribution is very unstable, and its regular crossing of the 90 per cent threshold reflects departures from exponentiability which render the interpretation of $b$-value variations meaningless. Each grey and coloured line corresponds to annual variations calculated after removal of a year, while black lines correspond to annual variations considering all years.

**Figure S9.** Annual jackknife of $b_{\text{mle}}$, window duration and Lilliefors $p$-value for $M_c = 3.0$, with a window size of 200 earthquakes, after removal of seismic crises. This figure is an expanded version of the top plots of Figs 6(a) and 6(b), which includes the duration of the moving window. Each grey and coloured line corresponds to annual variations calculated after removal of a year, while black lines correspond to annual variations considering all years.

**Figure S10.** Annual jackknife of $b_{\text{pos}}$, window duration and Lilliefors $p$-value for $M' \ge 0.0 M_c = 3.0$, with a window size of 100 magnitude differences, before removal of seismic crises. This figure is an expanded version of the bottom plots of Figs 6(a) and 6(b), which includes the duration of the moving window. Each grey and coloured line corresponds to annual variations calculated after removal of a year, while black lines correspond to annual variations considering all years.

**Figure S11.** Influence of surface loads extent (*d*) on modelled stress variations in central Nepal. At seismogenic depths, that is, 20 km, stress variations are mainly controlled by the surface loads located in the immediate vicinity.

**Figure S12.** Influence of the spacing of the loads on stresses. The redistribution of loads is essential to evaluate stresses at shallow depths correctly. When no redistribution is applied ($\alpha_{\text{load}} = 0.1^\circ$) computed stresses are highly location-dependent and grossly incorrect.

**Figure S13.** Stress calculated with the Boussinesq model for September 2003 at coordinates 84°E, 28°N as a function of depth for different surface loading spacings ($\alpha_{\text{load}}$) and smoothing ($\alpha$). The distance between surface loads must be at most half of the depth at which stresses are considered. For instance, a surface source separation of 0.05° corresponds to intersource distances of roughly 5 km and are needed to resolve to accurately compute stress at 10 km depth.

**Figure S14.** Comparison of spherical-harmonics based layered elastic model and Boussinesq elastic homogeneous half-space model.

DATA AVAILABILITY

GRACE/−FO data is available at https://doi.org/10.18715/IPG P.2023.Igquie56. The GIA model used to correct GRACE/−FO gravity fields is available at https://podaac.jpl.nasa.gov/data set/Te llus_GIA_L3_1-DEG_V1.0. Atmospherical and non-tidal loading GRACE/−FO correction models are available at https://poda ac.jpl.nasa.gov/data set/GRACE AOD1B GRAV_FGZ_RL06. The seismicity catalogue and its declustered versions at and above completeness magnitude in the study area, as well as the hydro logically derived stress estimates, and all results derived from these and used in this study are available on Zenodo.org at https://doi.org/10.5281/zenodo.15649006 with a Creative Com mons Attribution 4.0 International Licence. All computation and figures were obtained using PYTHON, PYGMT (Wessel et al. 2019; Tian et al. 2025) and ZMAP (Wiemer 2001).

REFERENCES

Ader, T.J. & Avouac, J.-P., 2013. Detecting periodicities and declustering in earthquake catalogs using the Schuster spectrum, application to Hi- malayan seismicity, *Earth planet. Sci. Lett.*, **377–378**, 97–105.

Ader, T.J. et al., 2012. Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan thrust: implica tions for seismic hazard, *J. geophys. Res.: Solid Earth*, **117**(B4). doi:10.1029/2011JB009071.

Adhikari, L.B. et al., 2021. Orogenic collapse and stress adjustments re-vealed by an intense seismic swarm following the 2015 Gorkha Earth-quake in Nepal, *Front. Earth Sci.*, 9. doi:10.3389/feart.2021.659937.

Adhikari, L.B. et al., 2023. Seismically active structures of the Main Hi- malayan Thrust revealed before, during and after the 2015 Mw 7.9 Gorkha earthquake in Nepal, *Geophys. J. Int.*, **232**(1), 451–471.

Amorèse, D., Grasso, J.-R. & Rydelek, P.A., 2010. On varying b-values with depth: results from computer-intensive tests for Southern California, *Geophys. J. Int.*, **180**(1), 347–360.

Beeler, N.M. & Lockner, D.A., 2003. Why earthquakes correlate weakly with the solid Earth tides: Effects of periodic stress on the rate and prob- ability of earthquake occurrence, *J. geophys. Res.: Solid Earth*, **108**(B8). doi:10.1029/2001JB001518.

Bettinelli, P., Avouac, J.-P., Flouzat, M., Bollinger, L., Ramillien, G., Rajaure, S. & Sapkota, S., 2008. Seasonal variations of seismicity and geodetic strain in the Himalaya induced by surface hydrology, *Earth planet. Sci. Lett.*, **266**(3), 332–344.

Bollinger, L., Avouac, J.-P., Cattin, R. & Pandey, M.R., 2004. Stress buildup in the Himalaya, *J. geophys. Res.: Solid Earth*, **109**(B11). doi:10.1029/2003JB002911.

Bollinger, L., Perrier, F., Avouac, J.-P., Sapkota, S., Gautam, U. & Tiwari, D.R., 2007. Seasonal modulation of seismicity in the Himalaya of Nepal, *Geophys. Res. Lett.*, **34**(8), 2006GL029192. doi:10.1029/2006GL029192.

Boussinesq, J., 1885. Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques: principalement au calcul des déformations et des pressions que produisent, dans ces solides, des efforts quelconques exercés sur une petite partie de leur surface ou de leur intérieur; mémoire suivi de notes étendues sur divers points de physique mathématique et d'analyse, Gauthier-Villars.

Burtin, A., 2005. Seismotectonics of the Himalayan Arc from Regional Seis- mogram Moment Tensor Inversion, PhD thesis, Oregon State University.

Burtin, A., Bollinger, L., Vergne, J., Cattin, R. & Nábělek, J., 2008. Spectral analysis of seismic noise induced by rivers: A new tool to monitor spatiotemporal changes in stream hydrodynamics, *J. geophys. Res.: Solid Earth*, **113**(B5). doi:10.1029/2007JB005034.
Chanard, K., Avouac, J.-P., Ramillien, G. & Genrich, J., 2014. Modeling deformation induced by seasonal variations of continental water in the Himalaya region: sensitivity to Earth elastic structure, *J. geophys. Res.: Solid Earth*, **119**(6), 5097–5113.

Chanard, K., Métois, M., Rebischung, P. & Avouac, J.-P., 2020. A warning against over-interpretation of seasonal signals measured by the Global Navigation Satellite System, *Nat. Commun.*, **11**(1), 1375. doi:10.1038/s41467-020-15100-7.

Chen, W.-P. & Yang, Z., 2004. Earthquakes beneath the Himalayas and Tibet: evidence for strong lithospheric mantle, *Science*, **304**(5679), 1949–1952.

Cochran, E.S., Vidale, J.E. & Tanaka, S., 2004. Earth tides can trigger shallow thrust fault earthquakes, *Science*, **306**(5699), 1164–1166.

Colledge, M., Aubry, J., Chanard, K., Pétrélis, F., Duverger, C., Bollinger, L. & Schubnel, A., 2023. Susceptibility of Microseismic Triggering to Small Sinusoidal Stress Perturbations at the Laboratory Scale, *J. geophys. Res.: Solid Earth*, **128**(4), e2022JB025583. doi:10.1029/2022JB025583.

Collettini, C., Barchi, M.R., De Paola, N., Trippetta, F. & Tinti, E., 2022. Rock and fault rheology explain differences between on fault and distributed seismicity, *Nat. Commun.*, **13**(1), 5627. doi:10.1038/s41467-022-33373-y.

Dal Zilio, L., van Dinther, Y., Gerya, T. & Avouac, J.-P., 2019. Bimodal seismicity in the Himalaya controlled by fault friction and geometry, *Nat. Commun.*, **10**(1), 48. doi:10.1038/s41467-018-07874-8.

Davison, C., 1938. *Studies on the Periodicity of Earthquakes*, T. Murby & Company.

Dieterich, J.H., 1994. A constitutive law for rate of earthquake production and its application to earthquake clustering, *J. geophys. Res.: Solid Earth*, **99**(B2), 2601–2618.

Dobslaw, H., Dill, R., Bagge, M., Klemann, V., Boergens, E., Thomas, M., Dahle, C. & Flechtner, F., 2020. Gravitationally Consistent Mean Barystatic Sea Level Rise From Leakage-Corrected Monthly GRACE Data, *J. geophys. Res.: Solid Earth*, **125**(11), e2020JB020923. doi:10.1029/2020JB020923.

Dublanchet, P., 2020. Stress-Dependent b Value Variations in a Heterogeneous Rate-and-State Fault Model, *Geophys. Res. Lett.*, **47**(13), e2020GL087434. doi:10.1029/2020GL087434.

Dutilleul, P., Johnson, C.W., Bürgmann, R., Wan, Y. & Shen, Z.-K., 2015. Multifrequential periodogram analysis of earthquake occurrence: An alternative approach to the Schuster spectrum, with two examples in central California, *J. geophys. Res.: Solid Earth*, **120**(12), 8494–8515.

Duverger, C., Mazet-Roux, G., Bollinger, L., Trilla, A.G., Vallage, A., Hernandez, B. & Cansi, Y., 2021. A decade of seismicity in metropolitan France (2010–2019): The CEA/LDG methodologies and observations, *BSGF - Earth Sci. Bull.*, **192**, 25. doi:0.1051/bsgf/2021014.

Fu, Y., Argus, D.F., Freymueller, J.T. & Heflin, M.B., 2013. Horizontal motion in elastic response to seasonal loading of rain water in the Amazon Basin and monsoon water in Southeast Asia observed by GPS and inferred from GRACE, *Geophys. Res. Lett.*, **40**(23), 6048–6053.

Fukao, Y. & Furumoto, M., 1985. Hierarchy in earthquake size distribution, *Phys. Earth planet. Inter.*, **37**(2), 149–168.

Gauer, L.-M., Chanard, K. & Fleitout, L., 2023. Data-driven gap filling and spatio-temporal filtering of the grace and grace-fo records, *J. geophys. Res.: Solid Earth*, **128**(5), e2022jb025561. doi:10.1029/2022JB025561.

Guérin-Marthe, S., Nielsen, S.B., Bird, R., Giani, S. & Di Toro, G., 2019. Earthquake nucleation size: evidence of loading rate dependence in Laboratory Faults, *J. geophys. Res.: Solid Earth*, **124**(1), 689–708.

Gutenberg, B. & Richter, C.F., 1944. Frequency of earthquakes in California, *Bull. seism. Soc. Am.*, **34**(4), 185–188.

Guttorp, P. & Hopkins, D., 1986. On estimating varying b values, *Bull. seism. Soc. Am.*, **76**(3), 889–895.

Hanks, T.C. & Kanamori, H., 1979. A moment magnitude scale, *J. geophys. Res.: Solid Earth*, **84**(B5), 2348–2350.

Heimisson, E.R. & Avouac, J.-P., 2020. Analytical prediction of seismicity rate due to Tides and Other Oscillating Stresses, *Geophys. Res. Lett.*, **47**(23), e2020GL090827. doi:10.1029/2020GL090827.

Heimisson, E.R. & Segall, P., 2018. Constitutive law for earthquake production based on rate-and-state friction: Dieterich 1994 Revisited, *J. geophys. Res.: Solid Earth*, **123**(5), 4141–4156.

Heki, K., 2003. Snow load and seasonal variation of earthquake occurrence in Japan, *Earth planet. Sci. Lett.*, **207**(1), 159–164.

Herrmann, M. & Marzocchi, W., 2020. Inconsistencies and Lurking Pitfalls in the Magnitude–Frequency Distribution of High-Resolution Earthquake Catalogs, *Seism. Res. Lett.*, **92**(2A), 909–922.

Hoste-Colomer, R., Bollinger, L., Lyon-Caen, H., Burtin, A. & Adhikari, L.B., 2017. Lateral structure variations and transient swarm revealed by seismicity along the Main Himalayan Thrust north of Kathmandu, *Tectonophysics*, **714–715**, 107–116.

Hoste-Colomer, R. et al., 2018. Lateral variations of the midcrustal seismicity in western Nepal: Seismotectonic implications, *Earth planet. Sci. Lett.*, **504**, 115–125.

Hsu, Y.-J., Kao, H., Bürgmann, R., Lee, Y.-T., Huang, H.-H., Hsu, Y.-F., Wu, Y.-M. & Zhuang, J., 2021. Synchronized and asynchronous modulation of seismicity by hydrological loading: A case study in Taiwan, *Sci. Adv.*, **7**(16), eabf7282. doi:10.1126/sciadv.abf7282.

Ide, S. & Aochi, H., 2005. Earthquakes as multiscale dynamic ruptures with heterogeneous fracture surface energy, *J. geophys. Res.: Solid Earth*, **110**(B11). doi:10.1029/2004JB003591.

Ide, S., Yabe, S. & Tanaka, Y., 2016. Earthquake potential revealed by tidal influence on earthquake size-frequency statistics, *Nat. Geosci.*, **9**(11), 834–837.

Ingleby, T., Wright, T.J., Hooper, A., Craig, T.J. & Elliott, John.R., 2020. Constraints on the geometry and frictional properties of the Main Himalayan Thrust Using Coseismic, Postseismic, and interseismic deformation in Nepal, *J. geophys. Res.: Solid Earth*, **125**(2), e2019JB019201. doi:10.1029/2019JB019201.

Johnson, C.W., Fu, Y. & Bürgmann, R., 2017. Stress models of the annual hydro atmospheric, thermal, and tidal loading cycles on California faults: perturbation of background stress and changes in seismicity, *J. geophys. Res.: Solid Earth*, **122**(12), 10 605–10 625.

Johnson, C.W., Fu, Y. & Bürgmann, R., 2020. Hydro atmospheric modulation of stress and seismicity on shallow faults in southern Alaska, *Earth planet. Sci. Lett.*, **530**, 115904. doi:10.1016/j.epsl.2019.115904.

Kagan, Y.Y., 2010. Earthquake size distribution: power-law with exponent $B = 12?$, *Tectonophysics*, **490**(1), 103–114.

Kaneko, Y., Nielsen, S.B. & Carpenter, B.M., 2016. The onset of laboratory earthquakes explained by nucleating rupture on a rate-and-state fault, *J. geophys. Res.: Solid Earth*, **121**(8), 6071–6091.

Landerer, F.W. et al., 2020. Extending the Global Mass Change Data Record: GRACE follow-on instrument and science data performance, *Geophys. Res. Lett.*, **47**(12), e2020GL088306. doi:10.1029/2020GL088306.

Lilliefors, H.W., 1969. On the Kolmogorov-Smirnov Test for the exponential distribution with Mean Unknown, *J. Am. Stat. Assoc.*, **64**(325), 387–389.

Main, I.G., Meredith, P.G. & Jones, C., 1989. A reinterpretation of the precursory seismic b-value anomaly from fracture mechanics, *Geophys. J. Int.*, **96**(1), 131–138.

Marzocchi, W. & Sandri, L., 2003. A review and new insights on the estimation of the b-value and its uncertainty, *Ann. Geophys.*, **46**(6). doi:

Marzocchi, W., Spassiani, I., Stallone, A. & Taroni, M., 2020. How to be fooled searching for significant variations of the b-value, *Geophys. J. Int.*, **220**(3), 1845–1856.

Merian, P., 1834. *Über die in Basel wahrgenommenen Erdbeben: nebst einigen Untersuchungen über Erdbeben im Allgemeinen*, Wieland.

Métivier, L., de Viron, O., Conrad, C.P., Renault, S., Diament, M. & Patau, G., 2009. Evidence of earthquake triggering by the solid earth tides, *Earth planet. Sci. Lett.*, **278**(3-4), 370–375.

Mizrahi, L., Nandan, S. & Wiemer, S., 2021. The effect of declustering on the size distribution of Main shocks, *Seism. Res. Lett.*, **92**(4), 2333–2342.

Pandey, M.R., Tandukar, R.P., Avouac, J.-P., Vergne, J. & Héritier, T., 1999. Seismotectonics of the Nepal Himalaya from a local seismic network, *J. Asian Earth Sci.*, **17**(5), 703–712.

Pantobe, L., Burtin, A., Chanard, K. & Komorowski, J.-C., 2024. Evolution of shallow volcanic seismicity in the hydrothermal system of La Soufrière de Guadeloupe following the April 2018 Mlv 4.1 earthquake, *J. Volc. Geoth. Res.*, **447**, 107989. doi:10.1016/j.jvolgeores.2023.107989.
Peltier, W.R., Argus, D.F. & Drummond, R., 2015. Space geodesy constrains ice age terminal deglaciation: The global ICE-6G.C (VM5a) model, *J. geophys. Res.: Solid Earth*, **120**(1), 450–487.

Peltier, W.R., Argus, D.F. & Drummond, R., 2018. Comment on “An Assessment of the ICE-6G.C (VM5a) Glacial Isostatic Adjustment Model” by Purcell et al., *J. geophys. Res.: Solid Earth*, **123**(2), 2019–2028.

Perrey, A., 1861. Sur la fréquence des Tremblements de terre relativement à l'âge de la lune pendant la seconde moitié du dix-huitième siècle, et sur la fréquence deu phénomène relativement au passage de la lune au méridien, *C. R. Séances Acad. Sci.*, **LII**, 146–151.

Pétrélis, F., Chanard, K., Schubnel, A. & Hatano, T., 2021. Earthquake sensitivity to tides and seasons: Theoretical studies, *J. Stat. Mech. Theor. Exp.*, **2021**(2), 023404. doi:10.1088/1742-5468/abda29.

Reasenberg, P.A., 1985. Second-order moment of central California seismicity, 1969–1982, *J. geophys. Res.: Solid Earth*, **90**(B7), 5479–5495.

Rivière, J., Lv, Z., Johnson, P.A. & Marone, C., 2018. Evolution of b-value during the seismic cycle: Insights from laboratory experiments on simulated faults, *Earth planet. Sci. Lett.*, **482**, 407–413.

Scholz, C.H., 1968. Experimental study of the fracturing process in brittle rock, *J. geophys. Res. (1896-1977)*, **73**(4), 1447–1454.

Scholz, C.H., 2015. On the stress dependence of the earthquake b value, *Geophys. Res. Lett.*, **42**(5), 1399–1402.

Schuster, A., 1897. On lunar and solar periodicities of earthquakes, *Proc. Roy. Soc. Lond.*, **61**(369-377), 455–465.

Shlien, S., 1972. Earthquake-Tide Correlation, *Geophys. J. Int.*, **28**(1), 27–34.

Sibson, R.H., 1992. Implications of fault-valve behaviour for rupture nucleation and recurrence, *Tectonophysics*, **211**(1), 283–293.

Steer, P., Simoes, M., Cattin, R. & Shyu, J.B.H., 2014. Erosion influences the seismicity of active thrust faults, *Nat. Commun.*, **5**(1), 5564. doi:10.1038/ncomms6564.

Stevens, V.L. & Avouac, J.-P., 2015. Interseismic coupling on the main Himalayan thrust, *Geophys. Res. Lett.*, **42**(14), 5828–5837.

Tan, Y.J., Waldhauser, F., Tolstoy, M. & Wilcock, W. S.D., 2019. Axial Seamount: Periodic tidal loading reveals stress dependence of the earthquake size distribution (b value), *Earth planet. Sci. Lett.*, **512**, 39–45.

Tapley, B.D., Bettadpur, S., Watkins, M. & Reigber, C., 2004. The gravity recovery and climate experiment: Mission overview and early results, *Geophys. Res. Lett.*, **31**(9). doi:10.1029/2004GL01992.

Tian, D. et al., 2025. PyGMT: A python interface for the generic mapping tools, Zenodo. doi:10.5281/zenodo.15071586.

Tinti, S. & Mulargia, F., 1986. On the estimation of the geometric distribution parameter in seismology, *Statistica*, **46**(2), 163–188.

Ueda, T., Kato, A., Johnson, C.W. & Terakawa, T., 2024. Seasonal modulation of crustal seismicity in Northeastern Japan Driven by Snow Load, *J. geophys. Res.: Solid Earth*, **129**(3), e2023JB028217. doi:10.1029/2023JB028217.

van der Elst, N.J., 2021. B-positive: a robust estimator of aftershock magnitude distribution in transiently incomplete catalogs, *J. geophys. Res.: Solid Earth*, **126**(2), e2020JB021027. doi:10.1029/2020JB021027.

Wang, K. & Fialko, Y., 2018. Observations and modeling of coseismic and postseismic deformation due to the 2015 M 7.8 Gorkha (Nepal) Earthquake, *J. geophys. Res.: Solid Earth*, **123**(1), 761–779.

Wang, W., Shearer, P.M., Vidale, J.E., Xu, X., Trugman, D.T. & Fialko, Y., 2022. Tidal modulation of seismicity at the Coso geothermal field, *Earth planet. Sci. Lett.*, **579**, 117335. doi:10.1016/j.epsl.2021.117335.

Wang, X., Wei, S. & Wu, W., 2017. Double-ramp on the Main Himalayan Thrust revealed by broadband waveform modeling of the 2015 Gorkha earthquake sequence, *Earth planet. Sci. Lett.*, **473**, 83–93.

Wessel, P., Luis, J.F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W.H.F. & Tian, D., 2019. The generic mapping tools version 6, *Geochem. Geophys. Geosyst.*, **20**(11), 5556–5564.

Wiemer, S., 2001. A software package to analyze seismicity: ZMAP, *Seism. Res. Lett.*, **72**(3), 373–382.

Wilcock, W. S.D., 2009. Tidal triggering of earthquakes in the Northeast Pacific Ocean, *Geophys. J. Int.*, **179**(2), 1055–1070.

Woessner, J. & Wiemer, S., 2005. Assessing the quality of earthquake catalogues: estimating the magnitude of completeness and its uncertainty, *Bull. seism. Soc. Am.*, **95**(2), 684–698.

Zaliapin, I., Gabrielov, A., Keilis-Borok, V. & Wong, H., 2008. Clustering analysis of seismicity and aftershock identification, *Phys. Rev. Lett.*, **101**(1), 018501. doi:10.1103/PhysRevLett.101.018501.

Zhan, Z. & Shearer, P.M., 2015. Possible seasonality in large deep-focus earthquakes, *Geophys. Res. Lett.*, **42**(18), 7366–7373.