# Seismic Velocity Recovery in the Subsurface: Transient Damage and Groundwater Drainage Following the 2015 Gorkha Earthquake, Nepal

Luc Illien, Christoph Sens-Schönfelder, Christoff Andermann, Odin Marc, Kristen L. Cook, Lok B. Adhikari, Niels Hovius

► To cite this version:

Luc Illien, Christoph Sens-Schönfelder, Christoff Andermann, Odin Marc, Kristen L. Cook, et al. Seismic Velocity Recovery in the Subsurface: Transient Damage and Groundwater Drainage Following the 2015 Gorkha Earthquake, Nepal. Journal of Geophysical Research: Solid Earth, 2022, 127, 10.1029/2021JB023402. insu-03619990

HAL Id: insu-03619990

https://insu.hal.science/insu-03619990

Submitted on 25 Mar 2022

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
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2021JB023402

**Key Points:**

* We estimate a recovery time scale (<1 yr) in seismic velocity changes after the Gorkha earthquake using ambient noise correlations

* Velocity recoveries are modeled with relaxation functions characterized by a constant maximum relaxation timescale that is peak ground velocity-independent

* We highlight a transient enhanced permeability from the velocity changes in the first ~6 months following the main shock

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

L. Illien,
lillien@gfz-potsdam.de

**Citation:**

Illien, L., Sens-Schönfelder, C., Andermann, C., Marc, O., Cook, K. L., Adhikari, L. B., & Hovius, N. (2022). Seismic velocity recovery in the subsurface: Transient damage and groundwater drainage following the 2015 Gorkha earthquake, Nepal. *Journal of Geophysical Research: Solid Earth*, 127, e2021JB023402. https://doi.org/10.1029/2021JB023402

Received 14 OCT 2021
Accepted 13 JAN 2022

© 2022 The Authors.
This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes.

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

# Seismic Velocity Recovery in the Subsurface: Transient Damage and Groundwater Drainage Following the 2015 Gorkha Earthquake, Nepal

Luc Illien¹,² ID, Christoph Sens-Schönfelder¹ ID, Christoff Andermann¹ ID, Odin Marc³ ID, Kristen L. Cook¹ ID, Lok B. Adhikari⁴ ID, and Niels Hovius¹,² ID

¹Helmholtz Centre Potsdam, German Research Centre for Geosciences (GFZ), Potsdam, Germany, ²Department of Geosciences, University of Potsdam, Potsdam, Germany, ³Géosciences Environnement Toulouse, UMR, CNRS/IRD/UPS, Observatoire Midi-Pyrénées, Toulouse, France, ⁴Department of Mines and Geology, Kathmandu, Nepal

**Abstract** Shallow earthquakes frequently disturb the hydrological and mechanical state of the subsurface, with consequences for hazard and water management. Transient post-seismic hydrological behavior has been widely reported, suggesting that the recovery of material properties (relaxation) following ground shaking may impact groundwater fluctuations. However, the monitoring of seismic velocity variations associated with earthquake damage and hydrological variations are often done assuming that both effects are independent. In a field site prone to highly variable hydrological conditions, we disentangle the different forcing of the relative seismic velocity variations $\delta v$ retrieved from a small dense seismic array in Nepal in the aftermath of the 2015 $M_w$ 7.8 Gorkha earthquake. We successfully model transient damage effects by introducing a universal relaxation function that contains a unique maximum relaxation timescale for the main shock and the aftershocks, independent of the ground shaking levels. Next, we remove the modeled velocity from the raw data and test whether the corresponding residuals agree with a background hydrological behavior we inferred from a previously calibrated groundwater model. The fitting of the $\delta v$ data with this model is improved when we introduce transient hydrological properties in the phase immediately following the main shock. This transient behavior, interpreted as an enhanced permeability in the shallow subsurface, lasts for ~6 months and is shorter than the damage relaxation (~1 yr). Thus, we demonstrate the capability of seismic interferometry to deconvolve transient hydrological properties after earthquakes from non-linear mechanical recovery.

**Plain Language Summary** Earthquake ground shaking damage the rocks in the subsurface of the Earth, altering their strength and their permeability. After the main shock, the rock properties slowly return to their pre-earthquake state, but the duration of this recovery is poorly constrained. One way to investigate these time-dependent changes is through the monitoring of seismic velocity inferred from ambient ground vibration recorded at seismic stations. Here, we constrain the evolution of seismic velocity following the large 2015 $M_w$ 7.8 Gorkha earthquake in Nepal, in a field site characterized by seasonal groundwater fluctuations. We find that the velocity recoveries after the main shock and the aftershocks can be modeled with the same recovery timescale, independently from the initial shaking intensity. This suggests that earthquakes of different sizes activate the same geological structures and mechanisms during the recovery phase. Thanks to the unique hydrological setting of our field site and a model that links seismic velocity and groundwater level, we also show that this change of rock properties after the main shock is accompanied by a transient change in hydrological properties, an observation inferred for the first time with seismic measurement.

## 1. Introduction

Following the passage of seismic waves, a wide range of transient effects have been observed near the Earth's surface, including increased landslide rates (Marc et al., 2015), enhanced permeability (Manga et al., 2012; Xue et al., 2013), and perturbations of frictional properties in fault zones (Pei et al., 2019). These observations suggest that earthquakes induce a lingering effect in the properties of near-surface rocks that may be linked to non-linear mesoscopic elasticity (NLME, e.g., Gassenmeier et al., 2016; Marc et al., 2021). This phenomenon is generally expressed by a drop in elastic moduli after a dynamic or static strain perturbation that is followed by a non-instantaneous recovery of these moduli. This recovery phase, also called relaxation or slow dynamics, is linear on a logarithmic time scale (Snieder et al., 2017) and can last anywhere from a few seconds (Shokouhi et al., 2017)
to several years (Brenguier et al., 2008; Gassenmeier et al., 2016). Because most subsurface materials display this behavior (Gliozi et al., 2018; Shokouhi et al., 2017), understanding the amplitudes and timescales of the damage and recovery process of NLME is important for post-earthquake hazard mitigation.

In the field, the study of slow dynamics has been particularly advanced by the development of seismic interferometry techniques that monitor relative seismic velocity changes $\delta v = dv/v$ in the subsurface over time. Observations of co-seismic velocity drop and subsequent recovery in epicentral areas now abound and have been obtained from seismic ambient noise correlations (Brenguier et al., 2008; Gassenmeier et al., 2016; Hobiger et al., 2014; Wegler & Sens-Schönfelder, 2007) or waveform deconvolution in boreholes (Nakata & Snieder, 2011; Sawazaki et al., 2009; Wu et al., 2010). However, constraints on the physical mechanisms responsible for NLME in the field and the prediction of its amplitudes, timescales, and associated effects have remained scare for several reasons. First, the spatially averaged nature of the observation techniques does not allow for the precise identification of the responsible relaxation process among the many post-seismic processes acting at all depths and scales within a perturbed substrate. This complexity has prompted seismologists to use exponential functions characterized by variable timescales to fit velocity recoveries caused by individual events (Gassenmeier et al., 2016; Hobiger et al., 2014; Qin et al., 2020) rather than using particular physical relaxation models constrained from laboratory experiments (Bittner & Popovics, 2021; Lieou et al., 2017; Ostrovsky et al., 2019). Although this empirical approach can facilitate comparison between events, the understanding and prediction of the wide range of different recovery timescales (from minutes to years) between studies and sometimes within the same epicentral area (Viens et al., 2018) are limited. Moreover, aftershocks may induce superposed damage and healing processes, which may affect the observed recovery time of the main shock (Sawazaki et al., 2018).

The effects of slow dynamics may be obscured by hydrological fluctuations (Illien et al., 2021; Kim & Lee, 2019; Sens-Schönfelder & Wegler, 2006), which can influence the seismic velocity. Monitoring of hydrologically induced velocity variations ($\delta v_H$) is often done under the assumption that hydrological changes and NLME are independent processes that can be superimposed such that the observed $\delta v$ signal is simply the sum of hydrological and NLME effects ($\delta v = \delta v_{NLME} + \delta v_H$). However, there is evidence that both effects are not independent. It has been shown that the hydrological conditions of hillslopes can alter the NLME-response to dynamic strain (Bontemps et al., 2020). Moreover, as mentioned above, transient hydrological behavior following co-seismic ground shaking has been widely reported in borehole measurements (Elkhoury et al., 2006; Shi et al., 2015; Xue et al., 2013) and streamflow (Wang et al., 2004), suggesting that the hydrological system is also impacted by the transient variation of material properties. For example, the opening of cracks, which is often used to explain coseismic velocity decreases, can also introduce a change in substrate permeability (Elkhoury et al., 2006; Xue et al., 2013). Lastly, the similarity between the seismic velocity recovery timescale (~50 days, Taira et al., 2015) and the duration of the stream discharge increase (Wang & Manga, 2015) observed after the 2014 South Napa earthquake suggests a strong link between relaxation-induced velocity changes and transient hydrological properties. Because of the complexity of both processes and their coupling, it has not yet been possible to document the shaking induced perturbation of the hydrological system by means of seismic interferometry.

To investigate the shaking induced variations of a hydrological system with seismic interferometry, we use a seismic-hydrological dataset from the Nepal Himalayas that (a) features strong hydrological forcing, (b) includes the recovery phase of a large crustal earthquake, and (c) is described by a calibrated hydrological model that connects precipitation input to seismic velocity variations (Illien et al., 2021). Our approach involves accurate observations of seismic velocity changes, correcting the velocity changes for NLME effects due to the seismic activity and finally investigating the ability of the hydrological model to describe the residual velocity changes during different phases of the main shock recovery. Our field site is located in the epicentral area of the 2015 M$_w$ 7.8 Gorkha Earthquake (Figure 1a), in the Bhote Koshi catchment in Nepal about 60 km north east of Kathmandu in the steep ridge and valley topography of the lesser Himalayas. The region experienced strong ground shaking (Wei et al., 2018), widespread landsliding (Roback et al., 2018), and numerous aftershocks (Adhikari et al., 2015). Due to a distinct wet and dry season in which ~80% of the annual precipitation occurs during the Indian Summer Monsoon between ~May and ~October (Bookhagen & Burbank, 2010; Brunello et al., 2020), the hydrological conditions at this site are highly variable. This combination of pronounced and well-constrained hydrological and seismic forcing makes our field site a suitable location to study the interplay of seismic damage and hydrology.

The paper is organized as follows: we present the data and the seismic interferometry technique used to estimate velocity changes in Section 2. Section 3 shows the corresponding raw velocity changes observed after the Gorkha
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Map of the study area. (a) Black solid lines show the isolines for the Gorkha co-seismic slip (in cm) from the inverted solution of Elliott et al. (2016) using INSAR data. Yellow dots account for aftershocks of magnitude >4 (Adhikari et al., 2015). Red stars show the epicenter of the Gorkha earthquake and its main aftershock (12 May 2015). Green square is the Bhote Koshi observatory and blue star is the water gauge for measuring stage height of the Bhote Koshi river. (b) Close-up on the Bhote Koshi observatory. Red triangles show the site where seismic stations are deployed. BK stands for Bhote Koshi river.</figcaption></figure>

earthquake and its aftershocks. In Section 4, we present and discuss our models used to compute synthetic $\delta v$ values based on models for damage and hydrology. Section 4.1 is devoted to the damage-induced variations $\delta v_{N-LME}$ in which we introduce a new approach to describe the effects of the Gorkha mainshock and its aftershocks in a consistent model whereas Section 4.2 explores the residuals of the damage-corrected $\delta v$ time series using the hydrological model of Illien et al. (2021). This allows us to assess transient variations of the hydrological system in the Bhote Koshi catchment following the Gorkha event.

## 2. Data and Methods for Estimating Seismic Velocity Changes

Three broadband seismic stations (three-components Trillium compact 120 s) were installed on the 6 June 2015, 42 days after the Gorkha main shock near Chaku village (Figure 1b) and recorded until the 23 October 2018. The seismic stations were installed on a bedrock terrace at a distance of ~100 m from each other to achieve highly resolved temporal averaging at the same location. The metasedimentary rocks of the terrace are covered by a layer of regolith and colluvium. Because our stations were deployed after the Gorkha main shock, we also used data from the Gumba station (Figure 1b) of the Nepalese Seismological Center to confirm that our field site experienced a co-seismic velocity drop nearby and is in a recovery phase. This station has a single component and is located at 4.3 km from our field site and 1,700 m higher the Chaku terrace. For Gumba station, we evaluated data from 1 January 2014 to 9 December 2015. Daily precipitation were also measured from a network of precipitation gauges set up in the Bhote Koshi observatory (see the Data Availability Statement Section). We note that no major landslides occurred in the vicinity of our seismic stations (Marc et al., 2019) which imply that observed velocity changes are unlikely to be caused by redistributions of surface materials.

### 2.1. Estimation of Daily Relative Seismic Velocity Changes

We use seismic ambient noise to monitor variations of seismic velocity in the subsurface (Sens-Schönfelder & Brenguier, 2019). To reduce the impact of high amplitude signals in the noise correlation process, we use the following pre-processing scheme: the seismic traces are trimmed to 1 hr segments, downsampled to 50 Hz (only for Chaku stations), and detrended. We filter Chaku stations in the 4–8 Hz frequency range and data from Gumba station in the 2–4 Hz range due to limited seismic energy at higher frequencies. We normalize Chaku amplitudes to 1 in the Fourier spectrum (spectral whitening) and perform single station cross correlation (SC method,HBiger et al., 2014), using
$$C_{k_1, k_2}(t_i, \tau) = \int_{t_i - T/2}^{t_i + T/2} \text{sgn}[X_{k_1}(t')] \cdot \text{sgn}[X_{k_2}(t' + \tau)] dt', \quad (1)$$

where $t_i$ is the time of the trace and $\tau$ is the lapse time of the correlation. $T$ is the length of the correlated noise segments determines the temporal resolution of the $\delta v$ time series. The sgn function represents the 1-bit normalization of the signal in which we set positive amplitudes to 1 and negative amplitudes to -1. $k_m$ stands for the different components $m = Z, N, E$ with $k_1 \neq k_2$ for SC. Because Gumba has one component only, we compute the autocorrelation of the vertical component $C_{zz}(t_i, \tau)$. Correlation functions are calculated with a time step of 1 hr before averaging them every 24 hr to obtain daily correlation functions (DCFs). We store all the DCFs in a correlation matrix, as shown in Figure S1.

We use the stretching technique (Sens-Schönfelder & Wegler, 2006) to estimate relative velocity variations. After a spatially homogeneous relative velocity change $\delta v = dv/v$ in the medium, the time delay $\delta \tau = dt/\tau$ can be observed in the DCFs coda with $\delta v = -\delta \tau$ where $\tau$ is the correlation lapse time and $dt$ is the absolute time shift of a coherent phase with travel time $\tau$. Depending on the daily velocity changes, the DCFs ($C(t_i, \tau)$) are stretched or compressed when compared to a long term average reference $\xi(\tau)$. To avoid the effects of a possible degradation of a unique reference when averaged over the whole time period (Sens-Schönfelder et al., 2014), we use multiple references $\xi_r(\tau)$ at the Chaku site by computing monthly references $\xi_r(\tau)$ with an overlap of 15 days (we illustrate the use of different references in Figure S2). For each of these references, we calculate the correlation coefficients $R_r(t_i, \epsilon_j)$ between stretched versions of the reference and the DCF such that

$$R_r(t_i, \epsilon_j) = \int_{\tau_1}^{\tau_2} C(t_i, \tau) \xi_r(\tau * (1 + \epsilon_j)) d\tau \quad (2)$$

where $\tau$ is the traveltime of waves in the DCF and $\epsilon_j$ indicates a set of stretch-values that are tested in the time window set by $[\tau_1, \tau_2]$. We define the length of the time window as follows: we skip four signal periods $T$, where one period corresponds to the lowest frequency of the bandpass filter we previously applied (here $T = 0.25$ s), before computing the stretching on a duration of 12 periods (corresponding window indicated on Figure S1). Introducing $\tau_1$ is necessary to avoid the use of early arrivals that are prone to changes in noise sources characteristics. All $R_r(t_i, \epsilon_j)$ values are stored in a similarity matrix.

For each reference, a first daily velocity measurement $\delta v_r(t_i)$ can be done by reading the amount of stretching $\epsilon_j$ that yields the daily maximum $R_r(t_i, \epsilon_j)$ value. Combining the measurements done with the $N$ various references, we stack all similarity matrices $R_r(t_i, \epsilon_j)$ after correcting for any average shifting ($\overline{\delta v_r}$) due to the velocity differences between the references (full method described in Sens-Schönfelder et al., 2014) following the relation

$$R(t_i, \epsilon_j) = \sum_{r=1}^{N} R_r(t_i, \epsilon_j) - \text{shift}(\overline{\delta v_r}). \quad (3)$$

$R(t_i, \epsilon_j)$ describes the daily velocity variations obtained from one combination of sensor components $k_1, k_2$. We applied this method to the three possible combinations (ZN, ZE, EN) for each of the Chaku stations. We finally stack the resulting nine $R(t_i, \epsilon_j)$ matrices (three stations with three combinations) and pick the $\epsilon_j(t)$ with the maximum $R(t_i, \epsilon_j)$ again. The final daily $\delta v(t_i)$ at the Chaku site is equal to this specific $\epsilon_j(t)$. For Gumba station, we use only one reference as the use of multiple references does not improve the retrieved $\delta v$ values.

## 2.2. Local Aftershocks Catalog and Estimation of Associated Velocity Changes $\delta v_A$

Aftershocks recorded after the Gorkha earthquake may bias the recovery timescale estimated after the main shock by inducing further velocity drops and recoveries. However, due to potentially large hydrological fluctuations at the daily timescale of the interferometric processing, it may be challenging to dissociate the effect of cumulative aftershocks from hydrologically induced velocity variations. To address this issue, dedicated velocity change measurements following local aftershocks were conducted at a finer temporal resolution. Despite aftershock catalogs being available for the Gorkha earthquake (Adhikari et al., 2015; Baillard et al., 2017), their relevance for our field site remains limited as they lack information about the local shaking at the Chaku site. Therefore, to estimate the cumulative effects of shaking due to the aftershocks on the velocities, we build a catalog based on the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Evidences for NLME at Chaku. Black dashed lines in plots (a), (b), and (d) indicate the date of the Gorkha earthquake. (a) Black dots show the raw daily $\delta v$ measured at Chaku. Colored solid lines display results using the hydrological model of Illien et al. (2021) with different initial conditions. (b) Blue lines show local daily precipitation. Black stars stand for the PGV of aftershocks that caused a seismic velocity drop at the Chaku site. (c) Shows the same data as (a), but with the $\delta v$ of each year plotted on top of each other. The 2015 velocity is in red. (d) $\delta v$ variations estimated from Gumba station.</figcaption></figure>

daily peak ground velocity (PGV) recorded at Chaku. We first retain days with PGV greater than $1e^{-4} \text{ m s}^{-1}$. In the field, this value is approximately an order of magnitude lower than the minimum excitation required to induce a detectable change in rock properties as reported in the literature (Elkhoury et al., 2006; Wu et al., 2010). To exclude potential spurious peaks due to local artifacts, we check if the corresponding signals were also recorded at another temporary station (Hindi station on Figure 1b) located at ~3 km from our site. Using this procedure, we pick 82 potential aftershocks.

To test whether these events triggered NLME, we perform single station cross correlations of the ambient noise centered around the 82 events using the same method described in Section 2.1, but with a 10-min interval for the estimation of $\delta v$. We find that 18 events triggered a seismic velocity drop that was observable at this resolution. We quantify the co-seismic velocity drops by taking the difference between the median $\delta v$ value of the 12 hr preceding the aftershocks (no detectable velocity drops occurred during this time span) and the median value of the first hour succeeding the events.

# 3. Seismic Velocity Changes

## 3.1. Evidence for Non-Linear Recovery After Gorkha Earthquake

In Figure 2a, we report the daily relative seismic velocity changes estimated at the Chaku site and the daily precipitation totals recorded at nearby precipitation gauges in Figure 2b. The Chaku $\delta v$ time series exhibits a clear annual cyclicity exerted by the climatic forcing with a consistent drop of up to 8% in measured $\delta v$ values during
the monsoon season. Because of these significant hydrology-induced velocity changes and our dataset starting at 42 days after the Gorkha main shock (dashed lines in Figure 2), the recognition of any non-hydrological component in the $\delta v$ time series is strongly overprinted. Nevertheless, several arguments pinpoint the presence of NLME recovery in our time series.

First, we report the velocity changes observed at the Gumba seismic station (Figure 2d) as general evidence for NLME in the study area. A clear velocity drop of ~5% is observed at the date of the Gorkha earthquake. We attribute the noisy nature of the measurements to the lack of averaging in the velocity retrieval at this station for which only a single component is available. For this reason and because of the limited data coverage after the main shock, we do not attempt to characterize the recovery phase following the main shock at this station. Nonetheless, the clear co-seismic drop shows that ground shaking during the Gorkha event has caused damage in the Bhote Koshi catchment that is likely followed by a phase of recovery of subsurface material properties.

A second observation pointing to NLME behavior comes from a comparison of the annual cycles in $\delta v$ as shown in Figure 2c. In 2016–2018, the mean annual $\delta v$ cycle peaked to approximately +4% at the end of the pre-monsoon season in May. At the same time of the year in 2015, a clear offset from this value was observed with $\delta v$ as low as approximately –1%. Despite our precipitation dataset starting the 6 June 2015, it is unlikely that this offset is caused by climatic conditions. Indeed, with the 2015 monsoon being rather weak compared to precipitation totals of other monsoons seasons (Figure S3), a dryer season would cause the 2015 $\delta v$ data to be relatively higher than in the other years. This was not observed.

Finally, the last argument indicating NLME processes comes from hydrological modeling. We previously showed that the seismic velocity at Chaku reflects the groundwater content of the substrate in the vicinity of seismic instruments (Illien et al., 2021). This can be shown using the precipitation data recorded at our field site (Figure 2b). For comparison, we report this model in Figure 2a. We consider two different initial conditions for $\delta v$ in our model: one using the initial observed $\delta v$ (green line in Figure 2a) and another using the expected $\delta v$ value based on observations from years 2016 to 2018 at this time of the year (red line in Figure 2a). Both synthetics show good agreement with the velocities from April 2016 to the end of the time series—the period in which the model was calibrated, assuming that the NLME effect should be negligible in comparison to the hydrological influence on $\delta v$. However, velocities in 2015 are largely overestimated by the hydrological model with an offset of ~4% at the start of the time series. We note that this mismatch is progressively reduced at later times and converge toward the hydrological calibration. This supports a significant second control on $\delta v$ during the observation period, in addition to the pervasive hydrological influence. Considering this list of arguments, we conclude that a recovery behavior due to NLME likely occurred at Chaku.

## 3.2. Seismic Velocity Drop and Recovery Induced by Single Aftershocks

We observe small velocity drops that are particularly visible during the first dry season of the Chaku dataset (starting ~November 2015 in Figure 2a). We attribute these drops to further dynamic strain perturbations induced by aftershocks. Figure 2b shows the PGV measurements corresponding to the aftershock catalog we described in Section 2.2. The occurrence of the velocity drops in the daily $\delta v$ time series agrees with the timing of the reported ground shaking.

Observed velocity drops range from 0.25% to 1.5% and appear to have a linear relationship with PGV values (0.25–1.3 cm s⁻¹; Figure S4) although with a moderate scatter ($R^2$ = 0.62). For events occurring during dry periods, a clear slow dynamics behavior is observed with a distinct nonlinear recovery in the following hours after the initial drop (Figures 3a–3d). We highlight the characteristic log-linear behavior by averaging the data at a 30-min resolution and showing the first 100 hr in $\delta v$ after the velocity drops in a log-linear plot (Figures 3e–3h. The fit of a log-linear function of the form $\delta v = s \log(t) + C$, typical of the NLME functional form (TenCate et al., 2000), gives a satisfactory representation of the velocities. To avoid the possible larger hydrological modulation of $\delta v$ at late recovery times, we will model aftershock effects considering only an early time span.
<figure><img src="image_5.png" /><figcaption>Figure 3. Velocity recoveries following aftershocks. (a–d) Shows the velocity obtained at a 10-min resolution with the red dots indicating the first 100 hr after the events. (e–h) Shows the close-up of the results in the first 100 hr after the events in a log-linear plot. Results are averaged at a 30-min resolution. The red lines depict the fit of a log-linear slope on the first 24 hr of relaxation.</figcaption></figure>

## 4. Modeling $\delta v$: Derivation and Implications

In this section, we develop and use models to fit the seismic velocity changes presented in Section 3 and discuss their implications. The classic approach to decompose seismic velocity changes $\delta v$ is a linear superposition of forcing that can be written as

$$ \delta v = \delta v_{\text{NLME}} + \delta v_{\text{H}} \quad (4) $$

where $\delta v_{\text{NLME}}$ are the velocity changes due to NLME and $\delta v_{\text{H}}$ are the hydrologically induced velocity changes. $\delta v_{\text{NLME}}$ can be further decomposed into two components representing the relaxations due to Gorkha ($\delta v_{\text{G}}$), and its aftershocks ($\delta v_{\text{A}}$). To go beyond the linear description of expression 3, which does not account for transient post-seismic hydrological behavior, we propose a modeling approach based on two iterations: we first model the effect of NLME using conventional exponential functions. This approach is compared to the use of universal relaxation functions $R(t)$ which are calibrated for the first time on field data and are characterized by constant relaxation timescales, independent from ground shaking amplitude. To avoid a contamination by strong hydrological variations in the fitting, we calibrate the functions $R(t)$ using the initial 24 hr $\delta v$ dynamics following aftershocks events.

In a second step, we remove the inferred $\delta v_{\text{NLME}}$ component from the $\delta v$ time series to obtain residuals that represent the hydrological induced variations $\delta v_{\text{H}}$ (Section 4.2). We test whether $\delta v_{\text{H}}$ is not only influenced by precipitation but also by seismic damage. Because the meteorological effect on $\delta v_{\text{H}}$ is well constrained by the model of Illien et al. (2021), we introduce a transient drainage parameter in this model to estimate $\delta v_{\text{H}}^*$, which represents the seismically forced part of the hydrological component.

### 4.1. Post-Seismic Relaxations

We first apply the classic approach to model the recovery as an exponential recovery of the moduli and show that despite having numerous parameters for each event, the model performance is insufficient. Therefore, we propose a new strategy that uses a universal relaxation function and allows the description of all aftershocks and the mainshock with one consistent model, facilitating a correction of the time series for NLME effects.

#### 4.1.1. Modeling $\delta v_{\text{NLME}}$ With Exponential Functions

In seismic interferometry studies (Gassenmeier et al., 2016; Hobiger et al., 2014; Qiu et al., 2020), the nonlinear recovery in seismic wave velocity $\delta v$ is often fitted with the following function:

$$ \delta v(t) = \delta v_0 \exp \left[ \frac{-(t - t_0)}{\tau} \right] + C \quad (5) $$
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Fitting of the aftershocks recoveries. On each plot (a-d), green lines show the best fit of the exponential function for each observed recoveries (black lines). Red lines show the best fitting model using the relaxation function with a constant maximum relaxation timescale ($\tau_{max}$ = 155 days).</figcaption></figure>

where $t_0$ is the time of the earthquake occurrence, $\delta v_0$ is the initial co-seismic velocity drop at the temporal observation scale, $\tau$ is a characteristic time scale of recovery, and $C$ is a permanent drop.

To estimate the three empirical parameters of the exponential model, we use the velocity changes computed during the first 24 hr following the four aftershocks presented in Section 3.2 (Figure 3). In this time-span, a clear drop-recovery signal with no apparent hydrological-induced variations is observed (Figure 3). Assuming that $C = 0$ for the small excitations caused by the aftershocks, we fit expression 5 to the four $\delta v$ time series to obtain the characteristic timescales for the aftershocks $\tau_A$. The recovery time constants range from $\tau_A = 1.18$ to 3.03 days (Figure 4). To demonstrate the performance of this model, we build two synthetic time-series for the velocity variations induced by all aftershocks $\delta v_A$ using these two end-member values (Figure 5a, full method in Text S1 in Supporting Information S1).

After removing the synthetic $\delta v_A$ from the full $\delta v$ Chaku time-series, we fit the residuals with Equation 5 to obtain the recovery time constant $\tau_G$ for the Gorkha earthquake. We find a best fitting model with $\tau_G = 198$ days and a confidence interval of 80 days $< \tau_G < 1,208$ days that includes all model solutions with a variance ratio above 95% (Figure S5 in Supporting Information S1). The value used for the aftershocks correction ($\tau_A$ of 1.18 or 3.03 days) does not influence the inferred $\tau_G$. Synthetic time-series corresponding to the joint effect from the mainshock and the aftershocks are in Figures 6a and 6b together with the data residuals after correction for $\delta v_{NLME}$. The strongest differences are observed in the early part of the recovery depending on the characteristic timescale $\tau_G$ chosen for the main shock. Despite using the longest time scale for aftershock recovery of $\tau_A = 3.03$ days, the recoveries

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Synthetic seismic velocities induced by aftershocks $\delta v_A$. (a) Models built with the two end-member values $\tau_A$ measured with the exponential functions. (b) Models built using superposition of the relaxation functions of models $R^{155}$, $R^{846}$, and $R^{250}$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6. NLME models built with exponential recovery functions. (a) Each curve indicates synthetic recoveries characterized by different $\tau_G$ within the 95% confidence interval of the best fitting model ($\tau_G$ = 198 days). We superposed on this curve the recoveries associated with the synthetic $\delta v_A$ time series ($\tau_A$ = 3.03 days). (b) Corresponding residuals from the models shown in (a). Light gray line show the raw data. Dashed line indicates the zoomed window for the plot shown in (c).</figcaption></figure>

seem to not be fully corrected between ~November 2015 and ~June 2016 (Figure 6c). This suggests that longer timescales of relaxation after aftershocks should be introduced to fully correct for $\delta v_A$.

### 4.1.2. Modeling $\delta v_{NLME}$ With a Universal Relaxation Function

A characteristic of NLME is that the functional form of the recovery process is linear on a logarithmic time scale (Figure 3). A very convenient way to model this behavior is provided by the universal relaxation function used by Snieder et al. (2017). In this framework, the relative seismic velocity changes are described by a relaxation function $R(t)$:

$$ \delta v(t) = \delta v_{ss} + s R(t - t_0) \quad (6) $$

where $\delta v_{ss}$ is the steady state value of $\delta v(t)$ and $s$ is a scaling factor. $R(t)$ is the relaxation function that represents a multitude of processes with characteristic timescales. These timescales are distributed between a lower bound $\tau_{min}$ and a maximum relaxation time $\tau_{max}$. This theory leads to a superposition of these exponential processes that is given by

$$ R(t) = \int_{\tau_{min}}^{\tau_{max}} \frac{1}{\tau} e^{-(t-t_0)/\tau} d\tau. \quad (7) $$

A justification of Equation 6 based on the Arrhenius law is given by Snieder et al. (2017) but we recall a few important properties of the relaxation function: $R(t)$ exhibits a logarithmic behavior between the bounds $\tau_{min}$ and $\tau_{max}$ and its value at $t=0$ is finite and determined by $R(0) = \ln(\tau_{max}/\tau_{min})$. The prefactor $1/\tau$ increases the contribution of the processes with the shortest relaxation times, which leads to a uniform distribution of barrier energies according to Arrhenius law. Figure 7 illustrates the influence of the parameters $\tau_{min}$ and $\tau_{max}$.

Because of the multi-scale character of the universal relaxation function, we can describe the effects of the weak aftershock perturbations and the strong perturbation induced by the main shock with the same relaxation times $\tau$ in Equation 7. As we observe logarithmic recovery from the earliest measurement in Figure 3, we fix the parameter $\tau_{min}$ to 1 hr corresponding to the observation timescale. In the lab, minimum relaxation times down to $10^{-2}$ s have been reported (Shokouhi et al., 2017; but these smaller timescales $\tau_{min}$ would not affect the model fit).

We construct three models for the NLME with the relaxation function (5). First, the recovery phases of the four aftershocks with the clear recoveries shown in Figure 3 are fitted by adjusting a single $\tau_{max}$ to minimize the cumulative squared residuals. This consists in (a) numerically integrating Equation 7 and (b) fitting Equation 6 to $\delta v_A$ by adjusting the scaling $s$ for each aftershock. The red lines in Figure 4 show the obtained data fit. The best fit is
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Sensitivity of the relaxation function *R*(*t*). (a) The different colors account for the different τ<sub>max</sub> indicated in the figure. τ<sub>min</sub> is fixed to 1 hr. (b) The different colors account for the different τ<sub>min</sub> indicated in the figure. τ<sub>max</sub> is fixed to 400 days.</figcaption></figure>

found with τ<sub>max</sub> = 155 days (misfit curve in Figure S6a in Supporting Information S1). We will refer to this model as R<sup>155</sup> where the superscript stands for the fitted maximum relaxation timescale τ<sub>max</sub>.

The second value for τ<sub>max</sub> is inferred by fitting the complete long term δν data for the recovery of the main shock (Figure 8). τ<sub>max</sub> = 846 days is the best estimate in this case (misfit curve, Figure S6b in Supporting Information S1). Finally, we estimate a third timescale τ<sub>max</sub>, combining the two previous measurements by stacking the misfit curves (Figure S6c in Supporting Information S1) corresponding to the fit of the four aftershocks (R<sup>155</sup>) and the misfit curve from the fit of the entire time series (R<sup>846</sup>). This combined estimate yields τ<sub>max</sub> = 250 days as the value minimizing the combined misfit.

We compute three different NLME models (R<sup>155</sup>, R<sup>846</sup>, and R<sup>250</sup>) characterized by the different τ<sub>max</sub> values (Figure 8). For the 18 aftershocks, we use the measured velocity drop values (Figure S4 in Supporting Information S1) to compute the value *s* in Equation 5 that scales the relaxation function *R*(*t*) and stack the resulting functions (Figure 5b). After removing the aftershock perturbations δν<sub>A</sub> from the Chaku δν time series, we use Equation 5 again to adjust the scaling of *R*(*t*) for the main shock recovery δν<sub>G</sub>. The total NLME-induced δν<sub>NLME</sub> from this procedure with its obtained residuals are shown in Figure 8. We note that in comparison with the exponential approach, the aftershocks induce a larger and long-lasting perturbation of δν (Figure 5) which better

<figure><img src="image_10.png" /><figcaption>**Figure 8.** NLME models built with relaxation functions. (a) Each curve indicates synthetic recoveries characterized by different maximum relaxation timescale τ<sub>max</sub>. (b) Corresponding residuals from the models shown in (a). Light gray line shows the raw data. Dashed line indicates the zoomed window for the plot shown in (c).</figcaption></figure>
describes the observed effects of aftershocks between ~November 2015 and ~June 2016 when compared to the time series with the exponential models (Figure 6c vs. Figure 8c).

### 4.1.3. Implications of the Universal Relaxation Function and the Modeling Approaches

We used two methods to model the effect of NLME on the estimated $\delta v$. The first approach, using simple exponential functions, yielded a poor correction of the aftershocks-induced velocity changes (Figure 6c), despite using a dedicated relaxation timescale for aftershocks ($\tau_{aft} \sim 1.18-3.03$ days). In the second approach, we calibrated the universal relaxation function $R(t)$ (Snieder et al., 2017) with the same maximum relaxation time $\tau_{max}$ for all aftershocks events and the main shock of the Gorkha earthquake (Figure 8). The fit using $R(t)$ better captures the effect induced by aftershocks in the first part of the year 2016 (Figure 8). This agreement can be explained by the sensitivity of the $R(t)$ function to long relaxation times (Figure 7a), even when fitted on the early part of the relaxation curve following the aftershocks. Because of the apparent superiority of the $R$ models in this manuscript and considering the lower degrees of freedom to characterize the relaxation timescales $\tau$, we favor this approach.

We note that both our modeling approaches rely on the assumption of a linear summation of each induced perturbation. If the summation is realistic, it means that the ability to predict NLME requires the knowledge of strain history and not only the current state of the system. At our field site, this is important because our dataset starts 25 days after the $M_w$ 7.3 main aftershock of the 12 May 2015 (Figure 1). We did not correct for this event or any aftershocks occurring between the 25 April 2015 and the 6 June 2015. Nevertheless, we predict that most of the NLME effects are contained within the first approximate year ($R^{155}$ and $R^{250}$, Figure 8), a value consistent with the inferred recovery of landslide rates in the Bhote Koshi (~1 yr; Marc et al., 2021). If we assume that our inferred $\delta v$ estimated at rather high frequency (4–8 Hz) is a good proxy for shallow subsurface damage, this comparison with landsliding shows that our model is realistic and does not support a longer effect for NLME, such as inferred on model $R^{846}$ (Figure 8).

Another advantage in using $R(t)$ rather than the purely empirical approach is that the relaxation function may be more informative on the physical mechanisms responsible for NLME. The theory leading to the function $R(t)$ is based on an Arrhenius-like law (Snieder et al., 2017), in which the maximum relaxation timescale is given by

$$ \tau_{max} = A \exp \left( \frac{E_a^{max}}{k_B T} \right) \qquad (8) $$

in which $A$ is a prefactor, $E_a^{max}$ is an activation energy, $k_B$ is the Boltzmann constant, and $T$ is the temperature. $E_a^{max}$ can be interpreted as the barrier energy that needs to be overcome to reach a lower energy state from a metastable state. This barrier may correspond to characteristic contacts that undergo a particular thermally activated process in the slow dynamics phase, for example, dislocation creep or rearrangement transitions in granular composites. We obtained a good correction of the $\delta v$ data by using the same $\tau_{max}$ for events with variable initial perturbations, from PGV of $10^{-3}-10^{-2}$ m s$^{-1}$ for aftershocks, and in the range of ~5 × $10^{-1}$ m s$^{-1}$ for the Gorkha earthquake (Wei et al., 2018). Following Equation 7, this means that the nature of the physical mechanisms corresponding to $E_a^{max}$ and responsible for the longest relaxation timescale is independent from the intensity of ground shaking. Therefore, the relaxation timescales $\tau$ controlling slow dynamics in the probed medium would rather be a function of the ambient conditions such as temperature (Bekele et al., 2017), fluid content (Bittner & Popovics, 2021), or pre-existing damage (Astorga et al., 2018; Lyakhovsky et al., 1997, 2009) while the size of the initial excitation would control the number of characteristic broken contacts (Ostrovsky et al., 2019). This interpretation has important implications for the prediction of NLME and suggests that by studying the response induced by small events, one may predict the damage timescales induced by large dynamic strains. The investigation of a constant maximum relaxation timescale $\tau_{max}$ after dynamic strain perturbations of variable sizes could open a new perspective on NLME-induced changes: a complex physical phenomenon but with a potential deterministic behavior. This potential independence of $\tau$ from the ground shaking amplitude could explain the scattered relation between these variables when tested in field data (Viens et al., 2018). Considering the complexity of the relaxation processes in the Earth surface, the simple picture of a constant $\tau_{max}$ need to be tested in future works.

Velocity changes estimated at Chaku with lower frequency bands (1–2 and 2–4 Hz) exhibit smaller variations when plotted against the changes we report in this study at 4–8 Hz (Figure S7 in Supporting Information S1). This comparison indicates that the dominating NLME mechanisms are likely to be concentrated in near surface materials where smaller perturbations can induce strong changes at shallow depths (Qin et al., 2020) due to lower
confining pressure and more compliant materials. A relevant process is the re-arrangement of grains in soft spots of the near surface materials (Lieou et al., 2017). The higher susceptibility to dynamic strain of superficial loosely packed layers (Sawazaki et al., 2018) across a range of ground shaking intensities could explain the good fit of the $R(t)$ function of $\delta v$ after variable excitations using a constant $\tau_{max}$. At depths, the long term relaxation may happen in larger geological structures such as the fracture network, which is likely to expand through a large span of crustal depths in the tectonic regime of the Himalayas (Molnar et al., 2007). The simple picture of constant $\tau_{max}$ may be altered in these deeper layers where a variety of mechanisms can be activated, such as micro-crack closure (Brantut, 2015; Meyer et al., 2021), creeping of asperities (Aharonov & Scholz, 2018), or pressure-dissolution (Yasuhara & Elsworth, 2008). These mechanisms are generally activated above a certain dynamic strain threshold, required to break contacts under larger confining pressure. This can be justified by Amonton's law in which macroscopic friction is load dependent but also by recent observations that at the nanoscale, chemical bonds responsible for frictional force increase with normal load (Tian et al., 2020). Therefore, a constant universal $\tau_{max}$ might not hold if one compares different frequency bands that probe larger depths: The observation of diverse relaxations in the entire crust (Wang et al., 2019), the influence of confining pressure on velocity recovery (Meyer et al., 2021), and the example of larger NLME-induced changes at depths in fault zone (Qiu et al., 2020) support this direction. A spectrum of relaxation timescales responsible for slow dynamics (Shokouhi et al., 2017) may be needed to characterize different depths at any field site.

## 4.2. Hydrological Perturbation $\delta v_H^*$ After the Gorkha Earthquake

### 4.2.1. Static and Transient Model for Hydrological Changes

In the previous section, we modeled the $\delta v_{NLME}$ component by building three different relaxation models ($R^{155}$, $R^{846}$, and $R^{250}$) characterized by different maximum relaxation timescales $\tau_{max}$. In this section, we study the residuals obtained from these models (green lines, Figures 8a–8c) and compare them to the hydrological model of Illien et al. (2021); Figures 9a–9c.

We observe that the initial seismic velocity in June 2015 for the time-series corrected by models $R^{155}$ and $R^{250}$ are now comparable to the $\delta v$ level estimated in the following years at the same period (between +2.5% and +3% in the month of June), a feature that was not observed in the raw data (Section 3.1; Figures 2a–2c). This observation suggests that the residuals mainly describe the hydrological component $\delta v_H$, at least after a correction for NLME by models $R^{155}$ and $R^{250}$. To test this hypothesis, we use the initial seismic velocity of the residuals to calibrate the initial groundwater level condition used in the hydrological model of Illien et al. (2021) and plot the corresponding modeled velocities (red lines; Figures 9a–9c) without changing the original hydrological parameters inferred from the previous study. For the three NLME models, the agreement between the velocity residuals $\delta v$ and the hydrological model is greatly improved in comparison with the raw velocity data shown in Figure 2a. Nevertheless, the hydrological model predicts lower velocities than the observed residuals at the start of the 2015 monsoon, still causing a visible offset in the early part of the time series (Figures 9a–9c). In the model, lower velocities correspond to higher groundwater levels in the subsurface. This indicates that our hydrological prediction that was based on the velocities of the following years (2016–2018) overestimates the groundwater storage in the 2015 monsoon.

In the model of Illien et al. (2021), the groundwater drainage efficiency is proportional to the height of the hydraulic head $h(t)$ through a simple scaling:

$$ \frac{dh}{dt} = -a_{ss}h(t) + f(P(t), \text{vadose}(t)) \quad (9) $$

where $a_{ss}$ is the constrained steady state decay parameter that represents the average hydrological properties in the aquifer. $f$ is a function of the precipitation input $P(t)$ and the saturation condition in the vadose zone. A full derivation of the model is available in Illien et al. (2021). We test whether changing the parameter $a$ in a transient fashion following the Gorkha earthquake leads to better prediction of the velocity in 2015. We assume that the parameter $a$ is time-dependent and obeys the following evolution

$$ a(t) = a_{ss} \left( 1 + D \exp \left[ \frac{-(t - t_{\text{Gorkha}})}{\tau_{\text{hydro}}} \right] \right). \quad (10) $$
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Hydrological models vs. residuals from the NLME relaxation models. (a–c) Residuals from the models $R^{155}$, $R^{846}$, and $R^{250}$ are plotted in black. Red lines indicate the model from Illien et al. (2021) with an initial condition based on the residuals. (d–f) The green lines show the best fitting transient decay parameter $a(t)$. (g)–(i) The green lines indicate the modified hydrological models with the transient decay parameter $a(t)$ shown in plots (d–f). Close-ups on the data and the inferred models in 2015 in (j–l).</figcaption></figure>

We introduce a transient perturbation of the groundwater drainage with $Da_{ss}$ being the initial perturbation of the decay parameter at the date of Gorkha ($t_{Gorkha}$) and $\tau_{hydro}$ being the characteristic timescale for the recovery toward $a_{ss}$. The chosen form for $a(t)$ can be interpreted as a more efficient drainage of the groundwater table at early times after the earthquake, that progressively recovers toward a constant hydrological behavior. This is motivated by the observation that ground shaking can temporally increase stream discharge (Manga et al., 2003) and permeability measured in wells (Lai et al., 2014; Xue et al., 2013).

We minimize a least-square criterion to find the best fit between the velocities modeled with our time-dependent hydrological model and the $\delta v$ residuals obtained after removing the $\delta v_{NLME}$ synthetics. We explore a range of parameters for scaling $D$ and $\tau_{hydro}$ from Equation 9. For each NLME correction ($R^{155}$, $R^{846}$, and $R^{250}$), we report the best fitting transient decay parameter $a(t)$ in Figures 9d–9f and the associated modeled velocity changes in Figures 9g–9i. For all cases, introducing a transient increase of the groundwater drainage improves the fitting of $\delta v$ in the monsoon of 2015 (Figures 9j–9l). We find that best fitting values for the timescale $\tau_{hydro}$ range from 20 to 76 days and are therefore consistently one order of magnitude shorter than the maximum relaxation timescale $\tau_{max}$ applied in the NLME models. To compare the six inferred $\delta v_H$ models (Figures 9a–9c and 9g–9i), we compute their variances (Figure 10a). When no transient drainage parameter $a(t)$ is introduced, the model corrected with $R^{846}$ has the highest measured variance ($\sigma^2 = 4.3 \times 10^{-5}$) in comparison with $R^{155}$ and $R^{250}$ (both models around $\sigma^2 = 2.9 \times 10^{-5}$). With the transient parameter $a(t)$, $R^{250}$ is clearly the best fitting model ($\sigma^2 = 2.3 \times 10^{-5}$) while $R^{846}$ and $R^{155}$ both reproduce less than 90% of the $R^{250}$ variance based on their variance ratio (Figure 10b). Moreover, introducing the transient decay parameter $a(t)$ considerably improved the variance of model $R^{250}$ by a margin of ~20%. To test the significance of the fit, we perform a F-test (Text S2 in Supporting Information S1) between the model $R^{250}$ with no transient hydrological parameter (Figure 9c) and the model $R^{250}$ with the addition of the two parameters $D$ and $\tau_{hydro}$ (Figure 9i). We find that the introduction of $a(t)$ is statistically significant at 95% of confidence interval. Finally, we also explore a range of models with $\tau_{hydro}$ ranging from $10^1$ to $10^3$ days and a dedicated relaxation time for the Gorkha earthquake $\tau_{max}^G$ ranging from $2 \times 10^1$ to $5 \times 10^3$ days while retaining
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Variance of the hydrological models. (a) Absolute variance of the models, subscript h indicate models with the introduction of the transient parameter *a*(*t*). (b) Same plot as (a) but normalized with the best misfit value of the model *R*²⁵⁰.</figcaption></figure>

$\tau_{\text{max}}^{A} = 155$ days for the aftershocks and optimizing all the other parameters. Fitting all models characterized by the recovery timescales in this parameter space and minimizing a least-square criterion, we found that the best model is found for $\tau_{\text{hydro}} = 41$ days and $\tau_{\text{max}} = 450$ days (misfit space in Figure S8 in Supporting Information S1). These values are similar to our inferred model $R^{250}$ with $\tau_{\text{hydro}} = 35$ days.

The velocities in the 2015 monsoon are therefore better described with a transient drainage parameter, suggesting that the relaxation processes following the Gorkha earthquake affected hydrological properties. The duration of this perturbation for the best fitting model $R^{250}$ is ~6 months (Figure 9f). We address the validity of this claim and its implications in section 4.2.2.

### 4.2.2. Monitoring of Transient Hydrological Properties With Seismic Interferometry

A number of methods has been used to study how dynamic strain influences hydrological properties such as amplitudes and phase analysis in wells levels (Elkhoury et al., 2006; Xue et al., 2013), measurement of stream discharge (Manga et al., 2003), or monitoring with stable isotopes (Hosono et al., 2020). In parallel, a growing community of seismologists now use seismic interferometry to constrain groundwater storage (Kim & Lekic, 2019; Lecocq et al., 2017) but no attempt has been made to address earthquake hydrology topics with such methods. As mentioned before, this is partially due to the challenging decomposition of the several processes that influence seismic velocity. Our seismic interferometry analysis opens a window for monitoring transient hydrological behaviors on an intermediate spatial scale between point-based well measurements and catchment averaged isotopic and discharge analyses. We showed in Figure 9 that a transient increase in the drainage efficiency of the groundwater table improves our description of the seismic velocity changes in the 2015 monsoon in the aftermath of the Gorkha earthquake. This progress in the $\delta v$ fitting was tested for its significance, given the two parameters (*D* and $\tau_{\text{hydro}}$) we added to the original model of Illien et al. (2021; *F*-test in Text S2 in Supporting Information S1). However, the confidence interval of this test need to be taken carefully as our non-linear hydrological model may not produce normally distributed residuals, which are essential for parametric statistic tests (Gao, 2007). More interferometric datasets that are influenced by hydrological and seismic events should be tested in the future for cross-validation of our parametrization of transient properties *a*(*t*). Nevertheless, additional tests with a linear recovery for parameter *a*(*t*) (Figure S9 in Supporting Information S1) do not improve the variance observed with our exponential parametrization of (Equation 9; Figure S10 in Supporting Information S1).

Given the absence of additional constraint on $\tau_{\text{hydro}}$ in our study, the physical assumptions in our model are still supported by existing observations, such as a long lasting increase of permeability observed in other mountainous areas (Hosono et al., 2020), or the permeability healing phenomena observed for ~1 yr after the Wenchuan earthquake (Xue et al., 2013) and other South Californian earthquakes (Elkhoury et al., 2006). To further support our finding, we plot in Figure 11 the best fitting transient decay parameter (*a*(*t*) from model $R^{250}$) and an independent river stage height dataset from Bahrabise gauge station, located ~13 km downstream from our field site (blue star in Figure 1a). We compute the precipitation derived from the Global Precipitation Measurement data, IMERGHH 6B (Huffman et al., 2019) in a square of 100 m² upstream of the gauge (footprint in Figure S11 in
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Data from the Bahrabise gauge. The black solid line represents the stage height of the Bahrabise river (location of the gauge in Figure 1a). The height is corrected for an offset caused by the July 2016 glacial outburst flood (Cook et al., 2018). Precipitation estimated in the area are in blue and are obtained from the Global Precipitation Measurement data, IMERGHH 6B (Huffman et al., 2019). The green line shows the transient decay parameter $a(t)$ of the best fitting relaxation model $R^{250}$ ($\tau_{\text{hydro}}$ = 35 days). The background red color illustrates the period with enhanced permeability after Gorkha. G indicates the date of the Gorkha earthquake.</figcaption></figure>

Supporting Information S1) as it offers a suitable averaged measure to compare with the river height. The stage height measurement displays a co-seismic increase in discharge, supporting a release of mountain groundwater due to ground shaking (Wang et al., 2004). Additionally, the stage height has a clear co-evolution with monsoon precipitation with steep increase of the stage height that is concomitant with the onset of strong precipitation. However, the river gauge sensitivity to precipitation in 2015 seems relatively buffered, especially at the start of the 2015 monsoon when the onset of intense precipitation does not cause significant increase in stage height. This behavior looks to fade away rather quickly within the 2015 monsoon where the second pulse of precipitation induces a clear response in the stage height. A more permeable landscape with groundwater fluxes traveling more efficiently downstream or toward deeper layers at early times after the Gorkha earthquake is a plausible interpretation. Remarkably, the best fitting transient decay parameter $a(t)$ ($\tau_{\text{hydro}}$ = 35 days) recovers simultaneously to this observation, therefore showing a good agreement with this scenario (green line, Figure 11).

There is a limited number of experimental studies that links NLME and the evolution of hydrological properties. In limestones subjected to inelastic axial strain, it has been shown that after deformation, the seismic velocity was recovering for a few days but the permeability remained constant after a permanent increase due to damage (Brantut, 2015). The study mainly interprets the healing of velocities as the closure of micro-cracks porosity while the tortuosity of the pores network, which is the main control on permeability at the microscale (Kachanov & Sevostianov, 2005), remained unchanged. In this case, there is no co-evolution of hydrological properties with the slow dynamics phase. However, fluid flow in the field is thought to be largely controlled by the macroporosity (Baechle et al., 2004) and discrete fractures (Talwani et al., 2007). Notably, measurements of seismic velocity and permeability along a laboratory rock fracture both exhibit a phase of recovery after dynamic stressing (Shokouhi et al., 2020). At our field site, the estimated healing in hydrological properties from our model hints that the $\delta v_H$ variations could be contained in the fracture network. Possible mechanisms for permeability recovery includes fracture aperture modulation by destruction/creation of contact interfaces (Shokouhi et al., 2020) or colloids re-clogging (Mays & Hunt, 2007). In our relaxation models, $\tau_{\text{hydro}}$ is constantly shorter than the maximum recovery timescale $\tau_{\text{max}}$ used to correct NLME (Figure 9). This discrepancy between $\tau_{\text{max}}$ and $\tau_{\text{hydro}}$ may be explained by the non-linear relation between fracture aspect ratio and permeability (Ebigbo et al., 2016) due to percolation threshold. Another hypothesis would be that the changes responsible for these timescales are contained in different porosity units (e.g., micropores, macropores, fractures, etc.). New approach for characterizing NLME in the field are needed to disentangle these scenarios. At greater crustal depths, changes in hydrological properties may influence fluid migrations and low-frequency events, as observed several months after the 2011 Tohoku-Oki earthquake (Wang et al., 2021).

Future work may address the opposite role of water on relaxation processes. Pore water is generally considered to reduce frictional properties of interfaces, therefore, raising the susceptibility to ground shaking (Brenguier
et al., 2014). However, water also controls the rate of recovery through chemical reactions and changes in activation energies (Brantut, 2015; Liu & Szlufarska, 2010). The actual impact of such processes on ground velocity retrieved by seismic interferometry remains to be assessed.

## 5. Conclusions

In this study, we estimated relative seismic velocity changes $\delta v$ from single station cross-correlations in the aftermath of the 2015 $M_w$ 7.8 Gorkha earthquake for a duration of ~3 yr. Using the same characteristic relaxation timescales after the main shock and all the aftershocks (best fitting model for $\tau_{max}$ = 250 days), we corrected for the NLME effect. We found that the velocity changes evolve toward background values until the 2016 monsoon which suggests that most of the subsurface damage is recovered during the first year after the main shock. With the hydro-seismological model of Illien et al. (2021), we fitted the residual $\delta v$ corrected for NLME and inferred a shorter relaxation timescale $\tau_{hydro}$ that we attributed to an enhanced permeability of the subsurface that recovers gradually for ~6 months during the 2015 monsoon.

Special attention should be given when subtracting earthquake-induced velocity changes to constrain background hydrology as a transient behavior may be hidden in the data. Therefore, seismic interferometry studies may need to go beyond the assumed superposition of contributions $\delta v = \delta v_{NLME} + \delta v_H$ as the relaxation processes may affect the hydrological properties of the subsurface. In our study, we calibrated the nonlinear recovery with the relaxations triggered by the aftershocks, hence without biasing hydrological-induced velocity variations that can possibly be affected by ground shaking. Because of the importance of hydrological properties for freshwater resources, initiation of hillslope hazards (Iverson, 2000), and the frictional properties of fault zones (Talwani et al., 2007), we encourage the use of seismic techniques to estimate the hydrological response to large earthquakes using dense seismic arrays and multiple frequency bands.

## Data Availability Statement

The precipitation time-series for the Bhote Koshi observatory can be found at the following doi:10.5880/gfz.4.6.2021.002. Seismic data are available at the doi:10.14470/ka7560056170.

## Acknowledgments

This work was supported by the GFZ HART program. The authors thank Bhairab Sitaula and the HART Nepal team for the maintenance and installation of the seismic stations. Open access funding enabled and organized by Projekt DEAL.

## References

Adhikari, L. B., Gautam, U. P., Koirala, B. P., Bhattarai, M., Kandel, T., Gupta, R. M., et al. (2015). The aftershock sequence of the 2015 April 25 Gorkha-Nepal earthquake. *Geophysical Journal International*, 203(3), 2119–2124. https://doi.org/10.1093/gji/ggv412

Aharonov, E., & Scholz, C. H. (2018). A physics-based rock friction constitutive law: Steady state friction. *Journal of Geophysical Research: Solid Earth*, 123(2), 1591–1614. https://doi.org/10.1002/2016JB013829

Astorga, A., Guéguen, P., & Kashima, T. (2018). Nonlinear elasticity observed in buildings during a long sequence of earthquakes. *Bulletin of the Seismological Society of America*, 108(3), 1185–1198. https://doi.org/10.1785/0120170289

Baechle, G., Weger, R., Eberli, G., Miami, U., & Massaferro, J.-l. (2004). *The role of macroporosity and microporosity in constraining uncertainties and in relating velocity to permeability in carbonate rocks*. Paper presented at the SEG Int'l Exposition and 74th Annual Meeting, Denver. Retrieved from https://www.researchgate.net/publication/240611190_The_role_ofMacroporosity_and_microporosity_in_constraining_uncertainties_and_in_relatin_g velocity_to_permeability_in-carbonate rocks

Baillard, C., Lyon-Caen, H., Bollinger, L., Rietbrock, A., Letort, J., & Adhikari, L. B. (2017). Automatic analysis of the Gorkha earthquake aftershock sequence: Evidences of structurally segmented seismicity. *Geophysical Journal International*, 209(2), 1111–1125. https://doi.org/10.1093/gji/ggx081

Bekele, A., Birgisson, B., Ryden, N., & Gudmarsson, A. (2017). Slow dynamic diagnosis of asphalt concrete specimen to determine level of damage caused by static low temperature conditioning. *AIP Conference Proceedings*, 1806(1), 080012. https://doi.org/10.1063/1.4974637

Bittner, J. A., & Popovics, J. S. (2021). Mechanistic diffusion model for slow dynamic behavior in materials. *Journal of the Mechanics and Physics of Solids*, 150, 104355. https://doi.org/10.1016/j.jmps.2021.104355

Bontemps, N., Lacroix, P., Larose, E., Jara, J., & Taipe, E. (2020). Rain and small earthquakes maintain a slow-moving landslide in a persistent critical state. *Nature Communications*, 11(1), 1–10. https://doi.org/10.1038/s41467-020-14445-3

Bookhagen, B., & Burbank, D. W. (2010). Toward a complete Himalayan hydrological budget: Spatiotemporal distribution of snowmelt and rainfall and their impact on river discharge. *Journal of Geophysical Research*, 115(F3). https://doi.org/10.1029/2009JF001426

Brantut, N. (2015). Time-dependent recovery of microcrack damage and seismic wave speeds in deformed limestone. *Journal of Geophysical Research: Solid Earth*, 120(12), 8088–8109. https://doi.org/10.1002/2015JB012324

Brenguier, F., Campillo, M., Hadzioioannou, C., Shapiro, N. M., Nadeau, R. M., & Larose, E. (2008). Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations. *Science*, 321(5895), 1478–1481. https://doi.org/10.1126/science.1160943

Brenguier, F., Campillo, M., Takeda, T., Aoki, Y., Shapiro, N. M., Briand, X., et al. (2014). Mapping pressurized volcanic fluids from induced crustal seismic velocity drops. *Science*, 345(6192), 80–82. https://doi.org/10.1126/science.1254073

Brunello, C. F., Andermann, C., Marc, O., Schneider, K. A., Comiti, F., Achleitner, S., et al. (2020). Annually resolved monsoon onset and withdrawal dates across the Himalayas derived from local precipitation statistics. *Geophysical Research Letters*, 47(23), e2020GL088420. https://doi.org/10.1029/2020GL088420
Cook, K. L., Andermann, C., Gimbert, F., Adhikari, B. R., & Hovius, N. (2018). Glacial lake outburst floods as drivers of fluvial erosion in the Himalaya. *Science*, **362**(6410), 53–57. https://doi.org/10.1126/science.aat4981

Ebigbo, A., Lang, P. S., Paluszny, A., & Zimmerman, R. W. (2016). Inclusion-based effective medium models for the permeability of a 3D fractured rock mass. *Transport in Porous Media*, **113**(1), 137–158. https://doi.org/10.1007/s11242-016-0685-z

Elkhoury, J. E., Brodsky, E. E., & Agnew, D. C. (2006). Seismic waves increase permeability. *Nature*, **441**(7097), 1135–1138. https://doi.org/10.1038/nature04798

Elliott, J. R., Jolivet, R., Gonzalez, P. J., Avouac, J. P., Hollingsworth, J., Searle, M. P., et al. (2016). Himalayan megathrust geometry and relation to topography revealed by the Gorkha earthquake. *Nature Geoscience*, **9**(2), 174–180. https://doi.org/10.1038/ngeo2623

Gao, J. (2007). Nonlinear time series: Semiparametric and nonparametric methods. MRPA Paper No. 39563, (pp. 1–237). Monash University. https://doi.org/10.1111/j.1467-9892.2009.00642.x

Gassenmeier, M., Sens-Schönfelder, C., Eulenfeld, T., Bartsch, M., Victor, P., Tilmann, F., et al. (2016). Field observations of seismic velocity changes caused by shaking-induced damage and healing due to mesoscopic nonlinearity. *Geophysical Journal International*, **204**(3), 1490–1502. https://doi.org/10.1093/gji/ggv529

Gliozzi, A. S., Scalerandi, M., Anglani, G., Antonaci, P., & Salini, L. (2018). Correlation of elastic and mechanical properties of consolidated granular media during microstructure evolution induced by damage and repair. *Physical Review Materials*, **2**(1), 1–13. https://doi.org/10.1103/PhysRevMaterials.2.013601

Hobiger, M., Wegler, U., Shiomi, K., & Nakahara, H. (2014). Single-station cross-correlation analysis of ambient seismic noise: Application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake. *Geophysical Journal International*, **198**(1), 90–109. https://doi.org/10.1093/gji/ggu115

Hosono, T., Yamada, C., Manga, M., Wang, C. Y., & Tanimizu, M. (2020). Stable isotopes show that earthquakes enhance permeability and release water from mountains. *Nature Communications*, **11**(1), 2776. https://doi.org/10.1038/s41467-020-16604-y

Huffman, G. J., Bolvin, D. T., Nelkin, E. J., & Tan, J. (2019). *IMERG technical documentation*. Retrieved from https://gpm.nasa.gov/sites/default/files/document_files/IMERG_doc_190909.pdf

Illien, L., Andermann, C., Sens-Schönfelder, C., Cook, K. L., Baidya, K. P., Adhikari, L. B., et al. (2021). Subsurface moisture regulates Himalayan groundwater storage and discharge. *AGU Advances*, **2**(2), e2021AV000398. https://doi.org/10.1029/2021AV000398

Iverson, R. M. (2000). Landslide triggering by rain infiltration. *Water Resources Research*, **36**(7), 1897–1910. https://doi.org/10.1029/2000WR900090

Kachanov, M., & Sevostianov, I. (2005). On quantitative characterization of microstructures and effective properties. *International Journal of Solids and Structures*, **42**, 309–336. https://doi.org/10.1016/j.ijsolstr.2004.06.016

Kim, D., & Lekic, V. (2019). Groundwater variations from autocorrelation and receiver functions. *Geophysical Research Letters*, **46**, 13722–13729. https://doi.org/10.1029/2019GL084719

Lai, G., Ge, H., Xue, L., Brodsky, E. E., Huang, F., & Wang, W. (2014). Tidal response variation and recovery following the Wenchuan earthquake from water level data of multiple wells in the nearfield. *Tectonophysics*, **619**, 115–122. https://doi.org/10.1016/j.TECTO.2013.08.039

Lecocq, T., Longuevergne, L., Pedersen, H. A., Brenguier, F., & Stammler, K. (2017). Monitoring ground water storage at mesoscale using seismic noise: 30 years of continuous observation and thermo-elastic and hydrological modeling. *Scientific Reports*, **7**, 14241. https://doi.org/10.1038/s41598-017-14468-9

Lieou, C. K., Daub, E. G., Ecke, R. E., & Johnson, P. A. (2017). Slow dynamics and strength recovery in unconsolidated granular Earth materials: A mechanistic theory. *Journal of Geophysical Research: Solid Earth*, **122**(10), 7573–7583. https://doi.org/10.1002/2017JB014131

Liu, Y., & Szlufarska, I. (2010). Effect of trace moisture on friction. *Applied Physics Letters*, **96**(10), 101902. https://doi.org/10.1063/1.3356222

Lyakhovsky, V., Ben-Zion, Y., & Agnon, A. (1997). Distributed damage, faulting, and friction. *Journal of Geophysical Research B: Solid Earth*, **102**(12), 27635–27649. https://doi.org/10.1029/97jb01896

Lyakhovsky, V., Hamiel, Y., Ampuero, J. P., & Ben-Zion, Y. (2009). Non-linear damage rheology and wave resonance in rocks. *Geophysical Journal International*, **178**(2), 910–920. https://doi.org/10.1111/j.1365-246X.2009.04205.x

Manga, M., Beresnev, I., Brodsky, E. E., Elkhoury, J. E., Elsworth, D., Ingebritsen, S. E., et al. (2012). Changes in permeability caused by transient stresses: Field observations, experiments, and mechanisms. *Reviews of Geophysics*, **50**(2). https://doi.org/10.1029/2011RG000382

Manga, M., Brodsky, E. E., & Boone, M. (2003). Response of streamflow to multiple earthquakes. *Geophysical Research Letters*, **30**(5). https://doi.org/10.1029/2002gl016618

Marc, O., Behling, R., Andermann, C., Turowski, J. M., Illien, L., Roessner, S., et al. (2019). Long-term erosion of the Nepal Himalayas by bedrock landsliding: The role of monsoons, earthquakes and giant landslides. *Earth Surface Dynamics*, **7**(1), 107–128. https://doi.org/10.5194/esurf-7-107-2019

Marc, O., Hovius, N., Meunier, P., Uchida, T., & Hayashi, S. (2015). Transient changes of landslide rates after earthquakes. *Geology*, **43**(10), 883–886. https://doi.org/10.1130/G36961.1

Marc, O., Sens-Schönfelder, C., Illien, L., Meunier, P., Hobiger, M., Sawazaki, K., et al. (2021). Toward using seismic interferometry to quantify landscape mechanical variations after earthquakes. *Bulletin of the Seismological Society of America*, **113**(3), 1631–1649. https://doi.org/10.1785/0120200264

Mays, D. C., & Hunt, J. R. (2007). Hydrodynamic and chemical factors in clogging by montmorillonite in porous media. *Environmental Science and Technology*, **41**(16), 5666–5671. https://doi.org/10.1021/es062009s

Meyer, G. G., Brantut, N., Mitchell, T. M., Meredith, P. G., & Plümper, O. (2021). Time dependent mechanical crack closure as a potential rapid source of post-seismic wave speed recovery: Insights from experiments in Carrara Marble. *Journal of Geophysical Research: Solid Earth*, **126**(4), 1–29. https://doi.org/10.1029/2020JB021301

Molnar, P., Anderson, R. S., & Anderson, S. P. (2007). Tectonics, fracturing of rock, and erosion. *Journal of Geophysical Research*, **112**(3), 1–12. https://doi.org/10.1029/2005JF000433

Nakata, N., & Snieder, R. (2011). Near-surface weakening in Japan after the 2011 Tohoku-Oki earthquake. *Geophysical Research Letters*, **38**(17). https://doi.org/10.1029/2011GL048800

Ostrovsky, L., Lebedev, A., Riviere, J., Shokouhi, P., Wu, C., Stuber Geesey, M. A., et al. (2019). Long-time relaxation induced by dynamic forcing in geomaterials. *Journal of Geophysical Research: Solid Earth*, **124**(5), 5003–5013. https://doi.org/10.1029/2018JB017076

Pei, S., Niu, F., Ben-Zion, Y., Sun, Q., Liu, Y., Xue, X., et al. (2019). Seismic velocity reduction and accelerated recovery due to earthquakes on the Longmenshan fault. *Nature Geoscience*, **12**(5), 387–392. https://doi.org/10.1038/s41561-019-0347-1

Qin, L., Ben-Zion, Y., Bonilla, L. F., & Steidl, J. H. (2020). Imaging and monitoring temporal changes of shallow seismic velocities at the Garner valley near Anza, California, following the M7.2 2010 El Mayor-Cucapah earthquake. *Journal of Geophysical Research: Solid Earth*, **125**(1), 1–17. https://doi.org/10.1029/2019JB018070
Qiu, H., Hillers, G., & Ben-Zion, Y. (2020). Temporal changes of seismic velocities in the San Jacinto Fault zone associated with the 2016 Mw 5.2 Borrego Springs earthquake. *Geophysical Journal International*, 220(3), 1536–1554. https://doi.org/10.1093/gji/ggz538

Roback, K., Clark, M. K., West, A. J., Zekkos, D., Li, G., Gallen, S. F., et al. (2018). The size, distribution, and mobility of landslides caused by the 2015 Mw7.8 Gorkha earthquake, Nepal. *Geomorphology*, 301, 121–138. https://doi.org/10.1016/j geomorph.2017.01.030

Sawazaki, K., Saito, T., & Shiomi, K. (2018). Shallow temporal changes in S wave velocity and polarization anisotropy associated with the 2016 Kumamoto earthquake sequence, Japan. *Journal of Geophysical Research: Solid Earth*, 123(11), 9899–9913. https://doi.org/10.1029/2018JB016261

Sawazaki, K., Sato, H., Nakahara, H., & Nishimura, T. (2009). Time-lapse changes of seismic velocity in the shallow ground caused by strong ground motion shock of the 2000 Western-Tottori earthquake, Japan, as revealed from coda deconvolution analysis. *Bulletin of the Seismological Society of America*, 99(1), 352–366. https://doi.org/10.1785/0120080058

Sens-Schönfelder, C., & Brenguier, F. (2019). Noise-based monitoring. In N. Nakata, L. Gualtieri, & A. Fichtner (Eds.), *Seismic ambient noise* (pp. 267–301). Cambridge University Press.

Sens-Schönfelder, C., Pomponi, E., & Peltier, A. (2014). Dynamics of Piton de la Fournaise volcano observed by passive image interferometry with multiple references. *Journal of Volcanology and Geothermal Research*, 276, 32–45. https://doi.org/10.1016/j JVOLGORES.2014.02.012

Sens-Schönfelder, C., & Wegler, U. (2006). Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia. *Geophysical Research Letters*, 33(21). https://doi.org/10.1029/2006GL027797

Shi, Z., Wang, G., Manga, M., & Wang, C. Y. (2015). Continental-scale water-level response to a large earthquake. *Geofluids*, 15(1–2), 310–320. https://doi.org/10.1111/gfl.12099

Shokouhi, P., Jin, J., Wood, C., Rivière, J., Madara, B., Elsworth, D., et al. (2020). Dynamic stressing of naturally fractured rocks: On the relation between transient changes in permeability and elastic wave velocity. *Geophysical Research Letters*, 47(1), 1–10. https://doi.org/10.1029/2019GL083557

Shokouhi, P., Rivière, J., Guyer, R. A., & Johnson, P. A. (2017). Slow dynamics of consolidated granular systems: Multi-scale relaxation. *Applied Physics Letters*, 111(25). https://doi.org/10.1063/1.5010043

Snieder, R., Sens-Schönfelder, C., & Wu, R. (2017) The time dependence of rock healing as a universal relaxation process, a tutorial. *Geophysical Journal International*, 208(1), 1–9. https://doi.org/10.1093/gji/ggw377

Taira, T., Brenguier, F., & Kong, Q. (2015). Ambient noise-based monitoring of seismic velocity changes associated with the 2014 Mw 6.0 South Napa earthquake. *Geophysical Research Letters*, 42(17), 6997–7004. https://doi.org/10.1002/2015GL065308

Talwani, P., Chen, L., & Gahalaut, K. (2007). Seismogenic permeability, $k_s$. *Journal of Geophysical Research: Solid Earth*, 112(7). https://doi.org/10.1029/2006JB004665

TenCate, J. A., Smith, E., & Guyer, R. A. (2000). Universal slow dynamics in granular solids. *Physical Review Letters*, 85(5), 1020–1023. https://doi.org/10.1103/PhysRevLett.85.1020

Tian, K., Li, Z., Liu, Y., Gosvami, N. N., Goldsby, D. L., Szlufarska, I., et al. (2020). Linear aging behavior at short timescales in nanoscale contacts. *Physical Review Letters*, 124(2), 26801. https://doi.org/10.1103/PhysRevLett.124.026801

Viens, L., Denolle, M. A., Hirata, N., & Nakagawa, S. (2018). Complex near-surface rheology inferred from the response of greater Tokyo to strong ground motions. *Journal of Geophysical Research: Solid Earth*, 123(7), 5710–5729. https://doi.org/10.1029/2018JB015697

Wang, C. Y., & Manga, M. (2015). New streams and springs after the 2014 Mw6.0 South Napa earthquake. *Nature Communications*, 6, 6–11. https://doi.org/10.1038/ncomms8597

Wang, C. Y., Wang, C. H., & Manga, M. (2004). Coseismic release of water from mountains: Evidence from the 1999 (Mw = 7.5) Chi-Chi, Taiwan, earthquake. *Geology*, 32(9), 769–772. https://doi.org/10.1130/G20753.1

Wang, Q. Y., Campillo, M., Brenguier, F., Lecointre, A., Takeda, T., & Hashima, A. (2019). Evidence of changes of seismic properties in the entire crust beneath Japan after the Mw 9.0, 2011 Tohoku-oki earthquake. *Journal of Geophysical Research: Solid Earth*, 124(8), 8924–8941. https://doi.org/10.1029/2019JB017803

Wang, Q. Y., Campillo, M., Brenguier, F., Lecointre, A., Takeda, T., & Yoshida, K. (2021). Seismic evidence of fluid migration in northeastern Japan after the 2011 Tohoku-oki earthquake. *Earth and Planetary Sciences Letters*, 563, 116894. https://doi.org/10.1016/j.epsl.2021.116894

Wegler, U., & Sens-Schönfelder, C. (2007). Fault zone monitoring with passive image interferometry. *Geophysical Journal International*, 168(3), 1029–1033. https://doi.org/10.1111/j.1365-246X.2006.03284.x

Wei, S., Chen, M., Wang, X., Graves, R., Lindsey, E., Wang, T., et al. (2018). The 2015 Gorkha (Nepal) earthquake sequence: I. Source modeling and deterministic 3D ground shaking. *Tectonophysics*, 722, 447–461. https://doi.org/10.1016/j.tecto.2017.11.024

Wu, C., Peng, Z., & Ben-Zion, Y. (2010). Refined thresholds for non-linear ground motion and temporal changes of site response associated with medium-size earthquakes. *Geophysical Journal International*, 182(3), 1567–1576. https://doi.org/10.1111/j.1365-246X.2010.04704.x

Xue, L., Li, H.-B., Brodsky, E. E., Xu, Z.-Q., Kano, Y., Wang, H., et al. (2013). Continuous permeability measurements record healing inside the Wenchuan earthquake fault zone. *Science*, 340(6140), 1555–1559. https://doi.org/10.1126/science.1237237

Yasuhara, H., & Elsworth, D. (2008). Compaction of a rock fracture moderated by competing roles of stress corrosion and pressure solution. *Pure and Applied Geophysics*, 165(7), 1289–1306. https://doi.org/10.1007/s00024-008-0356-2