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

HAL
open science

# Subsurface evolution of the seasonally frozen ground on the northeastern Tibetan Plateau from a perspective of seismic interferometry

Chang Yu, Qing-Yu Wang, Jianxin Ma, Huajian Yao

► To cite this version:

Chang Yu, Qing-Yu Wang, Jianxin Ma, Huajian Yao. Subsurface evolution of the seasonally frozen ground on the northeastern Tibetan Plateau from a perspective of seismic interferometry. Geophysical Journal International, 2025, 241, pp.308-325. 10.1093/gji/ggaf042. insu-05082248

HAL Id: insu-05082248

https://insu.hal.science/insu-05082248v1

Submitted on 23 May 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Subsurface evolution of the seasonally frozen ground on the northeastern Tibetan Plateau from a perspective of seismic interferometry

Chang Yu<sup>1</sup>, Qing-Yu Wang<sup>2,3</sup>, Jianxin Ma<sup>4</sup> and Huajian Yao<sup>1,5</sup>

Accepted 2025 January 24. Received 2025 January 21; in original form 2024 June 11

## SUMMARY

The Tibetan Plateau, a critical region influencing both local and global atmospheric circulation, climate dynamics, hydrology and terrestrial ecosystems, is undergoing climate-driven changes, including glacial retreat, permafrost thaw and groundwater changes. Despite its importance, implementing continuous and systematic observations has been challenging due to the area's high altitude and extreme climate conditions. In this context, seismic interferometry emerges as a cost-effective method for the continuous monitoring of subsurface structural changes driven by environmental factors and internal geophysical processes. We investigate subsurface evolution using four years of seismic data from nine stations on the northeastern Tibetan Plateau, by applying coda wave interferometry across multiple frequency bands. Our findings highlight seismic velocity changes within the frequency bands 5–10, 0.77–1.54, and 0.25–0.51 Hz, revealing depth-dependent seasonal and long-term changes. Near-surface and deeper strata exhibit similar seasonal patterns, with velocities increasing in winter and decreasing in summer driven by changes in hydrological processes, while intermediate ice-water phase strata show contrasting behaviour due to thermal elastic strain. Long-term trends suggest that the upper subsurface layer is affected by melting water and precipitation originating from Kunlun Mountains, whereas deeper layer reflect groundwater level variations influenced by climate change and human activities. This study provides insights into the environmental evolution of the Tibetan Plateau and its impact on managing local groundwater resources.

**Key words:** Interferometry; Time-series analysis; Hydrogeophysics; Coda waves; Seismic noise.

## 1 INTRODUCTION

Climate change presents an acute global challenge (Jorgenson et al. 2001; Taylor et al. 2013), with its impacts becoming increasingly pronounced, particularly in regions such as the Tibetan Plateau (Yang et al. 2019; Wang et al. 2021). Commonly referred to as the ‘Asia Water Tower’, the Tibetan Plateau is crucial in providing water resources for approximately two billion individuals (Li et al. 2022; Yao et al. 2022). The thawing of permafrost in this region, driven by climate change, could directly influence the regional hydrological processes and subsequently threaten the water resource security (Gao et al. 2022). Additionally, the thawing of frozen soil accelerates the release of carbon dioxide previously stored in the

frozen ground (Wang et al. 2020), exerting a pronounced effect on both regional and global climate patterns (Hodgkins et al. 2014). Gaining insights into the temporal variations of subsurface medium within these permafrost regions is essential to unravel the complex interactions among climate, ecosystems and hydrological processes (Viviroli et al. 2007; Bibi et al. 2019; Gao et al. 2022; Li et al. 2022). While traditional geological approaches involving drilling or observatories can furnish valuable data on subsurface medium changes, these approaches often grapple with significant economic barriers and spatial resolution constraints (Green et al. 2011; Biskaborn et al. 2015, 2019). Established geodetic techniques, including the Gravity Recovery and Climate Experiment (GRACE), the interferometric synthetic aperture radar and the global positioning

<footnote n="1"><sup>1</sup>Laboratory of Seismology and Physics of Earth's Interior, School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China. E-mail: hjyao@ustc.edu.cn</footnote>

<footnote n="2"><sup>2</sup>Institut Terre et Environnement de Strasbourg (ITES) UMR 7063, Université de Strasbourg, EOST/CNRS, Strasbourg 67084, France.
Email:qingyu.wang@unistra.fr</footnote>

<footnote n="3"><sup>3</sup>ISTerre, Université Grenoble Alpes, University Savoie Mont Blanc, CNRS, IRD, University Gustave Eiffel, Grenoble 38610, France</footnote>

<footnote n="4"><sup>4</sup>Earthquake Administration of Qinghai Province, Xining 810001, China</footnote>

<footnote n="5"><sup>5</sup>Mengcheng National Geophysical Observatory, University of Science and Technology of China, Mengcheng Anhui 233500, China</footnote>
system (GPS), are effective for monitoring surface deformation, but they fall short in providing direct depth-specific information (Chen & Liu 2016; Pan et al. 2018; Riel et al. 2018). In the context of complex mountainous landscapes and harsh climate conditions, traditional geophysical methods encounter distinct challenges in continuously monitoring subsurface variations at a relatively low cost.

Ambient noise-based seismic interferometry has emerged as an important advancement, providing direct insights into the dynamics of subsurface physical properties. The noise-based seismic monitoring approach uses coda waves from cross-correlation functions, derived from continuous seismic recordings. This method facilitates the tracking of temporal variations in seismic wave velocities, thereby shedding light on alterations in the elastic properties of subsurface medium. Empirical studies have shown that changes detected through this noise-based seismic method adeptly capture changes related to volcanic activities (Brenguier et al. 2008b, 2014; Wu et al. 2020), earthquakes (Brenguier et al. 2008a; Acarel et al. 2014; Wang et al. 2019) and environmental fluctuations (Hillers et al. 2014; Wang et al. 2017). Within permafrost regions, this approach uncovers distinct seasonal freeze-thaw cycles and permafrost degradation that arise from temperature fluctuations (James et al. 2019; Lindner et al. 2021). This approach enables continuous, cost-effective probing of subsurface evolution, providing complementary observations at various depths compared to traditional geological methods (Poli et al. 2020; Delouche & Stehly 2023).

In this study, we explore a unique data set from the Tibetan Plateau, where the seasonally frozen ground serves as a natural laboratory for probing subsurface evolution. We examine the effects of climatic changes on local subsurface evolution using seismic interferometry with nine stations deployed across the northeastern Tibetan Plateau. The analysis of seismic velocity changes reveals depth-dependent patterns, prompting an investigation into the primary factors governing these distinct seismic patterns. By integrating meteorological and hydrological data, our findings highlight the critical role that complex hydrothermal processes, influenced by climate change, play in driving subsurface deformation in the northeastern Tibetan Plateau.

# 2 DATA AND METHODS

## 2.1 Experimental setup

We analyse continuous ambient noise records from nine short-period stations located in the Nuomuhong region of the southeastern Qaidam Basin. The arid climate renders local economy and ecological security fragile and sensitive to global climate change (Song et al. 2014; Wang et al. 2014; Zhang et al. 2023). The seismic stations, forming a distinctive “V” pattern, lie at elevations ranging from 3200 to 3300 m and occupy strategic positions upstream of the alluvial fan. This setup provides a unique opportunity to understand subsurface evolutions in the context of global climate change.

Each of these stations is outfitted with uniform short-period instruments operating at a sampling frequency of 100 Hz. The instruments are mounted on concrete columns to enhance coupling with the ground and minimize interference from the complex local environmental conditions.

It is important to note that the stations occasionally experience unexpected data loss due to events such as signal transmission errors

(Fig. S1, Supporting Information). However, signal transmission errors do not always result in complete data loss. They can also lead to intermittent lapses in GPS time reception. Since GPS signals serve as highly accurate external reference clocks ensuring precise timing (Sens-Schönfelder 2008; Hirose & Ueda 2023), such interruptions can cause slight clock discrepancies between different stations. These clock errors may introduce positive or negative delays in the cross-correlations (Stehly et al. 2007; Xu & Song 2009). To mitigate the impact of these timing inaccuracies, we use autocorrelation functions derived from the vertical components of the stations in our analysis.

Seismic stations are nestled among local rivers that originate from the Kunlun Mountains and the Nuomuhong river is the most prominent river in this region (Fig. 1). The local river flow patterns exhibit a clear correlation with temperature fluctuations. The local rainy season typically spans from June to August, with precipitation increasing and evaporation decreasing as the elevation rises (Wang et al. 2008; Xiao et al. 2017). Moreover, given the pronounced evaporation rates (Fig. S2, Supporting Information) in the region, a significant proportion of the hydrological sustenance of these rivers is attributed to groundwater and surface runoff, which is facilitated by precipitation and melting water in the Kunlun Mountains (Wang et al. 2008; Zhang 2013; Xiao et al. 2017).

## 2.2 Data preparation

We segment the continuous data into daily intervals after removing both the mean value and the long-term trend from the data. We apply spectrum whitening (Shapiro et al. 2006) with running-absolute-average normalization across multiple frequency bands to reduce interference from noise sources (see Text S1 in the Supporting Information) and perform temporal normalization to mitigate the impact of abnormal signals and seismic events within frequency bands: 0.1–0.2, 0.2–0.5, 0.5–1, 1–10, and 10–30 Hz (Yao et al. 2006, 2011; Bensen et al. 2007; Zhang et al. 2018). We then combine signals to cover broad frequency bands ranging from 0.1–30 Hz.

We compute autocorrelation functions in one-hour windows with a 60 per cent overlap, then stack them daily to produce consistent autocorrelation results. To improve the signal-to-noise ratio of the autocorrelation functions, we apply the singular value decomposition method (Melo et al. 2013; Moreau et al. 2017) and retain eight singular values (Fig. S3, Text S2 and S3, Supporting Information). This parameter ensures optimal stability while preserving the desired temporal resolution (Figs S4–S8, Supporting Information). Additionally, to mitigate the influence of short-term signals and further improve data quality, we employ a time-weighted sliding average approach. Considering the similarities between neighbouring days, we opt for a time-weighted sliding average to smooth our autocorrelation functions. The formula for the weighted sliding average is as follows:

$$S_i(t) = \frac{1}{\sum_{k=i-N}^{k=i+N} w_k} \sum_{k=i-N}^{k=i+N} w_k s_k(t). \quad (1)$$

In eq. (1), $s_k(t)$ represents the daily correlation functions, $N$ represents the time length of the weighted running average and $w_k(t)$ is the weight assigned to each day. The weight is calculated as:

$$w_k = e^{-\frac{|k-i|}{L}}. \quad (2)$$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Map of the study area and station locations. (a) The locations of nine seismic stations (marked by triangles) are outlined by the square on the northeastern edge of the Tibetan Plateau. The dot represents the position of the meteorological station. Light lines depict local rivers, while the darker line outlines the primary local river, the Nomuhong River. (b) Zoom-in on the locations of seismic stations NMH1 to NMH9.</figcaption></figure>

Here, *L* represents the correlation length, which characterizes the similarity between different days. In this paper, we set *L* = 31. The choice of *N* has little effect on the final result given that the value of *N* is two times larger than *L* (Fig. S9, Supporting Information). The weighted sliding average method proves to be more effective in reducing the effects of transient fluctuations on the autocorrelation functions, thereby providing more stable results (Fig. S10, Supporting Information).

## 2.3 Monitoring of seismic velocity changes

We monitor seismic velocity changes using the wavelet cross-spectrum method (Mao *et al.* 2020). This method holds an advantage in analysing multiple frequency bands due to its optimal time-frequency resolution. We measure seismic phase differences using a 20 to 60 s window on autocorrelation functions, representing later arrivals and considered as multiple scattered waves, referred to as coda waves. These coda waves are less affected by non-uniformly distributed noise sources (Froment *et al.* 2010; Colombi *et al.* 2014). Additionally, they are sensitive to subtle changes in the medium, as the phase shift accumulates with increasing lapse time over correlation functions (Larose *et al.* 2006; Hadzioannou *et al.* 2009). As a supplementary strategy, we also assess seismic velocity changes using two other measurement methods: the stretching method (Lobkis & Weaver 2003; Sens-Schönfelder & Wegler 2006; Xie *et al.* 2023) and the moving window cross-spectrum (MWCS) method (Clarke *et al.* 2011). We set the sliding window length in MWCS to the 1.5 times the largest period of each frequency band, maintaining the 20 to 60 s time window range in both methods (Text S4, Supporting Information). Monitoring results are comparable across all three approaches, as shown in Figs S11–S12 (Supporting Information).

# 3 RESULTS

## 3.1 Time-series of seismic velocity changes

Across the three frequency bands, we observe pronounced and consistent seasonal changes in seismic wave velocities (Figs 2a–c) from all nine stations. We have excluded results from seismic station NMH1 due to its frequent data loss, particularly before 2016, resulting from signal transmission errors (Figs S1 and S13, Supporting Information). Additionally, we remove the results for days when more than two stations have missing data records.

The seismic monitoring results from 5–10 and 0.25–0.51 Hz bands display similar seasonal patterns, with velocities peaking in January and reaching their lowest values in July. In contrast, the 0.77–1.54 Hz band exhibits seasonal variations that oppose those observed in the other two bands.

Seasonal fluctuations in seismic velocity are typically influenced by local environmental factors such as temperature (Lindner *et al.* 2021), atmospheric pressure (Silver *et al.* 2007) and groundwater dynamics (Clements & Denolle 2018; Mao *et al.* 2022). Therefore, we collect meteorological data from a weather station near the seismic stations (Fig. 1) to help identify how environmental variables contribute to temporal variations within these layers. The environmental variables, depicted in Figs 2(d)–(f), include time-series of changes in air temperature, precipitation and air pressure, respectively. All of them display distinct annual patterns that likely contribute to the observed seasonal changes in seismic velocity.

Beyond the seasonal fluctuations, the results for the 5–10 Hz band indicate a significant overall decrease in seismic wave velocity, while the 0.25–0.51 Hz band reveals a subtle increase. Additionally, the intermediate frequency band of 0.77–1.54 Hz demonstrates a very subtle, long-term increasing trend.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time-series of seismic observations and meteorological recordings. Panels (a), (b) and (c) represent the seismic velocity changes with the wavelet cross-spectrum method in the frequency ranges of 5–10, 0.77–1.54 and 0.25–0.51 Hz, respectively. These changes are averaged over results with a standard deviation in a light shade from all seismic stations, excluding station NMH1. The dashed lines represent the linear trend of the velocity change results for each frequency band. (d) Air temperature smoothed with a 31-d sliding window; (e) daily precipitation; and (f) atmospheric pressure smoothed with a 31-d sliding window.</figcaption></figure>

In addition to environmental factors, variations in noise sources may disturb the observed seismic velocity changes (Zhan *et al.* 2013). To determine whether our results are affected by noise source variations, we investigate the daily power spectral density (PSD) and root mean square (rms) of seismic records at station NMH6 (Figs S14 and S15, Supporting Information). There are no obvious seasonal patterns or long-term variations in the PSD and rms within our analysed frequency bands (0.25–0.51, 0.77–1.54, 5–10 Hz). Additionally, coda waves are less sensitive to noise source variations (Hadziioannou *et al.* 2009; Colombi *et al.* 2014). We therefore are able to capture how the seismic velocity changes in the subsurface.

The distinct seasonal and long-term variations across frequency bands likely reveal depth-dependent characteristics associated with changes in the physical properties of the medium (Fig. 3b). The depth information of coda waves depends on a combination of both body and surface wave depth sensitivity (Obermann *et al.* 2013, 2016, 2019). However, determining this information remains

challenging due to the unknown parameters related to the frequency-dependent mean-free path and its evolution over lapse time (Wang *et al.* 2019; Barajas *et al.* 2022). The coda wave window, spanning 20–60 s, allows us to simplify the assumption that it is predominantly influenced by the sensitivity kernel of surface waves.

We did not detect obvious higher-mode Rayleigh wave energy from the cross-correlations, suggesting that higher-mode energies are not prominent in our study areas. Therefore, to estimate the dominant depth, we compute only the fundamental-mode sensitivity kernels of Rayleigh waves (Herrmann 2013) using the Crust1.0 model (Laske *et al.* 2013).

We normalize the sensitivity kernels by its absolute maximum value within its frequency band to emphasize the most prominent depth.

Fig. 3(b) shows that seismic velocity changes in the 5–10 Hz frequency band mainly reflect changes within the first 50 m, and the 0.77–1.54 Hz band reflects changes above approximately 500 m.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Estimated sensitivity kernel and subsurface structure schematic. It is important to note that we employ a combination of linear and logarithmic scales for better visualization: a linear scale of 0–10 m and a logarithmic scale beyond 10 m in panels (a) and (b). (a) Velocity model used in the study; (b) normalized sensitivity kernel for each frequency band; and (c) schematic cross-section of the rock types beneath the seismic stations.</figcaption></figure>

Seismic results from the lowest frequency band of 0.25–0.51 Hz are controlled by changes deeper than 1 km. Hence, we interpret the measurements obtained from the three frequency bands as indicative of seismic velocity variations within the shallow, intermediate and deep layers.

In addition to calculating depth sensitivity, we also benefit from previous local geological investigations (Wang *et al.* 2008; Zhang 2013). These studies have unveiled the existence of a seasonal active layer at depths of approximately 1–2 m (the first layer in Fig. 3c). This layer undergoes freeze–thaw cycles, significantly impacting water content in the shallow regions (Wang *et al.* 2007, 2008; Zhang 2013). Beneath this active layer, drilling data (Wang *et al.* 2008; Zhang 2013) has uncovered unconsolidated deposits consisting of Holocene and Upper Pleistocene alluvial layers, with a thickness of approximately 200 m (the second layer in Fig. 3c). These unconsolidated layers primarily comprise gravel pebbles, gravel, medium-fine sand and clay soil. The loose structure of this shallow layer allows for significant fluctuations in water content, as documented in previous research (Tan *et al.* 2021).

There is a clay layer beneath the unconsolidated layer and above the water table at ~ 500 m (Wang *et al.* 2008; Zhang 2013). The clay layer compose of Middle Pleistocene glaciofluvial deposits characterized by the intermittent presence of grey gravel, sand, clay

layers and occasional pebble clasts (the third layer in Fig. 3c). The thickness of the exposed deposits varies between approximately 50 and 100 m (Zhang 2013). These sediments, rich in clay content, have a compact structure with mild cementation, resulting in low permeability. In certain conditions, the clay layer can become desiccated and fissured above the water table, further impacting the layer's properties (Grisak & Cherry 1975; Ruland *et al.* 1991).

# 4 DISCUSSION

## 4.1 Seismic velocity changes in the shallow layer

Seismic velocity changes in the shallow layer (5–10 Hz) are predominantly driven by seasonal effects, with velocities increasing in winter and decreasing in summer. These changes are sensitive to variations in physical properties of medium within the first 50 m at depth (Fig. 3).

We first compare the time-series of seismic velocity with local meteorological data, as shown in Figs 2(a), (d), (e) and (f).

In the shallow medium, we observe a negative correlation between temperature and the relative change in seismic velocity. This suggests that thermoelastic strain, induced by temperature
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Seasonal variations in shallow seismic velocity compared with meteorological data. (a) Map of the study area's topography and local rivers, with the locations of nine seismic stations represented by the triangle outlined by the square. The dot represents the position of the meteorological station. Local rivers are shown with light lines, and the Nuomuhong river is outlined by a darker line, indicating the major water system in the study area. (b) Comparison of annual patterns of shallow seismic velocity changes with local temperature, precipitation and river flow data averaged over multiple years. Temperature and precipitation data are recorded by the weather station identified in panel (a). River flow data is sourced from Wang *et al.* (2008). The y-axis direction for shallow velocity change is inverted for clarity. (c) Distribution of correlation coefficients between annual shallow velocity changes and annual net precipitation (precipitation minus evaporation) derived from ERA5 over the study period (2015–2018). (d) Distribution of correlation coefficients between annual shallow velocity changes and annual snow melting water derived from ERA5 over the study period (2015–2018).</figcaption></figure>

variations, is unlikely to be the primary controlling factor for seismic velocity in this layer (Ben-Zion & Leary 1986; Hillers et al. 2015). This is because, if velocity changes were dominated by thermomechanic strain, we would expect a positive correlation with temperature accompanied by an approximate phase difference of ~T/8, where T represents the period (one year for seasonal effects). The loose structures in the shallow layer and the locations of the stations on sloping surfaces (Fig. 1) may contribute to structural instability, which could reduce the effects of temperature (Section 4.3).

We also examine the relationship between atmospheric pressure variations and seismic velocity changes. Considering the sensitivity of seismic velocity to atmospheric pressure changes, estimated to be $3.2 \times 10^{-8}$ (Olivier & Brenguier 2016), an annual atmospheric pressure amplitude of approximately 600 Pa (Fig. 2f) would result in velocity changes smaller than $2 \times 10^{-5}$, which is inconsistent with our observations. This analysis suggests that atmospheric pressure, too, is unlikely to be a dominant factor affecting shallow layer seismic velocity changes.

After excluding temperature fluctuations and atmospheric pressure variations as the primary drivers, we observe a notable negative

correlation between seismic velocity changes and local precipitation records (Figs 2a and e).

This suggests that changes in water content may negatively impact seismic velocity due to alterations in pore pressure (Wang et al. 2024) within the shallow medium (Fig. 4b).

The change in water content in the shallow medium can primarily be attributed to two factors. The first factor is direct precipitation, which directly impacts the moisture content of the shallow medium. However, localized high rates of evaporation limit the influence of local precipitation (Fig. S2, Supporting Information). The second factor involves runoff from sources such as the Nuomuhong river and other local rivers, facilitated by the high permeability of unconsolidated deposits in the shallow medium. The flow of the Nuomuhong river, the main surface water system around the stations, exhibits seasonal patterns similar to those of temperature and has an inverse relationship with annual shallow velocity changes (Fig. 4b). Given that runoff is predominantly recharged by precipitation and melting water in the Kunlun Mountains (Wang et al. 2008; Zhang 2013; Xiao 2018), similar to the Nuomuhong river (Fig. S16, Supporting Information), it can be inferred that runoff
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparative analysis of soil moisture and shallow seismic velocity variations, combined with long-term trends in local hydrological data near seismic stations. (a) Soil moisture variations at the average location of the array, with the shaded area representing the study period (2015–2018). The short line indicates the linear trend of soil moisture throughout this interval. Note that the direction of y-axis is inverted for clarity. (b) Comparison of shallow velocity change and soil moisture data over study period. Note that the direction of y-axis of soil moisture data is inverted for clarity. (c) Distribution of linear long-term trends in soil moisture during the study period. (d) Distribution of linear long-term trends in net precipitation (precipitation minus evaporation) throughout the study period. (e) Distribution of linear trends of melting water during the study period.</figcaption></figure>

near the seismic stations experiences similar seasonal variations, consequently affecting shallow seismic velocity changes.

To verify our assumptions, we gather monthly averaged hydrological data from the fifth generation of global climate and weather data set provided by the European Centre for Medium-Range Weather Forecasts (ERA5-Land), which offers a resolution of 0.1° × 0.1° (Muñoz Sabater 2019). We choose this data set for its comprehensive coverage and relatively high spatial resolution, essential for analysing seasonal variation in runoff patterns in high-altitude areas, where meteorological observations are typically insufficient due to extreme natural conditions.

We compute the correlation coefficients between the annual monthly average changes in shallow seismic velocity and the annual monthly average net precipitation (precipitation minus evaporation) as well as separately with snow melting water (Figs 4a and b). We map the correlation coefficients in Figs 4(c) and (d). We observe strong negative correlations with seasonal variations across most of the Kunlun Mountains around the study region. Conversely, the low-altitude areas located in the northern part of the seismic stations exhibit a different pattern, due to the reduced precipitation and melting water alongside increased evaporation (Fig. S16, Supporting Information), which falls outside the primary focus of

our study. Considering the similar seasonal patterns observed between the main local water system (the Nuomuhong river) and these two primary recharge sources (net precipitation and melting water in the Kunlun Mountains), we suggest that the runoff from these sources follows similar variations patterns and predominantly influences the water content in the shallow medium around the seismic stations. Given the high permeability of shallow medium, as illustrated in Fig. 3(c), we assume that the water content variation pattern within the top 50 m are approximately equivalent to those in the shallow soil moisture. Consequently, we use soil moisture data from ERA5-Land model to further validate our hypothesis.

The soil moisture data from ERA5 are derived from four distinct layers within the top 2.89 m of the subsurface. To ensure the relevance of this data to our study region, which cover a surface less than 10 km², we employ bilinear interpolation on the four nearest nodes around the average location of the seismic array. This process yields representative soil moisture values for each layer. Subsequently, we construct a weighted average of the soil moisture across these layers based on the thickness of each layer to derive a representative soil moisture value for our analysis (Zhong et al. 2023). We identify obvious seasonal and long-term patterns in the soil moisture variations (Fig. 5a). The soil moisture data exhibits a strong
correlation with the observed seismic velocity changes (Fig. 5b), thereby demonstrating significant relevance between variations in water content and shallow seismic velocity changes. We note better consistency between soil moisture in the high-altitude areas and our observed shallow velocity changes (Fig. S17, Supporting Information). This suggests that the shallow water content is significantly influenced by hydrological process in the Kunlun Mountains, with flow systems varying in different flow paths and sources. Solely relying on soil moisture data from the vicinity of the stations does not accurately reflect the variations in water content within the shallow layer. Employing more precise interpolation methods (Kormos et al. 2018; Casellas et al. 2020) or hydrological models (Razandi et al. 2015; Condon et al. 2020) may improve the accuracy of evaluating shallow water content in fine local scale around the seismic array. Furthermore, incorporating multiple types of data, such as direct sampling (Sun et al. 2024; Zou et al. 2024) and remote sensing (Xue et al. 2009; Bao et al. 2010), would be beneficial for future studies in the vicinity of the station areas. Our discussions suggest that water content variations, originating from the runoff from Kunlun Mountains, are the predominant factors influencing seismic velocity change in the shallow layer.

In addition to the observed seasonal variations, a subtle yet significant long-term decrease in seismic velocity (Fig. 2a) can be attributed to an incremental increase in the water content within the shallow subsurface. The analysis of local soil moisture data from 2015 to 2018 reveals a gradual upward trend (Figs 5a, b and Fig. S17, Supporting Information). This long-term increase in water content is primarily driven by two factors: enhanced precipitation and accelerated glacier melting (Hugonnet et al. 2021) in the vicinity of the stations and across the broader Kunlun Mountains area (Figs 5d, e and Fig. S25, Supporting Information). These changes are consistent with the warming trends of the Tibetan Plateau documented in previous studies (Yang et al. 2019; Luo et al. 2020; Zou et al. 2022), demonstrating a clear connection between global climate change and variations in shallow water content on the plateau.

Topography may also contribute to the long-term increase in water content by affecting the distribution in the shallow layer (Murphy et al. 2009; Pellet & Hauck 2017; Gnann et al. 2025). Mapping of the long-term trend of soil moisture (Fig. 5c) highlights the significant influence of topography. In low-altitude areas, particularly around the seismic stations, water redistribution along topographic gradients may amplify the observed increase in seismic wave velocity trends (Fan et al. 2019; Gnann et al. 2025). Additionally, the proximity of the stations to the Nuomuhong river, the major local water system, likely exacerbates this long-term trend, making it more pronounced. It is also noteworthy that in certain high-altitude areas, the trends of soil moisture have decreased, a phenomenon that can be attributed to topographical effects that facilitate water transport. Such spatial variations in moisture trends underscore the complex interplay between topography, climate change and hydrological dynamics in shaping the observed changes in seismic wave velocity.

## 4.2 Seismic velocity changes in the intermediate layer

Seasonal seismic velocity changes observed in the intermediate layer (0.77–1.54 Hz) peak around July and reach their lowest values around January. This layer displays sensitivity to variations roughly spanning from 200 to 500 m in depth, as illustrated in Fig. 3(b). This layer mainly contains the third layer of clay with compact structure

in Fig. 3(c). Unlike the shallow layer, the seismic velocity response of the intermediate layer does not display a negative correlation with precipitation. This lack of correlation indicates a decreased sensitivity of seismic velocity changes to local hydrological processes within this layer, a characteristic attributed to its compact, low-permeability nature.

In contrast, the direct comparison of time-series data has unveiled a significant correlation between changes in seismic velocity and air temperature, albeit with a slight phase lag (Figs 2 and S18, Supporting Information). This noteworthy correlation confirms observations of seasonal seismic velocity changes influenced by temperature fluctuations, a phenomenon associated with thermoelastic strain (Meier et al. 2010; Richter et al. 2014; Hillers et al. 2015). We therefore model the local thermoelastic strain (Tsai 2011) and simulate the synthetic seismic velocity changes $\frac{\Delta v}{v_{\text{syn}}}(t)$ using the following equations:

$$ \frac{\Delta v}{v_{\text{syn}}} (t) = -\frac{m}{\mu} A_t(t) e^{-k_t y} \cos(k_t x) (1 - 2v), \quad (3) $$

where $A_t(t)$ is given by

$$ A_t(t) = \frac{1 + v}{1 - v} k_t \alpha_{th} T_0 \sqrt{\frac{\kappa}{\omega}} e^{-\sqrt{\frac{\omega}{2\kappa}} y_b} \cos \left( \omega t - \sqrt{\frac{\omega}{2\kappa}} y_b - \frac{\pi}{4} \right). \quad (4) $$

$\mu$ and $v$ represent the shear modulus and Poisson's ratio, respectively. We set $v = 0.25$ and $\mu = 6.4 \times 10^9 \text{Pa}$ according to the layer composition in the study area (Wang et al. 2008). $m$ is the third-order Murnaghan constant, which is used to describe the nonlinear elasticity of materials (Murnaghan 1952; Norris 2024). Constraining the exact value of the third-order Murnaghan constant $m$ for geological structures is challenging. Only a few materials have been fully characterized in laboratory experiments, such as steel and sandstone (Hughes & Kelly 1953; Winkler & Liu 1996; Sarkar et al. 2003; Payan et al. 2009). Additionally, previous studies have reported that $m$ can vary significantly, within a range of approximately $\pm 10^{14} \text{Pa}$ (Tsai 2011) due to its strong dependence on rock type, temperature, pressure and other factors (Johnson & Rasolofosaon 1996; D'Angelo et al. 2008; Payan et al. 2009). We provide a rough estimate of $m \approx -3.84 \times 10^{13} \text{Pa}$ to fit the amplitude of observed velocity changes in the intermediate and deep layer. This estimate results in a ratio of $\frac{m}{\mu} = -6000$, which has a larger absolute value than the value of $-2000$ reported by Tsai (2011). This discrepancy may be attributed to the presence of relatively thick and loose sediment layers in the vicinity of our stations (Bao et al. 2017; Xiao et al. 2021), as well as the influence of water content in the shallow layers. Loose sediments and higher water saturation can significantly affect the elastic properties of the medium, leading to larger absolute values of $m$ (Johnson & Rasolofosaon 1996; Winkler & Liu 1996; Winkler & McGowan 2004).

$k_t$ denotes the horizontal wave number of temperature, related to dominant wavelength of local topography, set at $\frac{2\pi}{50 \text{ km}}$ through testing. $\alpha_{th}$ refers to linear thermal expansion coefficient, for which we have selected a value of $10^{-5}/^{\circ}\text{C}$ following Berger (1975). $T_0$ denotes the amplitude of annual temperature variations, estimated using a simple sinusoidal function with $T_0 = 13.35^{\circ}\text{C}$ (Fig. S19, Supporting Information). $\kappa$ is the thermal diffusivity, estimated at $10^{-6} \text{ m}^2 \text{s}^{-1}$ (Berger 1975; Weinert et al. 2020). $\omega$ corresponds to the angular frequency of temperature variations, with a chosen value of $2 \times 10^{-7} \text{s}^{-1}$ to match the main temperature seasonal pattern. $y_b$ represents the thickness of the loose layer that transmits temperature changes without transmitting stress, set based on the delay days between temperature and velocity changes results to $y_b = 0.2 \text{ m}$. $x$
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison of the observed velocity change result and the model fit result in the intermediate layer (0.77—1.54 Hz).</figcaption></figure>

indicates the location of the seismic stations, set to $x = 0$ to ensure that $\cos(k_x x) = 1$. $y$ indicates the depth of model, set to $y = 300$ m as an average estimation.

The synthetic seismic velocity changes ($\frac{\Delta v}{v_{\text{syn}}}(t)$) resulting from thermal strain reasonably align with the seasonal patterns (Fig. 6). However, we observe some differences between the observed velocity changes and the model fit results, which may be attributed to secondary influences originating from hydrological processes, random data loss, and local stress field variations. For example, desiccated and fissured clay layer can amplify effect of hydrological processes (Grisak & Cherry 1975; Ruland et al. 1991).

This study did not fully account for the effects of groundwater flow, which can impact subsurface temperature profiles. Groundwater discharge is generally warmer in winter and cooler in summer compared to surface water (Bilby 1984; Cunjak & Power 1986; Power et al. 1999; Kurylyk et al. 2019), potentially causing a “negative thermal offset” that contributes to the observed phase lag. Furthermore, a long-term decreasing trend in groundwater levels (see Section 4.4) may affect both groundwater temperature and discharge rates (Chu et al. 2008; Kurylyk et al. 2013), potentially contributing to the discrepancies between our model and the observed velocity variations in the intermediate layer since January 2018.

We observe some rapid fluctuations, especially around 2017 January, which may arise from the combined effects of unexpected data loss and short-period events. Unexpected data loss can degrade data quality and influence the final results. Local short-period events, which have not been fully recognized, may also cause changes in subsurface stress and lead to sharp velocity variations (Text S5, Figs S1 and S20–21, Supporting Information).

We hence deduce that the seismic velocity changes in the intermediate layer are predominantly influenced by thermal strain due to its compact nature, resulting in the observed contrasting seasonal pattern. Beyond the seasonal effects, the long-term trend of seismic velocity changes in this layer is fairly flat (Fig. 2b). There is a correlation between long-term trends of temperature records from the local meteorological station (Fig. 1) and seismic velocity changes (Fig. S18, Supporting Information). Given the constraints

of the data length, further exploration is needed to elucidate how the long-term changes in temperature are influencing the seismic velocity in the intermediate layer.

## 4.3 Different thermo-hydrological mechanisms in the shallow and intermediate layers

As discussed above, seismic velocity variations in the shallow layer are primarily controlled by hydrological processes, whereas those in the intermediate layer are dominated by thermal elastic strain. This distinction arises from the markedly different physical and hydrological properties of these two layers, which influence how their seismic velocities respond to temperature fluctuations (Sens-Schönfelder & Larose 2008; Mao et al. 2019; Shen et al. 2024).

In our study area, the stations lie on an upstream slope of an alluvial fan, where the shallow medium primarily consists of unconsolidated deposits. These loose materials are characterized by abundant cracks and high permeability, making the shallow layer particularly susceptible to changes in water content. Such variations not only alter pore pressure but also enhance crack connectivity, thereby significantly modifying the shallow medium and reducing its overall stability (Wu et al. 2010; Yin et al. 2022). Furthermore, freeze-thaw cycles in the seasonal active layer (Fig. 3c) amplify the effects of water content changes, further weakening weaken the soil structure.

We infer that these complex geological processes, spanning high crack density, increased permeability and variations in water content, substantially alter the physical properties of the shallow medium, weakening its strength and reducing its sensitivity to thermal elastic strain. Consequently, water content changes emerge as the primary driver of velocity variations in the shallow layer, as pore pressure effects outweigh potential influences from thermal strain.

In contrast, the intermediate layer consists of Middle Pleistocene glaciofluvial deposits with relatively compact structures and mild cementation, resulting in low permeability. Additionally, higher confining pressures at greater depths limit crack development. Under these conditions, variations in water content have minimal influence on both pore pressure and the elastic modulus. Consequently, thermal elastic strain emerges as the primary driver of
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Seismic velocity changes and water storage from GRACE mascon products. (a) Time-series of liquid water equivalent thickness from CSR mascon solutions, with thick short lines representing the linear trends for distinct time periods (2002–2012, 2013–2018 and 2019–2023). The shaded area represents the study period (2015–2018). (b) Time-series of liquid water equivalent thickness from JPL mascon solutions. (c) Comparison between monthly averages of the seismic velocity changes in the deep layer and the liquid water equivalent thickness, after removing trends and mean values from both JPL and CSR mascon solutions during our study period (2015–2018). The y-axes for liquid water equivalent thickness in panels (a), (b) and (c) are inverted for direct comparison. (d) Comparison of the observed velocity change result and the model fit result in the deep layer (0.25–0.51 Hz). (e) Mapping of the long-term trend of terrestrial water storage from JPL solutions. (f) Mapping of the long-term trend of terrestrial water storage from CSR solutions.</figcaption></figure>

seismic velocity changes in the intermediate layer. Hence, we suggest that distinct geological compositions and processes, including different water content and crack density, lead to contrasting sensitivities of seismic velocity variations to environmental factors in the shallow and intermediate layers (Winkler & McGowan 2004).

## 4.4 Seismic velocity changes in the deep layer

Seismic velocity changes in the deep layer (0.25–0.51 Hz) also exhibit strong seasonal effects, with the highest values occurring around January and the lowest values around July, similar but not identical to the pattern observed in the shallow layer. The sensitivity depth associated with this layer primarily ranges from 1 to 6 km (Fig. 3b). Based on the previous drilling data conducted in the low altitude areas near the seismic stations (Wang et al. 2008), we estimate that the depth of the water table around the stations is approximately 500 m (Fig. 3c) with a simple linear regression (Fig. S22, Supporting Information). Despite the error of the linear regression, it is reasonable to conclude that the main sensitivity range of this layer is located below the water table depth, where the

stress changes resulting from thermal strain have dissipated. The replenishment of local groundwater primarily originates from precipitation and the melting water in the Kunlun Mountains (Wang et al. 2008; Bibi et al. 2019). Given the consistent seasonal variation patterns of these two main recharge sources discussed above (Figs 4c and d) and rapid water circulation (Tan et al. 2021; Zhang et al. 2023), we suggest that the seasonal variations of groundwater play a significant role in the seismic velocity changes in the deep layer.

Due to the harsh natural conditions around study region, direct groundwater observations are challenging. Therefore, we use GRACE mascon products from the Jet Propulsion Laboratory (JPL; Watkins et al. 2015; Wiese et al. 2023) and the Center for Space Research (CSR; Save et al. 2016; Save 2020) to further verify our assumptions. Given the significant contribution of ground water storage to the terrestrial water storage in the vicinity of the study region (Xiang et al. 2016), the GRACE-derived liquid water equivalent thickness can be approximated to the groundwater storage.

Except that there is a temporal observational gap in GRACE data around 2018 (Landerer et al. 2020), the spatial bilinear interpolated liquid water equivalent thickness shows a consistent seasonal
pattern comparable with velocity change in the deep layer (Figs 7a–c). When water storage peaks each summer around July, seismic velocity changes tend to reach their minimum.

Therefore, we suggest that seasonal groundwater variation plays a significant role in the seismic velocity changes observed in the deep layer.

We evaluate pore pressure changes from groundwater changes using the following equation (Tsai 2011):

$$ \frac{\Delta v}{v_{\text{syn}}} (t) = \frac{m}{\mu} A_h(t) e^{-k_h y} \cos(k_h x) (1 - 2\nu). \quad (5) $$

$A_h(t)$ is given by

$$ A_h(t) = \frac{(1 + \nu)(1 - 2\nu)}{E(1 - \nu)} k_h \alpha p_0 \times \sqrt{\frac{\kappa_{hy}}{\omega}} e^{\frac{\pi}{4} - \frac{1}{2}\cot^{-1}\kappa_{hy}} \cos\left(\omega t - \frac{1}{2}\cot^{-1}\kappa_{hy}\right), \quad (6) $$

where $E$ is Young's modulus. $\alpha$ is the Biot-Willis coefficient. $\kappa_{hy}$ is the hydraulic diffusivity. Following the value choices in Tsai (2011), we set $E = 1.6 \times 10^{10}$ Pa, $\alpha = 0.7$, and $\kappa_{hy} = 0.4 \text{ m s}^{-2}$. $p_0$ represents the annual amplitude of pore pressure variations from the groundwater. $k_h$ denotes the horizontal wave number of groundwater variations and we set it to the same value as $k_t$, which is $\frac{2\pi}{50 \text{ km}}$.

Groundwater variations tend to increase with altitude in many high mountain areas (Somers & McKenzie 2020), consistent with previous geological prospecting findings in the low-altitude areas (Wang et al. 2008; Zhang 2013) close to the study region. Drilling data at an altitude of approximately 2900 m reveals groundwater fluctuations up to 4.4 m. Considering that seismic stations are located at around 3300 m and are in close proximity to rivers, we have estimated that the annual groundwater fluctuates around the seismic stations approximate 5 m, equivalent to $p_0 = 2.5 \times 10^4$ Pa. Assuming that the local groundwater changes follow a similar seasonal pattern to that of temperature (Figs 2 and 4), we fit a simple sinusoidal function to temperature to model the annual groundwater changes. Since our model primarily focuses on fitting smooth seasonal variations, it is inevitable that some local rapid changes between the actual groundwater variations and the model cannot be simulated.

We model the synthetic velocity changes from the direct elastic loading effects using the following equation (Tsai 2011):

$$ \frac{\Delta v}{v_{\text{syn}}} (t) = \frac{m}{\mu} A_e(t) e^{-k_h y} \cos(k_h x) (1 - 2\nu), \quad (7) $$

where $A_e(t)$ is defined as following:

$$ A_e(t) = \frac{(1 + \nu)\phi p_0}{E} \cos \omega t. \quad (8) $$

$\phi$ denotes the porosity, which we have set to 0.2 based on the results of previous groundwater model in the Qaidam Basin (Hao et al. 2018; Xiao 2018).

We consider the effects of thermal elastic strain (eqs 3 and 4), pore pressure changes (eqs 5 and 6), and direct elastic loading (eqs 7 and 8) to model the seasonal variations in seismic velocity at a depth of 3500 m, which is used as an average estimation (Fig. S26, Supporting Information). The model demonstrates a robust fit with seasonal patterns (Fig. 7d). This supports the hypothesis that changes in seismic velocity within the deep layer are predominantly influenced by changes in groundwater levels. Discrepancies in the annual patterns between the model and seismic observations may arise from inaccuracies in the model parameters we have chosen and from a lack of consideration of local rapid changes. Additionally,

the absence of rapidly changing information can also affect the accuracy of the model fit.

Moreover, the slight upward trend observed in velocity changes within the deep layer (Fig. 2c) suggests a potential decline in local groundwater levels. This trend is evident in both the GRACE Mascon products (Figs 7a and b) during our study period. There is a significant increasing trend during 2002–2012. However, a turning point occurred around 2013 and the groundwater levels around seismic stations reveals a decreasing trend after 2013, encompassing our study period (Fig. 7). It is important to note that the y-axes for liquid water equivalent thickness in Figs 7(a)–(c), are inverted for direct comparison. This decrease in groundwater levels during 2013–2018, common across the study region (Figs 7e and f), aligns with the long-term trend observed in the Qaidam Basin (Fig. S23, Supporting Information), as reported in previous studies (Bibi et al. 2019; Zou et al. 2022). Climate may have complex effects on changes in groundwater level in the study region. During the summertime (average of June, July and August), the North Atlantic Oscillation index exhibits a strong, abrupt change around 2013 (Fig. S24, Supporting Information). This change likely has significant influences on atmospheric circulation and affects summer precipitation on the Tibetan Plateau (Liu et al. 2021), potentially contributing to the observed turning point in long-term groundwater trends around 2013. Previous studies also suggest that the El Niño-Southern Oscillation may greatly influence groundwater levels on the Tibetan Plateau (Gao et al. 2022; Xie et al. 2022; Chen et al. 2024). Additionally, the melting of ice and snow due to increasing temperatures will temporarily increase groundwater storage. However, the long-term water loss resulting from glacier retreat (Fig. S25, Supporting Information) may not be fully compensated by the increase in precipitation (Hugonnet et al. 2021; Wang et al. 2023; Zhang et al. 2023). Moreover, human activities like agriculture and industry, which heavily rely on local groundwater due to the alpine arid climate, significantly influence the long-term groundwater variations (Li et al. 2019; Zhao et al. 2021). We suggest that the observed long-term decrease in groundwater storage during 2013–2018 is responsible for the subtle increasing trend in seismic velocities within the deep layer over our study period. Considering the limited areas studied and the uncertainty in the model fit, further research is needed to gain insight into the complex hydrological processes under climate change at different depths on the Tibetan Plateau.

## 4.5 Illustration of subsurface processes

Following a thorough evaluation in the preceding discussions, we outline the overall subsurface physical processes by Fig. 8. During the summer season, there is a notable increase in precipitation and melting water within the Kunlun Mountains. These climatic phenomena exert a significant influence on the groundwater levels within the deeper layer and the water content in the shallow layer, leading to an increase in both. Consequently, the heightened water content leads to increased pore pressures and the application of direct elastic hydrologic loads, ultimately causing a reduction in seismic velocity.

The impact of hydrological processes is mitigated by the presence of Middle Pleistocene glaciofluvial deposits in the intermediate layer. These deposits are characterized by their compactness and low permeability. The primary factor affecting this layer is thermal strain induced by the rising temperatures, which correspondingly results in an increase in seismic velocity in summer.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Illustration of the subsurface processes. (a) Physical processes during summer. (b) Physical processes during winter.</figcaption></figure>

On the contrary, during the winter season, there is a reduction in precipitation and melting water, leading to an increase in seismic velocity in both the shallow and deep strata. Simultaneously, temperature fluctuations play a key role in diminishing seismic velocity within the intermediate layer due to thermal strain.

Furthermore, the observed long-term decrease in seismic velocity (Fig. 2a) in the shallow layer suggests an increase in water content, while the observed increase trend in the deep layer (Fig. 2c) corresponds to a reduction in groundwater levels during the study period. This observed discrepancy in long-term trends can be attributed to variations in hydrological factors across multiple depth strata, such as flow pathways, cycling depths, discharge locations and geological structures, as illustrated in Fig. 8.

## 5 CONCLUSIONS

We use unique high-altitude seismic recordings to monitor seismic velocity changes over time and unveil alterations in the physical properties of the subsurface medium within the Nuomuhong region of the northeastern Tibetan Plateau. Through the integration of local meteorological data, ERA5 climate re-analysis data set and GRACE products, we identify the dominant influence of hydrological dynamic processes, primarily driven by precipitation and melting water from the Kunlun Mountains, on the observed seasonal velocity variations, except within the intermediate layer. The intermediate layer exhibits distinct velocity changes attributed to thermal strain induced by the presence of compact and poorly permeable Middle Pleistocene glaciofluvial deposits at these depths. The contrasting

trends observed in the shallow and deep layers suggest a decoupling in water circulation processes across different depths.

Our findings highlight the potential of using ambient noise as an effective means for continuously monitoring hydrothermal dynamics in frozen ground environments. In the face of intricate environmental conditions, seismic interferometry emerges as a robust tool with a wide range of applications, from ecosystem conservation and water resource management to enhancing our understanding of the impacts of climate change on subsurface processes.

## ACKNOWLEDGMENTS

We gratefully acknowledge the financial support of the National Key Research and Development Program of China (Grant No. 2021YFC3000704), and the National Natural Science Foundation of China (Grant No. 42488301). The authors want to gratefully acknowledge the Editor, Dr Anne Obermann, as well as Dr Xiaozhuo Wei and an anonymous reviewer for their valuable comments and suggestions, which helped to significantly improve this manuscript.

## SUPPORTING INFORMATION

Supplementary data are available at GJIRAS online.

**Text S1:** Details of spectral whitening.

**Text S2:** Singular value decomposition.

**Text S3:** Detailed discussion about SVD.

**Text S4:** Results with MWCS and stretching.
**Text S5:** Possible reasons for transient seismic velocity changes in the intermediate layer.

**Figure S1:** The data coverage of the seismic data used in this study. (a) The daily percentage of usable data. (b) The daily number of stations with recorded data.

**Figure S2:** Long-term meteorological records. The figure presents the extended meteorological records collected from the meteorological station (see Figure 1 in the main text for the station's location). The dataset comprises daily measurements, with the exception of precipitation, which is averaged annually. The shaded areas indicate the time range covered by the seismic records used in this article. The linear long-term trends of the meteorological are depicted by straight lines in each subplot. (a) Temperature variations; (b) Evaporation rates; (c) Annual averaged precipitation.

**Figure S3:** The distribution of average singular values from the decomposition of auto-correlation matrices corresponding to various stations.

**Figure S4:** Comparison of auto-correlation functions before and after SVD. The figure illustrates a comparison of the auto-correlation functions of station NMH7, both prior to and after SVD. Missing data is represented in dark blue. (a) and (b) depict the coda waves from the negative and positive portions of the original auto-correlation functions for station NMH7 within the time window of 20 s to 100 s; (c) and (d) show the auto-correlation functions after SVD for station NMH7 within the same time window.

**Figure S5:** Comparison of NMH7's velocity change results before and after SVD in the frequency range of 0.77—1.54 Hz. The figure presents a comparison of NMH7's results before and after SVD, focusing on the frequency range of 0.77 Hz to 1.54 Hz. The light green lines depict the results before SVD, while the purple lines represent the results after SVD. (a) Relative velocity variation; (b) Error in velocity variation; (c) Average correlation coefficients between daily correlation functions and the reference correlation function.

**Figure S6:** Comparison of auto-correlations from station NMH7 in the 0.77-1.54 Hz frequency band after filtering and SVD. (a) Daily auto-correlation functions (right side) after filtering. (b) Daily auto-correlation functions (right side) after SVD. (c) Stacked auto-correlation functions (right side) after filtering, with standard deviation shown in light green. (d) Stacked auto-correlation functions (right side) after SVD, with standard deviation shown in light green.

**Figure S7:** The travel times differences and SNR of auto-correlations from station NMH7 in the frequency band of 0.77 Hz-1.54 Hz across lapse time. (a) The SNR of stacked CCF after filtering along the lapse time. (b) The SNR of stacked CCF after SVD along the lapse time. (c) The travel time differences between the auto-correlations in the August 20, 2016, and the stacked auto-correlations (reference auto-correlations after filtering). (d) The travel time differences between the auto-correlations in the August 20, 2016, and the stacked auto-correlations (reference autocorrelations after SVD).

**Figure S8:** Comparison of NMH7's velocity change results with filtering and SVD in the frequency range of 0.77—1.54 Hz. The light green lines represent the results with filtering, while the purple lines represent the results with SVD. (a) Relative velocity variation; (b) Error in velocity variation; (c) Average correlation coefficients between daily correlation functions and the reference correlation function.

**Figure S9:** The comparison of different choices of the parameter N in the weighted sliding average method in different frequency bands. The blue, red, orange, purple, green, light blue lines represent the results with N=61, 91, 181, 361, 721 and the full dataset,

respectively. (a) 5—10 Hz; (b) 0.77—1.54 Hz; (c) 0.25—0.51 Hz.

**Figure S10:** Comparison of average results using a sliding average with a time window of L = 61 days versus a weighted sliding average with L = 31 days in three frequency bands. The orange and blue curves depict the results obtained from the sliding average and the weighted sliding average, respectively. (a) 5—10 Hz; (b) 0.77—1.54 Hz; (c) 0.25—0.51 Hz.

**Figure S11:** Comparison of relative seismic velocity change results using three different methods at multiple stations across varying frequency bands. The red, purple, and orange solid curves represent the results from the wavelet cross-spectrum, moving window cross-spectrum (MWCS), and stretching methods, respectively. (a) Frequency band 5–10 Hz; (b) Frequency band 0.77–1.54 Hz; (c) Frequency band 0.25–0.51 Hz.

**Figure S12:** Comparison of relative seismic velocity change results using three different methods at multiple stations across varying frequency bands. The red, purple, and orange solid curves represent the results from the wavelet cross-spectrum, moving window cross-spectrum (MWCS), and stretching methods, respectively. The red, purple, and orange dashed curves represent the long-term trends of the results from the wavelet cross-spectrum, MWCS, and stretching methods, respectively. (a) Frequency band 5–10 Hz; (b) Frequency band 0.77–1.54 Hz; (c) Frequency band 0.25–0.51 Hz.

**Figure S13:** Time series of seismic observations including data from station NMH1 and meteorological recordings. (a), (b), and (c) represent the seismic velocity changes in the frequency ranges of 5—10 Hz, 0.77—1.54 Hz, and 0.25—0.51 Hz, respectively. These changes are averaged over results with standard deviation in light blue shade from all seismic stations including station NMH1. The orange dashed lines represent the linear trend of the velocity change results for each frequency band. (d) Air temperature smoothed with a 31-day sliding window; (e) Daily precipitation; (f) Atmospheric pressure smoothed with a 31-day sliding window.

**Figure S14:** The noise spectrum of station NMH6. (a) Daily PSD of station NMH6. (b) The average PSD of station NMH6.

**Figure S15:** Daily variation of the seismic records of station NMH6. (a) Daily root mean square (RMS) of station NMH6. (b) Daily standard deviation (STD) of station NMH6.

**Figure S16:** The spatial distribution and altitudinal variation of annual average of net precipitation (precipitation minus evaporation) and snow melting water. (a) The annual average of net precipitation (precipitation minus evaporation) around the stations. The red line represents the transect for the profiles in (c), (d), and (e). The orange triangle represents the location of the seismic stations. (b) The annual average of melting water around the stations. (c) The altitude profile along the transect AB (referenced in (a) and (b)). The red triangle represents the location of stations. (d) The profile of net precipitation (precipitation minus evaporation) along the transect AB. (e) The melting water profile along the transect AB.

**Figure S17:** The comparison between shallow velocity changes and soil moisture variations in the Kunlun Mountains. (a) The distribution of correlation coefficients between soil moisture and observed shallow seismic velocity changes. The purple square (95.8°E–96.39°E, 35.1°N–36.2°N) represents the area averaged in (b) and (c). (b) The average variations of soil moisture within the area marked by purple squares in (a). The blue shaded area represents our study period. The green line represents the linear trend in soil moisture during our study period. Note that the direction of the y-axis is inverted. (c) The comparison between shallow velocity changes and soil moisture data within the area highlighted by the
purple square in (a). The blue line represents the observed shallow velocity changes, and the purple line represents the soil moisture data.

**Figure S18:** The comparison between intermediate velocity changes and smoothed temperature, focusing on their long-term trends. (a) The normal comparison using their actual values. The solid purple line represents the seismic velocity changes in the intermediate layers. The purple dot-dashed line represents the linear trend of seismic velocity changes in the intermediate layers. The solid orange line represents the smoothed temperature. The orange dashed line represents the linear trend of the smoothed temperature. (b) The comparison after normalization, where both the seismic velocity changes and temperature data are scaled to a uniform range for a more direct comparison of their trends. The representation of the lines remains consistent with that in (a).

**Figure S19:** Comparison of synthetic and observed temperatures. The orange line represents the observed temperature records from the meteorological station during the study period. The blue line illustrates the synthetic temperature used in modeling thermoelastic effects as described in Section 4.2 of the main text.

**Figure S20:** The average seismic velocity changes in the frequency band of 0.77–1.54 Hz under different conditions. (a) Results retaining only days with full data. (b) Average results after interpolating all missing days’ data for the stations.

**Figure S21:** Local events recorded at stations NMH1–NMH3 on January 1st. (a)–(c) Possible local “repeating events” recorded at stations NMH1–NMH3.

**Figure S22:** Partial historical groundwater table levels derived from previous drilling activities in the downstream regions relative to our monitoring stations. The red circles represent the water depth measurements from these previous drillings. The green line represents the linear regression line fitted to the groundwater table data. The purple shaded area maps the altitude of our station areas. The drilling data are obtained from Wang et al. (2008).

**Figure S23:** Groundwater storage changes in the Qaidam Basin, as reported by Zou et al. (2022). The blue shaded area indicates the time range covered by the seismic records used in this study.

**Figure S24:** Analysis of temporal variations in GRACE mascon products from JPL and CSR at the average locations of the stations, alongside the summer NAO index. (a) Time series of liquid water equivalent thickness derived from CSR mascon solutions, with green lines representing the linear trends for distinct periods (2002–2012,2013–2018 and 2019–2023). The blue shaded area represents the study period (2015-2018). (b) Time series of liquid water equivalent thickness derived from JPL mascon solutions. Note that the y-axes for liquid water equivalent thickness in (a) and (b), are inverted for direct comparison. (c) The time series of the summer NAO index is represented by the orange line.

**Figure S25:** Average annual glacier mass change rate within the area bounded by 95.5°E – 97°E and 35°N – 36.5°N. The data are provided by Hugonnet et al. (2021) with a resolution 0.5° × 0.5°.

**Figure S26:** Models of different effects contributing to seismic velocity changes in the deep layer. The blue line represents the seismic velocity changes in the deep layer, and the black line represents the modeled effects. (a) Model of the thermoelastic strain effect. (b) Model of the poroelastic (pore pressure) effect. (c) Model of the direct elastic hydrologic loading effect. (d) Final model combining all effects.

**Figure S27:** Relative velocity change results after randomly removing twenty days of data 50 times in different frequency bands. The solid line represents the result without random data removal, while the shaded area represents the standard deviation range of the

results with random data removal. (a) 5–10 Hz; (b) 0.77–1.54 Hz; (c) 0.25–0.51 Hz.

**Figure S28:** Comparison of seismic velocity change results using different time windows. The blue line represents the results for the window 20 s–60 s (the window used in the main text). The orange line represents the results for the window 25 s–65 s. The yellow line represents the results for the window 30 s–70 s. The purple line represents the results for the window 35 s–75 s. The green line represents the results for the window 40 s–80 s. (a) Results in the frequency band of 5–10 Hz. (b) Results in the frequency band of 0.77–1.54 Hz. (c) Results in the frequency band of 0.25–0.51 Hz.

**Figure S29:** Average seasonal seismic velocity variations across seven frequency bands (5–10 Hz; 3–5 Hz; 1.54–3 Hz; 0.77–1.54 Hz; 0.51–0.77 Hz; 0.25–0.51 Hz; 0.2–0.25 Hz). The values are normalized by each band’s respective maximum absolute value to account for the large differences in amplitude across different layers.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## DATA AVAILABILITY

The daily autocorrelation functions of nine stations and average relatively seismic velocity change results (Yu et al. 2025) used in this article can be found at https://doi.org/10.5281/zenodo.14690410. The river location data provided by the National Catalogue Service for Geographic Information can be accessed at https://www.webmap.cn/commres.do?method= result100W. The code for the wavelet cross-spectrum method can be accessed at https://github.com/shujuanmao/dt-wavelet. The meteorological data provided by the National Cryosphere Desert Data Center (Ye 2022) is available at https://doi.org/10.12072/ncdc.Hydro.db2446.2022. ERA5-Land monthly averaged data can be found at https://doi.org/10.24381/cds.68d2bb30. CSR GRACE/GRACE-FO RL06.2 Mascon Solutions (Save et al. 2016; Save 2020) can be found at https://www2.csr.utexas.edu/grace/RL06_mascons.html. JPL GRACE and GRACE-FO Mascon Ocean, Ice, and Hydrology Equivalent Water Height Coastal Resolution Improvement Filtered Release 06.1 Version 03 (Watkins et al. 2015; Wiese et al. 2023) can be found at http://grace.jpl.nasa.gov/data/get-data/jpl_global_mascons/. Surface wave sensitivity kernels are computed by the Computer Programs in Seismology (Herrmann 2013), can be accessed at https://www.eas.slu.edu/eqc/ComputerProgramsSeismology/index.html. The North Atlantic Oscillation index comes from the National Oceanic and Atmospheric Administration https://www.ncei.noaa.gov/access/monitoring/nao/. Glacier mass changes from Hugonnet et al. (2021) were provided by Theia https://doi.org/10.6096/13. Most figures in this study are plotted by Matlab R2021a (http://www.mathworks.com/), matplotlib (Hunter 2007) and GMT (Wessel et al. 2019).

## REFERENCES

* Acarel, D., Bulut, F., Bohnhoff, M. & Kartal, R., 2014. Coseismic velocity change associated with the 2011 van earthquake (m7.1): crustal response to a major event, *Geophys. Res. Lett.*, **41**(13), 4519–4526.

* Bao, J., Wang, Y., Song, C., Feng, Y., Hu, C., Zhong, S. & Yang, J., 2017. Cenozoic sediment flux in the qaidam basin, northern tibetan plateau, and implications with regional tectonics and climate, *Global Planet. Change*, **155**, 56–69.
Bao, Q., Liu, Y., Shi, J. & Wu, G., 2010. Comparisons of soil moisture datasets over the tibetan plateau and application to the simulation of asia summer monsoon onset, *Adv. Atmos. Sci.*, **27**, 303–314.

Barajas, A., Margerin, L. & Campillo, M., 2022. Coupled body and surface wave sensitivity kernels for coda-wave interferometry in a three-dimensional scalar scattering medium, *J. Geophys. Int.*, **230**(2), 1013–1029.

Ben-Zion, Y. & Leary, P., 1986. Thermoelastic strain in a half-space covered by unconsolidated material, *Bull. seism. Soc. Am.*, **76**(5), 1447–1460.

Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M. & Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, *J. Geophys. Int.*, **169**(3), 1239–1260.

Berger, J., 1975. A note on thermoelastic strains and tilts, *J. Geophys. Res.*, **80**(2), 274–277.

Bibi, S., Wang, L., Li, X., Zhang, X. & Chen, D., 2019. Response of groundwater storage and recharge in the qaidam basin (tibetan plateau) to climate variations from 2002 to 2016, *J. Geophys. Res.: Atmos.*, **124**(17–18), 9918–9934.

Bilby, R.E., 1984. Characteristics and frequency of cool-water areas in a western washington stream, *J. Freshw. Ecol.*, **2**(6), 593–602.

Biskaborn, B.K., Lanckman, J.-P., Lantuit, H., Elger, K., Streletskiy, D.A., Cable, W.L. & Romanovsky, V.E., 2015. The new database of the global terrestrial network for permafrost (gtn-p), *Earth Syst. Sci. Data*, **7**(2), 245–259.

Biskaborn, B.K. et al., 2019. Permafrost is warming at a global scale, *Nat. Commun.*, **10**(1), 264, doi:10.1038/s41467-018-08240-4.

Brenguier, F., Campillo, M., Hadzioioannou, C., Shapiro, N.M., Nadeau, R.M. & Larose, E., 2008a. Postseismic relaxation along the san andreas fault at parkfield from continuous seismological observations, *Science*, **321**(5895), 1478–1481.

Brenguier, F., Shapiro, N.M., Campillo, M., Ferrazzini, V., Duputel, Z., Coutant, O. & Nercessian, A., 2008b. Towards forecasting volcanic eruptions using seismic noise, *Nat. Geosci.*, **1**(2), 126–130.

Brenguier, F., Campillo, M., Takeda, T., Aoki, Y., Shapiro, N., Briand, X., Emoto, K. & Miyake, H., 2014. Mapping pressurized volcanic fluids from induced crustal seismic velocity drops, *Science*, **345**(6192), 80–82.

Casellas, E., Bech, J., Veciana, R., Miró, J.R., Sairouni, A. & Pineda, N., 2020. A meteorological analysis interpolation scheme for high spatial-temporal resolution in complex terrain, *Atmos. Res.*, **246**, 105–103.

Chen, B. & Liu, Z., 2016. Global water vapor variability and trend from the latest 36 year (1979 to 2014) data of ECMWF and NCEP reanalyses, radiosonde, GP, and microwave satellite, *J. Geophys. Res.: Atmos.*, **121**(19), 11442–11462.

Chen, T., Pan, Y., Ding, H., Jiao, J. & He, M., 2024. Investigating terrestrial water storage variation and its response to climate in southeastern Tibetan plateau inferred through space geodetic observations, *J. Hydrol.*, **640**, 131742, doi:10.1016/j.jhydrol.2024.131742.

Chu, C., Jones, N.E., Mandrak, N.E., Piggott, A.R. & Minns, C.K., 2008. The influence of air temperature, groundwater discharge, and climate change on the thermal diversity of stream fishes in southern Ontario watersheds, *Can. J. Fish. Aquat. Sci.*, **65**(2), 297–308.

Clarke, D., Zaccarelli, L., Shapiro, N.M. & Brenguier, F., 2011. Assessment of resolution and accuracy of the Moving Window Cross Spectral technique for monitoring crustal temporal variations using ambient seismic noise, *J. Geophys. Int.*, **186**(2), 867–882.

Clements, T. & Denolle, M.A., 2018. Tracking groundwater levels using the ambient seismic field, *Geophys. Res. Lett.*, **45**(13), 6459–6465.

Colombi, A., Chaput, J., Brenguier, F., Hillers, G., Roux, P. & Campillo, M., 2014. On the temporal stability of the coda of ambient noise correlations, *C. R. Geosci.*, **346**(11), 307–316.

Condon, L.E., Atchley, A.L. & Maxwell, R.M., 2020. Evapotranspiration depletes groundwater under warming over the contiguous united states, *Nat. Commun.*, **11**(1), 873, doi:10.1038/s41467-020-14688-0.

Cunjak, R.A. & Power, G., 1986. Winter habitat utilization by stream resident brook trout (salvelinus fontinalis) and brown trout (salmo trutta), *Can. J. Fish. Aquat. Sci.*, **43**(10), 1970–1981.

D'Angelo, R., Winkler, K. & Johnson, D., 2008. Three wave mixing test of hyperelasticity in highly nonlinear solids: sedimentary rocks, *J. acoust. Soc. Am.*, **123**(2), 622–639.

Delouche, E. & Stehly, L., 2023. Seasonal seismic velocity variations measured using seismic noise autocorrelations to monitor the dynamic of aquifers in greece, *J. Geophys. Res. Solid Earth*, **128**(12), e2023JB026759, doi:10.1029/2023JB026759.

Fan, Y. et al., 2019. Hillslope hydrology in global change research and earth system modeling, *Water Resour. Res.*, **55**(2), 1737–1772.

Froment, B., Campillo, M., Roux, P., Gouedard, P., Verdel, A. & Weaver, R.L., 2010. Estimation of the effect of nonisotropically distributed energy on the apparent arrival time in correlations, *Geophysics*, **75**(5), SA85-SA93.

Gao, H., Han, C., Chen, R., Feng, Z., Wang, K., Fenicia, F. & Savenije, H., 2022. Frozen soil hydrological modeling for a mountainous catchment northeast of the qinghai–tibet plateau, *Hydrol. Earth Syst. Sci.*, **26**(15), 4187–4208.

Gnann, S., Baldwin, J.W., Cuthbert, M.O., Gleeson, T., Schwanghart, W. & Wagener, T., 2025. The influence of topography on the global terrestrial water cycle, *Rev. Geophys.*, **63**(1), e2023RG000810, doi:10.1029/2023RG000810.

Green, T.R., Taniguchi, M., Kooi, H., Gurdak, J.J., Allen, D.M., Hiscock, K.M., Treidel, H. & Aureli, A., 2011. Beneath the surface of global change: impacts of climate change on groundwater, *J. Hydrol.*, **405**(3–4), 532–560.

Grisak, G. & Cherry, J., 1975. Hydrologic characteristics and response of fractured till and clay confining a shallow aquifer, *Can. Geotech. J.*, **12**(1), 23–43.

Hadzioioannou, C., Larose, E., Coutant, O., Roux, P. & Campillo, M., 2009. Stability of monitoring weak changes in multiply scattering media with ambient noise correlation: laboratory experiments, *J. acoust. Soc. Am.*, **125**(6), 3688–3695.

Hao, Q., Lu, C., Zhu, Y., Xiao, Y. & Gu, X., 2018. Numerical investigation into the evolution of groundwater flow and solute transport in the eastern Qaidam basin since the last glacial period, *Geofluids*, **2018**, 9260604, doi:10.1155/2018/9260604.

Herrmann, R.B., 2013. Computer programs in seismology: an evolving tool for instruction and research, *Seismol. Res. Lett.*, **84**(6), 1081–1088.

Hillers, G., Campillo, M. & Ma, K.F., 2014. Seismic velocity variations at TCDP are controlled by MJO driven precipitation pattern and high fluid discharge properties, *Earth planet. Sci. Lett.*, **391**, 121–127.

Hillers, G., Ben-Zion, Y., Campillo, M. & Zigone, D., 2015. Seasonal variations of seismic velocities in the San Jacinto fault area observed with ambient seismic noise, *J. Geophys. Int.*, **202**(2), 920–932.

Hirose, T. & Ueda, H., 2023. Continuous monitoring of instrumental clock errors at 50 volcanoes in japan based on seismic interferometry, *Earth Planets Space*, **75**(1), 50, doi:10.1186/s40623-023-01798-4.

Hodgkins, S.B., Tfaily, M.M., McCalley, C.K., Logan, T.A., Crill, P.M., Saleska, S.R., Rich, V.I. & Chanton, J.P., 2014. Changes in peat chemistry associated with permafrost thaw increase greenhouse gas production, *Proc. Natl. Acad. Sci. USA*, **111**(16), 5819–5824.

Hughes, D.S. & Kelly, J.L., 1953. Second-order elastic deformation of solids, *Phys. Rev.*, **92**(5), 1145, doi:10.1103/PhysRev.92.1145.

Hugonnet, R. et al., 2021. Accelerated global glacier mass loss in the early twenty-first century, *Nature*, **592**(7856), 726–731.

Hunter, J.D., 2007. Matplotlib: a 2-D graphics environment, *Comput. Sci. Eng.*, **9**(3), 90–95.

James, S.R., Knox, H.A., Abbott, R.E., Panning, M.P. & Screaton, E.J., 2019. Insights into permafrost and seasonal active-layer dynamics from ambient seismic noise monitoring, *J. Geophys. Res.: Earth Surface*, **124**(7), 1798–1816.

Johnson, P. & Rasolofosaon, P., 1996. Nonlinear elasticity and stress-induced anisotropy in rock, *J. Geophys. Res. Solid Earth*, **101**(B2), 3113–3124.

Jorgenson, M.T., Racine, C.H., Walters, J.C. & Osterkamp, T.E., 2001. Permafrost degradation and ecological changes associated with a warming climate in central alaska, *Clim. Change*, **48**(4), 551–579.
Kormos, P.R. et al., 2018. 31 years of hourly spatially distributed air temperature, humidity, and precipitation amount and phase from reynolds critical zone observatory, *Earth Syst. Sci. Data*, **10**(2), 1197–1205.

Kurylyk, B.L., Bourque, C.-A. & MacQuarrie, K.T., 2013. Potential surface temperature and shallow groundwater temperature response to climate change: an example from a small forested catchment in east-central New Brunswick (Canada), *Hydrol. Earth Syst. Sci.*, **17**(7), 2701–2716.

Kurylyk, B.L., Irvine, D.J. & Bense, V.F., 2019. Theory, tools, and multidisciplinary applications for tracing groundwater fluxes from temperature profiles, *Wiley Interdiscip. Rev.: Water*, **6**(1), e1329, doi:10.1002/wat2.139.

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.

Larose, E., De Rosny, J., Margerin, L., Anache, D., Gouedard, P., Campillo, M. & Van Tiggelen, B., 2006. Observation of multiple scattering of khz vibrations in a concrete structure and application to monitoring weak changes, *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.*, **73**(1), 016609, doi:10.1103/PhysRevE.73.016609.

Laske, G., Masters, G., Ma, Z. & Pasyanos, M., 2013. Update on crust 1.0—a 1-degree global model of earth’s crust, *Geophys. Res. Abstr.*, **15**(15), EGU2013–2658.

Li, H., Mao, D., Li, X., Wang, Z. & Wang, C., 2019. Monitoring 40-year lake area changes of the Qaidam basin, Tibetan plateau, using landsat time series, *Remote Sensing*, **11**(3), 343, doi:10.3390/rs11030343.

Li, X., Long, D., Scanlon, B.R., Mann, M.E., Li, X., Tian, F., Sun, Z. & Wang, G., 2022. Climate change threatens terrestrial water storage over the tibetan plateau, *Nat. Clim. Change*, **12**(9), 801–807.

Lindner, F., Wassermann, J. & Igel, H., 2021. Seasonal freeze–thaw cycles and permafrost degradation on Mt. Zugspitze (German/Austrian Alps) revealed by single-station seismic monitoring, *Geophys. Res. Lett.*, **48**(18), e2021GL094659, doi:10.1029/2021GL094659.

Liu, Y., Chen, H. & Hu, X., 2021. The unstable relationship between the precipitation dipole pattern in the tibetan plateau and summer nao, *Geophys. Res. Lett.*, **48**(4), e2020GL091941, doi:10.1029/2020GL091941.

Lobkis, O.I. & Weaver, R.L., 2003. Coda-wave interferometry in finite solids: recovery of P-to-S conversion rates in an elastodynamic billiard, *Phys. Rev. Lett.*, **90**(25), 254302.

Luo, S., Wang, J., Pomeroy, J.W. & Lyu, S., 2020. Freeze–thaw changes of seasonally frozen ground on the Tibetan Plateau from 1960 to 2014, *J. Clim.*, **33**(21), 9427–9446.

Mao, S., Campillo, M., van Der Hilst, R.D., Brenguier, F., Stehly, L. & Hillers, G., 2019. High temporal resolution monitoring of small variations in crustal strain by dense seismic arrays, *Geophys. Res. Lett.*, **46**(1), 128–137.

Mao, S., Mordret, A., Campillo, M., Fang, H. & van der Hilst, R.D., 2020. On the measurement of seismic traveltime changes in the time-frequency domain with wavelet cross-spectrum analysis, *J. Geophys. Int.*, **221**(1), 550–568.

Mao, S., Lecointre, A., van der Hilst, R.D. & Campillo, M., 2022. Space-time monitoring of groundwater fluctuations with passive seismic interferometry, *Nat. Commun.*, **13**(1), 4643, doi:10.1038/s41467-022-32194-3.

Meier, U., Shapiro, N.M. & Brenguier, F., 2010. Detecting seasonal variations in seismic velocities within Los Angeles basin from correlations of ambient seismic noise, *J. Geophys. Int.*, **181**(2), 985–996.

Melo, G., Malcolm, A., Mikesell, D. & van Wijk, K., 2013. Using SVD for improved interferometric green’s function retrieval, *J. Geophys. Int.*, **194**(3), 1596–1612.

Moreau, L., Stehly, L., Boué, P., Lu, Y., Larose, E. & Campillo, M., 2017. Improving ambient noise correlation functions with an SVD-based Wiener filter, *J. Geophys. Int.*, **211**(1), 418–426.

Muñoz Sabater, J., 2019. ERA5-land monthly averaged data from 1950 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS), doi:10.24381/eds.68d2bb30.

Murnaghan, F., 1952. Finite deformation of an elastic solid, *Bull. Am. Math. Soc*, **58**, 577–579.

Murphy, P., Ogilvie, J. & Arp, P., 2009. Topographic modelling of soil moisture conditions: a comparison and verification of two models, *Eur. J. Soil Sci.*, **60**(1), 94–109.

Ritzwoller, M.H. & Feng, L., 2019. Overview of pre- and post-processing of ambient-noise correlations, in Seismic Ambient Noise, eds, Nakata, N., Gualtieri, L. & Fichtner, A, Cambridge University Press, Cambridge, UK, pp. 144–187, doi: 10.1017/9781108264808.

Norris, A.N., 2024. Finite-amplitude waves in solids, in *Nonlinear Acoustics*, pp. 259–273, Springer-Verlag, Berlin.

Obermann, A., Planès, T., Larose, E., Sens-Schönfelder, C. & Campillo, M., 2013. Depth sensitivity of seismic coda waves to velocity perturbations in an elastic heterogeneous medium, *J. Geophys. Int.*, **194**(1), 372–382.

Obermann, A., Planès, T., Hadzioannou, C. & Campillo, M., 2016. Lapse-time-dependent coda-wave depth sensitivity to local velocity perturbations in 3-D heterogeneous elastic media, *J. Geophys. Int.*, **207**(1), 59–66.

Obermann, A., Planès, T., Larose, E. & Campillo, M., 2019. 4-D imaging of subsurface changes with coda waves: numerical studies of 3-D combined sensitivity kernels and applications to the Mw 7.9, 2008 Wenchuan earthquake, *Pure appl. Geophys.*, **176**(3), 1243–1254.

Olivier, G. & Brenguier, F., 2016. Interpreting seismic velocity changes observed with ambient seismic noise correlations, *Interpretation*, **4**(3), SJ77-SJ85.

Pan, Y., Shen, W.-B., Shum, C. & Chen, R., 2018. Spatially varying surface seasonal oscillations and 3-D crustal deformation of the Tibetan plateau derived from GPS and grace data, *Earth planet. Sci. Lett.*, **502**, 12–22.

Payan, C., Garnier, V., Moysan, J. & Johnson, P., 2009. Determination of third order elastic constants in a complex solid applying coda wave interferometry, *Appl. Phys. Lett.*, **94**(1), 011904, doi:10.1063/1.3064129.

Pellet, C. & Hauck, C., 2017. Monitoring soil moisture from middle to high elevation in switzerland: set-up and first results from the somomount network, *Hydrol. Earth Syst. Sci.*, **21**(6), 3199–3220.

Poli, P., Marguin, V., Wang, Q., d'Agostino, N. & Johnson, P., 2020. Seasonal and coseismic velocity variation in the region of l'aquila from single station measurements and implications for crustal rheology, *J. Geophys. Res. Solid Earth*, **125**(7), e2019JB019316, doi:10.1029/2019JB019316.

Power, G., Brown, R. & Imhof, J., 1999. Groundwater and fish—insights from northern north America, *Hydrol. Process.*, **13**(3), 401–422.

Razandi, Y., Pourghasemi, H.R., Neisani, N.S. & Rahmati, O., 2015. Application of analytical hierarchy process, frequency ratio, and certainty factor models for groundwater potential mapping using gis, *Earth Sci. Inform.*, **8**, 867–883.

Richter, T., Sens-Schönfelder, C., Kind, R. & Asch, G., 2014. Comprehensive observation and modeling of earthquake and temperature-related seismic velocity changes in northern Chile with passive image interferometry, *J. Geophys. Res. Solid Earth*, **119**(6), 4747–4765.

Riel, B., Simons, M., Ponti, D., Agram, P. & Jolivet, R., 2018. Quantifying ground deformation in the Los Angeles and Santa Ana coastal basins due to groundwater withdrawal, *Water Resour. Res.*, **54**(5), 3557–3582.

Ruland, W., Cherry, J. & Feenstra, S., 1991. The depth of fractures and active ground-water flow in a clayey till plain in southwestern Ontario, *Groundwater*, **29**(3), 405–417.

Sarkar, D., Bakulin, A. & Kranz, R.L., 2003. Anisotropic inversion of seismic data for stressed media: Theory and a physical modeling study on berea sandstone, *Geophysics*, **68**(2), 690–704.

Save, H., 2020. CSR GRACE and GRACE-FO RL06 Mascon Solutions v02, doi:10.15781/cqg9-nh24.

Save, H., Bettadpur, S. & Tapley, B.D., 2016. High-resolution csr grace rl05 mascons, *J. Geophys. Res. Solid Earth*, **121**(10), 7547–7569.

Sens-Schönfelder, C., 2008. Synchronizing seismic networks with ambient noise, *J. Geophys. Int.*, **174**(3), 966–970.

Sens-Schönfelder, C. & Larose, E., 2008. Temporal changes in the lunar soil from correlation of diffuse vibrations, *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.*, **78**(4), 045601, doi:10.1103/PhysRevE.78.045601.

Sens-Schönfelder, C. & Wegler, U., 2006. Passive image interferometry and seasonal variations of seismic velocities at merapi volcano, Indonesia, *Geophys. Res. Lett.*, **33**(21), doi:10.1029/2006GL027797.

Shapiro, N.M., Ritzwoller, M.H. & Bensen, G.D., 2006. Source location of the 26 sec microseism from cross-correlations of ambient seismic noise, *Geophys. Res. Lett.*, **33**(18), doi:10.1029/2006GL027010.
Shen, Z., Yang, Y., Fu, X., Adams, K.H., Biondi, E. & Zhan, Z., 2024. Fiber-optic seismic sensing of vadose zone soil moisture dynamics, *Nat. Commun.*, **15**(1), 6432, doi:10.1038/s41467-024-50690-6.

Silver, P.G., Daley, T.M., Niu, F. & Majer, E.L., 2007. Active source monitoring of cross-well seismic travel time for stress-induced changes, *Bull. seism. Soc. Am.*, **97**(1B), 281–293.

Somers, L.D. & McKenzie, J.M., 2020. A review of groundwater in high mountain environments, *Wiley Interdiscip. Rev.: Water*, **7**(6), e1475, doi: 10.1002/wat2.1475.

Song, C., Huang, B., Richards, K., Ke, L. & Hien Phan, V., 2014. Accelerated lake expansion on the Tibetan plateau in the 2000s: induced by glacial melting or other processes?, *Water Resour. Res.*, **50**(4), 3170–3186.

Stehly, L., Campillo, M. & Shapiro, N., 2007. Traveltime measurements from noise correlation: stability and detection of instrumental time-shifts, *J. Geophys. Int.*, **17**(1), 223–230.

Sun, J. et al., 2024. Soil water components control plant water uptake along a subalpine elevation gradient on the eastern qinghai-tibet plateau, *Agric. For. Meteorol.*, **345**, 109827, doi:10.1016/j.agrformet.2023.109827.

Tan, H., Zhang, Y., Rao, W., Guo, H., Ta, W., Lu, S. & Cong, P., 2021. Rapid groundwater circulation inferred from temporal water dynamics and isotopes in an arid system, *Hydrol. Process.*, **35**(6), e14225, doi:10.1002/hyp.14225.

Taylor, R.G. et al., 2013. Ground water and climate change, *Nat. Clim. Change*, **3**(4), 322–329.

Tsai, VC., 2011. A model for seasonal changes in GPS positions and seismic wave speeds due to thermoelastic and hydrologic variations, *J. Geophys. Res. Solid Earth*, **116**(B4), doi:10.1029/2010JB008156.

Viviroli, D., Dürr, H.H., Messerli, B., Meybeck, M. & Weingartner, R., 2007. Mountains of the world, water towers for humanity: Physiology, mapping, and global significance, *Water Resour. Res.*, **43**(7), doi:10.1029/2006WR005653.

Wang, G., Wang, Y., Li, Y. & Cheng, H., 2007. Influences of alpine ecosystem responses to climatic change on soil properties on the Qinghai-Tibet plateau, China, *Catena*, **70**(3), 506–514.

Wang, Q.-Y., Brenguier, F., Campillo, M., Lecointre, A., Takeda, T. & Aoki, Y., 2017. Seasonal crustal seismic velocity changes throughout Japan, *J. Geophys. Res. Solid Earth*, **122**(10), 7987–8002.

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, *J. Geophys. Res. Solid Earth*, **124**(8), 8924–8941.

Wang, Q.-Y., Cui, X., Frank, W.B., Lu, Y., Hirose, T. & Obara, K., 2024. Untangling the environmental and tectonic drivers of the Noto earthquake swarm in Japan, *Sci. Adv.*, **10**(19), eado1469, doi:10.1126/sciadv.ad01469.

Wang, S., Liu, F., Zhou, Q., Chen, Q. & Liu, F., 2021. Simulation and estimation of future ecological risk on the Qinghai-Tibet plateau, *Sci. Rep.*, **11**(1), 17603, doi:10.1038/s41598-021-96958-5.

Wang, T., Yang, D., Yang, Y., Piao, S., Li, X., Cheng, G. & Fu, B., 2020. Permafrost thawing puts the frozen carbon at risk over the Tibetan plateau, *Sci. Adv.*, **6**(19), eaaz3513, doi:10.1126/sciadv.eaz3513.

Wang, T., Yang, D., Yang, Y., Zheng, G., Jin, H., Li, X., Yao, T. & Cheng, G., 2023. Unsustainable water supply from thawing permafrost on the Tibetan plateau in a changing climate, *Sci. Bull.*, **68**(11), 1105–1108.

Wang, X., Yang, M., Liang, X., Pang, G., Wan, G., Chen, X. & Luo, X., 2014. The dramatic climate warming in the Qaidam basin, northeastern Tibetan plateau, during 1961-2010, *Int. J. Climatol.*, **34**(5), 1524–1537.

Wang, Y. et al., 2008. *Investigation and Assessment of the Groundwater Resources and Their Environmental Issues in the Qaidam Basin (in Chinese)*, Beijing Geology Press.

Watkins, M.M., Wiese, D.N., Yuan, D.-N., Boening, C. & Landerer, F.W., 2015. Improved methods for observing earth's time variable mass distribution with grace using spherical cap mascons, *J. Geophys. Res. Solid Earth*, **120**(4), 2648–2671.

Wegler, U., Lühr, B.-G., Snieder, R. & Ratdomopurbo, A., 2006. Increase of shear wave velocity before the 1998 eruption of Merapi volcano (Indonesia), *Geophys. Res. Lett.*, **33**(9), doi:10.1029/2006GL025928.

Weinert, S., Bär, K. & Sass, I., 2020. Thermophysical rock properties of the crystalline Gonghe basin complex (northeastern Qinghai-Tibet-plateau, China) basement rocks, *Environ. Earth Sci.*, **79**(3), 77, doi:10.1007/s12665-020-8808-9.

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

Wiese, D.N., Yuan, D.-N., Landerer, F.W. & Watkins, M.M., 2023. JPL GRACE and GRACE-FO Mascon ocean, ice, and hydrology equivalent water height coastal resolution improvement (CRI) filtered release 06.1 version 03, NASA Physical Oceanography Distributed Active Archive Center, doi:10.5067/TEMSC-3JC63.

Winkler, K.W. & Liu, X., 1996. Measurements of third-order elastic constants in rocks, *J. acoust. Soc. Am.*, **100**(3), 1392–1398.

Winkler, K.W. & McGowan, L., 2004. Nonlinear acoustoelastic constants of dry and saturated rocks, *J. Geophys. Res. Solid Earth*, **109**(B10), doi:10.1029/2004JB003262.

Wu, S.-M., Lin, F.-C., Farrell, J., Shiro, B., Karlstrom, L., Okubo, P. & Koper, K., 2020. Spatiotemporal seismic structure variations associated with the 2018 Kīlauea eruption based on temporary dense geophone arrays, *Geophys. Res. Lett.*, **47**(9), e2019GL086668, doi:10.1029/2019GL086668.

Wu, Y., & 2010. Mechanism of action of cracks water on rock landslide in rainfall, *J. Central South Univ. Technol.*, **17**(6), 1383–1388.

Xiang, L., Wang, H., Steffen, H., Wu, P., Jia, L., Jiang, L. & Shen, Q., 2016. Groundwater storage changes in the tibetan plateau and adjacent areas revealed from grace satellite gravity data, *Earth planet. Sci. Lett.*, **449**, 228–239.

Xiao, X., Cheng, S., Wu, J., Wang, W., Sun, L., Wang, X. & Wen, L., 2021. Shallow seismic structure beneath the continental China revealed by P-wave polarization, Rayleigh wave ellipticity and receiver function, *J. Geophys. Int.*, **225**(2), 998–1019.

Xiao, Y., 2018. Groundwater circulation patterns and its change trend in southern Qaidam basin, northwest China (in Chinese), *Ph.D. thesis*, China University of Geosciences, Beijing.

Xiao, Y., Shao, J., Cui, Y., Zhang, G. & Zhang, Q., 2017. Groundwater circulation and hydrogeochemical evolution in Nomhon of Qaidam Basin, northwest China, *J. Earth Syst. Sci.*, **126**(26), doi:10.1007/s12040-017-0800-8.

Xie, F., Larose, E., Wang, Q. & Zhang, Y., 2023. In-situ monitoring of rock slope destabilization with ambient seismic noise interferometry in southwest china, *Eng. Geol.*, **312**, 106922, doi:10.1016/j.enggeo.2022.106922.

Xie, J., Xu, Y.-P., Guo, Y., Wang, Y. & Chen, H., 2022. Understanding the impact of climatic variability on terrestrial water storage in the Qinghai-Tibet plateau of China, *Hydrol. Sci. J.*, **67**(6), 963–978.

Xu, Z.J. & Song, X., 2009. Temporal changes of surface wave velocity associated with major sumatra earthquakes from ambient noise correlation, *Proc. Natl. Acad. Sci. USA*, **106**(34), 14207–14212.

Xue, X., Guo, J., Han, B., Sun, Q. & Liu, L., 2009. The effect of climate warming and permafrost thaw on desertification in the Qinghai–Tibetan plateau, *Geomorphology*, **108**(3-4), 182–190.

Yang, M., Wang, X., Pang, G., Wan, G. & Liu, Z., 2019. The Tibetan Plateau cryosphere: observations and model simulations for current status and recent changes, *Earth Sci. Rev.*, **190**(December 2018), 353–369.

Yao, H., van der Hilst, R.D. & de Hoop, M.V., 2006. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis - I. Phase velocity maps, *J. Geophys. Int.*, **166**(2), 732–744.

Yao, H., Gouédard, P., Collins, J.A., McGuire, J.J. & van der Hilst, R.D., 2011. Structure of young east pacific rise lithosphere from ambient noise correlation analysis of fundamental- and higher-mode Scholte-Rayleigh waves, *C. R. Geosci.*, **343**(8–9), 571–583.

Yao, T. et al., 2022. The imbalance of the Asian water tower, *Nat. Rev. Earth Environ.*, **3**(10), 618–632.
Ye, S., 2022. Meteorological data set of six stations in Golmud river basin and its vicinity (1956-2018), National Cryosphere Desert Data Center, doi:10.12072/ncdc.Hydro.db2446.2022.

Yin, T., Yin, J., Wu, Y., Yang, Z., Liu, X. & Zhuang, D., 2022. Water saturation effects on the mechanical characteristics and fracture evolution of sandstone containing pre-existing flaws, *Theor. Appl. Fract. Mech.*, **122**, 103605, doi:10.1016/j.tafmec.2022.103605.

Yu, C., Wang, Q.-Y., Ma, J. & Yao, H., 2025. *Nuomuhong daily auto-correlation functions and average relatively seismic velocity change results*, Zenodo, doi:10.5281/zenodo.14690410.

Yuan, C., Bryan, J. & Denolle, M., 2021. Numerical comparison of time-, frequency- and wavelet-domain methods for coda wave interferometry, *J. Geophys. Int.*, **226**(2), 828–846.

Zhan, Z., Tsai, V.C. & Clayton, R.W., 2013. Spurious velocity changes caused by temporal variations in ambient noise frequency content, *J. Geophys. Int.*, **194**(3), 1574–1581.

Zhang, J., 2013. Groundwater resource evaluation of Xiangride-Nuomuhong piedmont plain (in Chinese), *Master's thesis*, China University of Geosciences, Beijing.

Zhang, Y. et al., 2018. 3-D crustal shear-wave velocity structure of the Taiwan strait and Fujian, SE China, revealed by ambient noise tomography, *J. Geophys. Res. Solid Earth*, **123**(9), 8016–8031.

Zhang, Y., Tan, H., Cong, P., Shi, D., Rao, W. & Zhang, X., 2023. Isotopic variations in surface waters and groundwaters of an extremely arid basin and their responses to climate change, *Hydrol. Earth Syst. Sci.*, **27**(21), 4019–4038.

Zhao, W., Lin., Y. Z., Zhou., P.P., Wang., G.C., Dang., X.Y. & Gu., X.F., 2021. Characteristics of groundwater in northeast Qinghai-Tibet plateau and its response to climate change and human activities: a case study of Delingha, Qaidam basin, *China Geol.*, **4**(3), 377–388.

Zhong, Z. et al., 2023. Disentangling the effects of vapor pressure deficit on northern terrestrial vegetation productivity, *Sci. Adv.*, **9**(32), eadf3166, doi:10.1126/sciadv.adf3166.

Zou, D. et al., 2024. Estimation of permafrost ground ice to 10 m depth on the Qinghai-Tibet plateau, *Permafr. Periglac. Process.*, **35**(3), 423–434.

Zou, Y. et al., 2022. Solid water melt dominates the increase of total groundwater storage in the Tibetan plateau, *Geophys. Res. Lett.*, **49**(18), e2022GL100092, doi:10.1029/2022GL100092.