Study on Response Process and Time Delay Effect of Groundwater Dynamic in Northeastern Margin of Tibetan Plateau

Abstract

Under the background of drastic global climate change, the evolution law of groundwater resources in the northeastern margin of the Tibetan Plateau presents new characteristics, and the groundwater is gradually becoming more complicated, diversified and disordered. In this study, cross-correlation analysis, wavelet analysis and cross-wavelet transform were used to explore the response mechanism and time delay effect of groundwater, exploitation amount, rainfall and surface runoff in the northeastern margin of the Tibetan Plateau. The results show that the groundwater depth increased with the increase in the exploitation amount and decreased with the increase in the rainfall-to-exploitation ration and the surface runoff-to-exploitation ratio from 1980 to 2020. On the annual scale, groundwater, rainfall and surface runoff had a strong cohesiveness period of 12a. On the monthly scale, groundwater lagged behind rainfall and surface runoff for 3 months and 2 months, respectively. The above conclusions provide a scientific theoretical basis for deepening the characteristics of groundwater endowment, the evolution law of water cycle elements and the theory of the hydrological cycle in the northeastern margin of the Tibetan Plateau.

1. Introduction

Groundwater has multiple attribute functions such as resource, environment and ecology, and it is an indispensable geological resource in production, life and ecology all over the world [1,2]. The Tibetan Plateau, known as the “headwaters of the Three Rivers” and the “Great Water Tower of Asia”, has abundant reserves of groundwater resources. The average amount of groundwater resources over the years is about 1500 × 10⁸ m³, accounting for 17% of the total groundwater resources in China [3,4]. Abundant groundwater not only plays an important role in water supply in the sustainable economic development of the Tibetan Plateau, but also plays an important strategic role in ensuring water resource security in China and southeast Asia. Known as the “third pole” of the Earth, the Tibetan Plateau is the initiating area of climate change in Asia and a sensitive area of ecological environment in response to nature and human beings [5,6]. Since the 1970s, the climatic characteristics of the Tibetan Plateau have undergone a significant shift from warm and dry to warm and wet, showing a trend of a continuous rise in temperature and a significant increase in rainfall, resulting in significant changes in the spatial and temporal distribution pattern of groundwater on the Plateau [7,8,9]. Xining is located in the northeast margin of the Tibetan Plateau and the Huangshui River valley basin. The formation in this area is Neogene red mudstone rich in soluble salt, and coupled with low rainfall, the surface water quality is poor [10]. The groundwater in the Huangshui River and its tributary river valley has become an important source of water for the economic development of local towns and people’s livelihood with its excellent water quality, abundant water volume and convenient and economical exploitation conditions [11]. In recent years, with the in-depth implementation of Western development in the new era and the “One Belt, One Road” strategy, the pace of local urban development has been accelerating. Under the background of global climate change, the groundwater in the northeastern margin of the Tibetan Plateau presents new characteristics, and is gradually becoming more complicated, diversified and disordered [12].

Research on groundwater response mechanisms has previously been carried out. The United Nations Educational, Scientific and Cultural Organization (UNESCO) established the multi-phase International Hydrological Program (IHP) in 1950. After 1980, a series of large-scale scientific research projects such as the Global Energy and Water Cycle Experiment (GEWEX) and the Regional Groundwater System Research Program (RGSPR) were launched, which promoted research on the evolution law and driving mechanism of groundwater resources at different scales. In 2005, the International Association of Hydrological Sciences (IAHS) emphasized the need to strengthen the research on the mechanism of changes in the hydrological cycle under the influence of natural changes and human activities. Since then, scholars from various countries have attached great importance to the stability and sustainability of regional groundwater systems, and have successively carried out research on the dynamic evolution of groundwater and related driving factors, such as the global area [13,14], the United States [15,16,17], South Korea [18], South Africa [19,20], and China [21,22]. There are many research methods on groundwater response and time delay, mainly focusing on multivariate statistical methods [23], numerical simulation [24,25], and wavelet analysis [26,27]. Relevant scholars have developed and innovated many technical methods and means based on time and place, reasonably revealing the response process and time delay effect of groundwater under the influence of natural changes and human activities. Cross-wavelet transform (CWT) can accurately identify the correlation and lag time between two time series with the advantage of localization features in the time domain and frequency domain. In recent years, more and more scholars have applied cross-wavelet transform to analyze the correlation and lag time between two groups of signals, such as precipitation and spring discharge [28], rainfall and dipole [29] and rainfall and runoff [30].

Due to the topographical and geomorphological differences of the Tibetan Plateau, the concealability of groundwater and the difficulty in obtaining data, early groundwater monitoring data are missing, and the research on groundwater resources was carried out late. Research on groundwater response mechanisms has been difficult, with few research results. The cross-wavelet transform method has yet to be applied and verified in the groundwater of the Qinghai–Tibet Plateau. Therefore, to promote ecological protection and devise a high-quality development strategy in the Yellow River Basin, it is of great theoretical and practical significance to study the response mechanism and time delay effect of groundwater in terms of natural driving factors in the northeastern margin of the Tibetan Plateau for realizing sustainable groundwater utilization. Based on this, the Xining region was taken as the research area in this work, and the response process and time delay effect of groundwater levels from 1980 to 2020 were determined using cross-correlation analysis, wavelet analysis and cross-wavelet transform. This paper aims to (1) explore the response process of the groundwater and driving factors; (2) clarify the periodic evolution law of the groundwater and driving factors by using continuous wavelet transform; (3) reveal the time delay effect of groundwater on driving factors by using cross-wavelet transform. The results of this work can provide a scientific theoretical basis for the safety of groundwater resources in the Tibetan Plateau and accurate decision making by the government, and are expected to provide examples and new ideas for groundwater research in other regions.

2. Materials and Methods

2.1. Study Area and Data Source

Xining is the city with the largest area, the most concentrated population and the most developed economy on the Tibetan Plateau (in Figure 1a). Its geographical coordinates are 100°58′48″–102°01′26″ east longitude and 36°24′40″–37°03′33″ north latitude. The research area is about 700 km², with jurisdiction over five districts and two counties. The study area is at the junction of three natural geographical areas, namely, the inland arid area in northwest China, the humid monsoon area in east China, and the alpine region of the Tibetan Plateau. It is also a sensitive ecological environment area responding to nature and human beings [31]. With typical semi-arid plateau continental climate characteristics (mainly cold and dry), the average annual temperature is 5.9 °C, and the average annual evaporation is 1708.4 mm. The average annual precipitation is 393.6 mm, and the annual distribution of precipitation is uneven, with 70% of the rainfall mainly concentrated from July to September [32]. Xining city is located in the middle and upper reaches of the Huangshui River basin, which distributes 56 rivers including the Huangshui River (Figure 1b). The average annual discharge of the Huangshui River is 32.2 m³/s, and the annual discharge in dry and rainy years is 12.5 m³/s and 63.3 m³/s, respectively. The distribution of the river flow is uneven across the year, and the flood season mainly occurs from July to September, accounting for about 70% of the annual runoff [33]. The Quaternary system is widely developed in the study area, mainly including the Middle Pleistocene (Q₂), Upper Pleistocene (Q₃) and Holocene (Q₄), and the underlying stratum is horizontal Tertiary mudstone with sandstone laterite (in Figure 1c). The Quaternary loose rock porewater is widely distributed in the plain area of the Huangshui Valley and the high terraces of the hilly areas on both sides, which is the main type of groundwater in the study area. The groundwater is rich in quantity in the middle and upper reaches of the Huangshui River and its larger tributaries, such as Beichuan River, Xinachuan River and Nanchuan River. Most of these waters are extremely fresh, with a hydrochemical type of HCO₃-Ca and TDS of 150–500 mg/L, which are the main water sources for exploitation and utilization in Xining [34].

The groundwater dynamic monitoring work in Xining has already begun, with the monitoring object being the groundwater depth information of Quaternary loose rock porewater. The groundwater monitoring well is located in the urban area of Xining, as shown in Figure 1b. The monitoring well continuously measures the depth of phreatic groundwater by automatically collecting water 24 times per day. Taking Xining groundwater dynamics as the research object, this paper analyzes the response mechanism and time delay effect of the groundwater dynamic using the average groundwater level data from 1980 to 2020. The groundwater monitoring data are mainly provided by the Qinghai Provincial Geological Environment Monitoring Station, the Qinghai Provincial Ecological Environment Monitoring Center and the Hydrogeological Environmental Geological Survey Center of China Geological Survey. The data from meteorological and hydrology stations are taken from the Qinghai Hydrology and Water Resources Monitoring Center.

2.2. Research Methods

2.2.1. Cross-Correlation Function

Cross-correlation function (CCF) refers to the correlation degree of two groups of sequence variables at any two moments. The correlation degree of the delay time of different sequences can be obtained by calculating the time difference dislocation movement, which can describe the relationship between the internal correlation degree and the relative lag between sequences [35]. The time corresponding to the first maximum correlation coefficient (peak) represents the time that the input signal lags behind the output signal. Given the homogeneous time series signals X_t and Y_t (t = 0, ±1…, ±d, d is the delay time), the cross-correlation coefficient r of X_t at time t and Y_t at time t + d is defined as [36,37]:

$$cc{f}_d=f\left(X_t,Y_{t+d}\right)=r_{X_t{Y}_{t+d}}=\frac{{\displaystyle \sum \left(X_t-\overline{X_t}\right)\left(Y_{t+d}-\overline{Y_{t+d}}\right)}}{\sqrt{{\displaystyle \sum {\left(X_t-\overline{X_t}\right)}^2{\displaystyle \sum {\left(Y_{t+d}-\overline{Y_{t+d}}\right)}^2}}}}$$

(1)

2.2.2. Wavelet Analysis

The wavelet analysis method has been developed as the main analysis method of signal analysis and time scale recognition, and is widely used in image processing [38], signal diagnosis [39], hydrometeorological sequence recognition [40,41], and other fields. Time domain and frequency domain are the basic forms of time series. Among them, time-domain analysis and frequency-domain analysis have accurate time and frequency positioning capabilities, respectively. Wavelet analysis is transformed from the Fourier formula, which can simultaneously reflect the time-domain and frequency-domain characteristics of non-stationary sequences, analyze the internal fine structure and extract the hidden regularity. Its time-frequency localization and multi-level resolution function can reveal the details of local changes (high frequency) and the overall evolution trend (low frequency) of the sequence. Groundwater data are a continuous and non-stationary time series with seasonal and periodic characteristics of multiple time scales. Therefore, the complex Morlet continuous wavelet transform is selected as the wavelet function in this paper to avoid the false oscillation of real wavelet transform coefficients and make the analysis results more accurate and reliable. The complex continuous wavelet coefficient function is

$$\begin{array}{c}Wf\left(a,b\right)={\left|a\right|}^{-1/2}\Delta t{\displaystyle \sum _{k-1}^{N}f\left(k\mathsf{\Delta }t\right)}\overline{\psi }\left(\frac{k\Delta t-b}{a}\right)\\ {\displaystyle {\int }_{-\infty }^{+\infty }\psi \left(t\right)dt=0}\end{array}$$

(2)

where f(t) is the original sequence; ψ(t) is the basis wavelet function; a is the scale factor, representing the time period length of the sequence; b is the time factor, reflecting the translation of the sequence in time; t is the time interval; W_f(a,b) is the Morlet wavelet transform coefficient, which represents the correlation strength between the original sequence and the wavelet function at the scale of a and the displacement of b.

The coefficients obtained by the complex Morlet wavelet transform are in complex form. The isoline of the real part of the coefficient can reflect the energy intensity information in different phases at different time scales. When the scale factor a is the same, the variation process of the wavelet coefficient with time t represents the evolution law of the time series at this scale. The positive and negative change in the real part of the coefficient represents the alternating change rule of the time series. The positive value of the coefficient corresponds to more periods of the series, the zero value corresponds to more and less transitional periods, and the negative value corresponds to fewer periods [42]. The magnitude of the modulus square of the wavelet system reflects the oscillation strength of the signal at different timescales and the energy distribution at a specific scale. The integral of this value along with the scale factor a on the time factor b is called the wavelet square difference function, and its expression is

$$Var(a)={{\displaystyle {\int }_{-\infty }^{\infty }\left|W_f(a,b)\right|}}^{2}db$$

(3)

The wavelet square difference can describe the oscillation intensity of the time series in different principal periods, and can reflect the distribution law of the principal period and wavelet coefficient with the change in scale. The larger the variance is, the more prominent the periodicity on this scale. The maximum value indicates the strongest periodic oscillation here, and the corresponding scale a is the first main period.

2.2.3. Cross-Wavelet Transform

Cross-wavelet transform (CWT) and wavelet transform coherence (WTC) can not only diagnose the similarity between signals and obtain the high and low value regions of their mutual oscillations in the time–frequency domain, but can also analyze the delay time in different frequency scales, thus revealing the same change period regions in the different time scales of the two groups of sequences [43]. Among them, cross-wavelet transform can analyze the phase relationship of the high-energy spectrum region, while the wavelet coherence can analyze the phase relationship in the low-energy spectrum region. Therefore, it is possible to intuitively judge the change details, phase differences, local features and time lag effects of the correlation between the two groups of sequences in the time–frequency domain with different time scales [44]. The calculation principle is to combine complex Morlet wavelets with cross-spectral analysis to carry out cross-wavelet transform and wavelet coherence. The cross-wavelet spectrum of the two time series is defined as:

$$W_n^{xy}\left(s\right)=W_n^x\left(s\right)W_n^{y\ast }\left(s\right)$$

(4)

where $W_n^{y\ast }\left(s\right)$ is the W* complex conjugate; $W_n^{xy}\left(s\right)$ is the complex number; and $\left|W_n^{xy}\left(s\right)\right|$ is the cross-wavelet power spectrum. The complex angle of W^xy in the complex plane is the phase relation between the sequence x and y in the time–frequency domain.

By comparing the crossover wavelet power spectrum with the red noise standard spectrum, the periodic correlation between the two sequences in different time scales can be obtained [45]. The red noise curve is simulated by the first-order autoregressive process (AR1), and its standard spectral distribution is as follows:

$$D\left(\frac{\left|W_n^x\left(s\right)W_n^{y\ast }\left(s\right)\right|}{{\sigma }_x{\sigma }_y}<p\right)=\frac{z_v\left(p\right)}{v}\sqrt{P_k^x{P}_k^y}$$

(5)

Wavelet coherence is the local correlation coefficient of two groups of time series in the time–frequency domain, which is used to measure the degree of correlation between the series. The Monte Carlo method is used to test the significance of the coefficient [46,47]. It is defined as:

$$R_n^2\left(S\right)=\frac{{\left|S\left(s^{-1}{w}_n^{xy}\left(s\right)\right)\right|}^2}{S\left(s^{-1}{\left|W_n^x\left(s\right)\right|}^2\right)\cdot S\left(s^{-1}{\left|W_n^y\left(s\right)\right|}^2\right)}$$

(6)

where $R_n^2\left(S\right)$ is the local correlation coefficient and S is the smoothing operator.

3. Results and Discussion

3.1. Identification of the Main Driving Factors on Groundwater

Natural factors (rainfall, evaporation, temperature, surface runoff) and human factors (population number, gross regional domestic product, construction land area, agricultural acreage, groundwater exploitation quantity) were selected as nine driving factors. The grey correlation method was used to analyze the correlation degree between the driving factors and groundwater depth during the period 1980–2020, as shown in Figure 2 and

Table 1. Summary of the correlation degree of groundwater depth and driving factors.

Driving Factor		Groundwater Depth	Order
Natural factors	Rainfall	0.816	2
	Evaporation	0.623	5
	Temperature	0.373	9
	Surface runoff	0.736	3
Human factors	Population number	0.709	4
	Gross regional domestic product	0.429	8
	Construction land area	0.539	6
	Agricultural acreage	0.506	7
	Exploitation quantity	0.862	1

. As can be seen from Figure 2 and Table 1, the overall correlation was ranked as exploitation quantity (0.862) > rainfall (0.816) > surface runoff (0.756) > population number (0.709) > evaporation (0.623) > construction land area (0.539) > agricultural acreage (0.506) > gross regional domestic product (GDP) (0.429) > temperature (0.373). Among them, the correlation degree of exploitation quantity, rainfall and surface runoff were more than 0.75, indicating that these three factors had a strong correlation with the groundwater depth, further reflecting that the dynamic of groundwater buried depth is influenced by both human and natural factors. The correlation degree of population number, evaporation and groundwater depth is lower than 0.7, and the correlation degree is relatively low. Groundwater is the main source of water supply in Xining city. The increase in the population will inevitably lead to the increase in water demand, so the impact of population on groundwater is mainly reflected in the amount of groundwater extraction. Xining city is a typical semi-arid plateau continental climate, and its evaporation is more intense. When the depth of the groundwater table is less than 5 m, the groundwater is greatly affected by evaporation. However, the buried depth of the groundwater level in this area is about 11 m (>5 m), so evaporation is not considered the main influencing factor of groundwater dynamics. The correlation degree of construction land area, agricultural acreage, GDP and temperature is less than 0.6, which can ignore its influence on the interannual change in groundwater depth. To sum up, the exploitation quantity, rainfall and surface runoff are the main driving factors for the dynamic change in the groundwater depth in Xining city, in which the extraction amount is the main factor and the rainfall and surface runoff are the secondary factors.

3.2. Response Process of Groundwater to Main Driving Factors

3.2.1. Response Process of Groundwater to Exploitation

Human exploitation activity is a direct factor in groundwater dynamic change in the study area. Figure 3 shows the relation curve between groundwater depth and exploitation amount. As can be seen from Figure 3a, the variation trend of groundwater depth and exploitation volume in Xining from 1980 to 2020 is relatively consistent, and the groundwater depth increases with the increase in the extraction volume, showing an overall trend of enlarging–decreasing–enlarging–decreasing. Its evolution process is divided into four stages: (1) The large-scale exploitation of groundwater resources in the early stage of the construction of Xining increased from 0.4 × 10⁸ m³ in 1980 to 2.2 × 10⁸ m³ in 2001, resulting in an increase in groundwater depth from 10.7 m in 1980 to 16.2 m in 2001. (2) After 2001, the groundwater exploitation gradually decreased to 1.8 × 10⁸ m³ in 2005, and the groundwater depth gradually decreased to 13.9 m. (3) After 2006, the recovery increased to 2.3 × 10⁸ m³/a in 2012, and the groundwater depth also increased to 14.6 m. (4) Since then, the mining yield has been gradually decreasing, and it remains at about 2.0 × 10⁸ m³/a after 2016, and the groundwater depth decreased to 13.6 m in 2020. According to Figure 3b, there is a positive correlation between groundwater depth and exploitation amount from 1980 to 2020. When the exploitation volume is small, the slope of the power function curve is steep. With the increase in exploitation, the groundwater depth showed a power function increasing trend, but the increasing speed slowed down. With the gradual increase in the extraction amount, the dot started to deviate from the curve, indicating that the dynamic evolution of groundwater tends to be complicated in the process of urbanization.

3.2.2. Response Process of Groundwater to Rainfall

(1)

Response process of groundwater evolution to annual rainfall

Rainfall is one of the important sources of groundwater recharge in Xining. Based on the monitoring data from 1980 to 2020, the variation rules of groundwater depth and rainfall are drawn, as shown in Figure 4. As can be seen from Figure 4a, the rainfall shows a cyclical increasing trend, increasing from 274.8 mm in 1980 to 516.5 mm in 2020. Under natural conditions, the groundwater depth varies periodically with rainfall. However, the burial depth increases first and then decreases, from 10.7 m in 1980 to 16.2 m in 2001 and then to 13.6 m in 2020. This is mainly because under the influence of factors such as extraction amount, groundwater depth no longer presents a seasonal change with rainfall. According to the columnar distribution of rainfall anomalies, the rainfall in rainy years was abundant, groundwater recharge increased, the groundwater depth was in a decreasing stage, and increased rainfall had a positive effect on groundwater recovery, such as in 2003, 2007 and 2019. In dry years, the rainfall was relatively poor, groundwater recharge decreased, the groundwater depth was at an increasing stage, and the decrease in rainfall had a negative effect on groundwater recovery, such as in 1982 and 1991. Therefore, rainfall is still an important factor affecting the dynamic change in the groundwater in Xining. With the continuous progress of urbanization, the groundwater in Xining is greatly affected by human exploitation activities, and the response degree of groundwater depth to rainfall is gradually weakened. First, the increase in the impervious interface of the underlying surface reduces the rainfall infiltration coefficient and, thus, the effective recharge of the groundwater. Second, due to the increase in the groundwater depth caused by human exploitation activities, rainfall needs to make up the water deficit in the very thick aeration zone before it can recharge the groundwater.

In order to explore the influence of rainfall on groundwater dynamic change characteristics under mining conditions, the correlation between the groundwater depth and the ratio of rainfall to exploitation quantity was further analyzed. According to Figure 4b, there is a negative correlation between groundwater depth and the rainfall-to-extraction ratio from 1980 to 2020. When the extraction amount is constant, the groundwater depth decreases in a power function with the increase in rainfall, and the decreasing speed slows down. When the rainfall-to-extraction quantity ratio is small, the slope of the power function curve is steeper, indicating that the influence of rainfall on groundwater dynamics in dry years is greater than that in rainy years. When the rainfall-to-extraction ratio increases, the groundwater depth dispersion becomes dispersed, indicating that the groundwater evolution mechanism becomes more complicated in rainy years.

(2)

Response process of groundwater evolution to monthly rainfall

The variation of rainfall in long time series is a superimposed process of mutation, tendency and periodicity. From the perspective of climate factors, the range of groundwater dynamic change depends on the superimposed effect of multiple rainfall events, and there is generally a linear quantitative relationship between monthly groundwater depth and rainfall. However, under objective conditions such as human factors, the change in groundwater depth often lags behind the rainfall. According to the variation curve of monthly rainfall in Figure 5, the rainfall presents a significant periodic variation feature. The rainy season is from July to September each year with relatively abundant rainfall, accounting for 56.2–77.4% of the total annual rainfall. The highest rainfall of 175.6 mm occurred in August 2007. The dry season is from October to June of the next year with small amounts of rainfall. The minimum rainfall is 0, which tends to occur in January and December each year. According to the changes in the groundwater depth and rainfall of different hydrogeological units shown in Figure 5, the dynamic change characteristics of groundwater depth obviously lag behind the rainfall. For example, the peak of the groundwater depth in washland occurred in August 1985, September 1989, September 1996, September 2007 and September 2018, and the peak depth of the terrace occurred in September 1985, October 1989, October 1996, October 2007 and October 2018, while the corresponding peak rainfall occurred in July 1985, August 1989, July 1996, July 2007 and July 2018, respectively. Therefore, the groundwater dynamics of Xining showed the characteristics of lagging 0–4 months behind the rainfall.

It can be seen from the correlation between the monthly groundwater depth and the rainfall of this month and the previous month in

Table 2. Correlation between monthly groundwater depth and rainfall.

Subarea	Groundwater Depth in Washland (m)		Groundwater Depth in Terrace (m)
Month	This Month	Last Month	This Month	Last Month
1	0.357	0.384	−0.394	0.435
2	0.338	0.354	0.345	0.385
3	−0.541	0.392	0.347	0.443
4	−0.573	0.488	−0.464	0.424
5	0.562	0.552	0.478	−0.528
6	−0.751	−0.646	0.573	−0.587
7	−0.722	−0.678	−0.719	−0.795
8	−0.867	−0.717	−0.745	−0.809
9	−0.878	−0.824	−0.817	−0.884
10	−0.528	−0.404	0.408	−0.740
11	0.419	0.496	0.367	0.464
12	0.369	0.379	0.342	0.399

that there is a significant negative correlation between the groundwater depth and the rainfall from July to September, while the correlation is weak from October to June of the next year, indicating that the groundwater depth is more sensitive to the response of rainfall in the rainy season. According to the correlation coefficients of different hydrogeological units, the groundwater depth in the floodplain (near the river zone) area has the strongest correlation with rainfall, and the response to rainfall is the most sensitive, followed by the terrace. In addition, the correlation between the groundwater depth and rainfall in the same month is stronger than that of the previous month in the flood plain area, while the correlation between the groundwater depth and rainfall of this month is weaker than that of the previous month in the terrace area, indicating that the lag effect of groundwater on rainfall in the terrace area is more obvious. The main reason is that the terraces are mostly hard impervious interfaces such as construction land or roads under the process of urbanization, which leads to the decrease and delay in effective rainfall recharge. In addition, the influence of human exploitation activities makes the groundwater depth increase, resulting in a thicker aeration zone and more water retention, resulting in a delayed groundwater recharge cycle. Therefore, the rainfall in the study area has a certain recharge lag in terms of the groundwater depth.

3.2.3. Response Process of Groundwater to Surface Runoff

Valley seepage recharge is also an important source of groundwater resources in Xining, especially near the river and floodplain area. River hydrological characteristics directly affect the dynamic change in groundwater. There is a close and frequent hydraulic relationship between groundwater and the Huangshui River in the valley plain area of Xining, and the aquifer lithology is dominated by sand, gravel and egg, with strong permeability. Therefore, surface runoff characteristics have a direct influence on groundwater dynamics. Figure 6 shows the relationship between groundwater depth and surface runoff. As can be seen from Figure 6a, the surface runoff from 1980 to 2020 increased, decreased, and then increased again from 13.6 m³/s in 1980 to 63.3 m³/s in 1989, decreased to 16.6 m³/s in 2001, and then increased to 61.0 m³/s in 2020. The variation trend of surface runoff and groundwater depth is basically consistent. From 1980 to 1989, the surface runoff continued to increase, but the groundwater depth continued to increase, indicating that the rainfall and surface water infiltration were insufficient to balance groundwater overexploitation. Although the discharge of the Huangshui River decreased from 1990 to 2000, the groundwater depth showed a decreasing trend due to the reduction in groundwater exploitation. From 2000 to 2020, the runoff continued to increase, and the groundwater depth also showed a small increasing trend. According to the columnar distribution of runoff anomalies, surface water is relatively rich in rainy years, and the amount of surface leakage recharge increases accordingly. Due to the large extraction amount, the groundwater depth continues to decrease, but the increase in surface runoff slows down the increase in the groundwater depth.

As can be seen from Figure 6b, the groundwater depth and discharge-to-exploitation ratio in the Huangshui River from 1980 to 2020 are negatively correlated. When the extraction amount is constant, the groundwater depth presents a power function decreasing trend with the increase in runoff, and the decreasing speed slows down. When the ratio of runoff to extraction is small, the slope of the power function curve is steep, indicating that the influence of surface runoff on groundwater dynamics is greater in dry years than in rainy years. When the runoff-to-extraction ratio increases, the scattered points of the groundwater depth in the figure increasingly deviate from the fitting curve, indicating that the evolution of groundwater becomes more complex with the increase in surface runoff. It is worth noting that the scatter of the surface runoff ratio is more dispersed than that of the rainfall ratio in Figure 6b, and the slope of its power function is lower than that of rainfall, indicating that surface runoff is also one of the important factors affecting the dynamic change in groundwater in Xining, but the influence degree is weaker than that of rainfall. With the intensification of human exploitation, the response degree of groundwater depth to rainfall and surface runoff gradually weakened.

It can be seen from the surface runoff data from January 1980 to December 2020 in Figure 7 that the flow of the Huangshui River presents significant seasonal variation characteristics. The high flow period occurred in from July to September every year with abundant water, accounting for 45.1–75.8% of the annual runoff, and the maximum discharge value was 121.5 m³/s in August 2018. The low flow period from October to June of the next year was small and stable, and the minimum flow was 3.4 m³/s in April 1980. At the same time, it can be seen from the groundwater depth curve that the groundwater depth of the floodplain and terrace also presents an annual periodic change rule, and decreases with the increase in discharge. It is worth noting that the dynamic variation of groundwater depth presents a time-lag characteristic, obviously lagging behind the variation curve of runoff of the Huangshui River. For example, the peak groundwater depth of the floodplain occurred in August 1983, September 1989, September 1997, September 2013 and September 2018, and the peak groundwater depth of the terrace occurred in September 1983, October 1989, October 1997, October 2013 and October 2018. The corresponding peak surface flow was in July 1983, July 1989, July 1997, July 2013 and July 2018. Therefore, the groundwater dynamics in Xining showed a synchronization that lagged behind the surface runoff for 1–4 months. According to the comparison between the groundwater depth curves of the floodplain and terrace in the figure, the annual maximum drop in groundwater depth can reach 3.8 m in the floodplain monitoring wells and 2.6 m in the terrace monitoring wells far from the river, indicating that the farther the distance from the river, the smaller the fluctuation in the groundwater depth with the exploitation amount, and the more delayed the change time.

3.3. Periodic Evolution of Groundwater and Main Driving Factors

3.3.1. Variation on Cone of Depression

It can be seen from the distribution of the real-part wavelet contours in Figure 8a that the groundwater depth in Xining has the phenomenon of different quasi-periodic oscillations at various time scales. There are three evident time scales for the evolution of groundwater depth dynamics from 1980 to 2020 with 5–7a, 9–14a and 17–25a: (1) The groundwater depth at 5–7a and smaller time scales fluctuates violently, and the periodic changes are local. Before 2010, the rule of alternating change was relatively clear. After that, the regular change disappeared, the periodic oscillation gradually became smooth, and the signal performance was messy and unstable. (2) The time scale from 9–14a is uniformly distributed in the time domain and has a significant global feature, with three complete alternating changes of positive/negative value. The three periods of depth reduction are 1980–1986, 1993–1999 and 2005–2011, respectively, with the gravity oscillating center of 1983, 1996 and 2008. The three increasing periods of groundwater depth were 1986–1993, 1999–2005 and 2011–2017 respectively, and the epicenters of shocks were 1990, 2002 and 2014, respectively. The negative dotted line is not completely closed after 2017, indicating that the groundwater depth was in the decreasing period from 2017 to 2020. Therefore, there is a high probability that the groundwater level will rise in the short term in the future. (3) There were two regional alternations between positive and negative values on the 17–25-year time scale. The two decreasing periods of groundwater depth were 1980–1985 and 1996–2007, and the oscillating center of gravity was 1980 and 2002, respectively. The two increasing periods of groundwater depth were from 1985 to 1997 and from 2007 to 2017, and the oscillating barycenter was distributed in 1992 and 2012, respectively. It is worth noting that the periodic evolution law of 17–25a after 2017 is consistent with that of 9–14a. The groundwater depth was in a decreasing period from 2017 to 2020, and the negative isoline did not completely close until 2020, which fully indicates that no matter what kind of time scale, the short-term inland water level will continue to rise in the future. To sum up, the evolution of groundwater depth in Xining has local changes in time domain and multi-level time scale structure characteristics. The groundwater depth in the study area does not exist in a specific and invariable periodic scale. However, with different time scales, the variation cycle of the groundwater depth also changes correspondingly, which shows that large time scales and small time scales are directly nested and contain each other. The decrease or increase in groundwater depth in several 9–14a time scales is mainly reflected in the rise or fall of the groundwater level in 17–25a time scale.

According to Formula (3), the wavelet square difference graph is solved and drawn. The periodic scale corresponding to each variance peak represents the main period of the time scale, as shown in Figure 8b. As can be seen from the figure, there are two significant peaks in the time series, among which 12a corresponds to the maximum peak, indicating that the groundwater depth signal oscillation is the strongest in the time scale of 9–14a; 21a corresponds to the second peak value, indicating that the signal oscillation at the time scale of 17–25a is the second; 6a corresponds to the third peak, with the weakest signal oscillation and negligible periodic effect. Therefore, the main cycle of controlling the periodic evolution process of groundwater dynamics in Xining from 1980 to 2020 is 12a, and because the shock intensity of the first peak is much greater than that of the second and other peaks, the first major cycle of groundwater depth is 12a. The 12a main cycle structure of the groundwater dynamic in Xining is mainly the superimposed effect of the geological tectonic environment, hydrogeological conditions and recharge and discharge relationship of the Tibetan Plateau.

3.3.2. Periodic Evolution of Key Drivers

The Morlet wavelet transform method in Section 3.3.1 was used to calculate the wavelet coefficients of rainfall and surface runoff, and the real-part contour map of the wavelet was drawn, as shown in Figure 9. As can be seen from Figure 9a, there are two obvious oscillating cycles of rainfall in Xining on different time scales with 9–14a and 17–25a. In the two cycles, there were two periods of positive/negative value alternating circulation in 17–25a and four periods of positive/negative value alternating circulation in 9–14a, with complete isoline closure and significant global characteristics. It can be seen from Figure 9b that the signal oscillation intensity on the scale of 9–14a is the largest, which is far greater than that on the scale of 17–25a. Therefore, the rainfall shows the interannual variation with the first main cycle of 12a. In addition, the isoline of wavelet coefficients in 2020 corresponding to scales 9–14a in Figure 9a is not completely closed, indicating that the future short-term rainfall is in an increasing trend. According to Figure 10a, the surface runoff has a time scale of 9–14a and 17–25a, and the wavelet coefficients of the two scales have two and four times the regional alternating changes of positive/negative value, respectively. It can be determined from Figure 10b that the peak value of the maximum signal oscillation intensity corresponds to 12a, indicating that the first main cycle of surface runoff is 12a. Therefore, the isolines corresponding to scales 9–14a in Figure 10a are not fully closed, indicating that the increasing trend of surface runoff will continue until after 2020. In conclusion, rainfall and surface runoff have the same periodic evolution process.

3.4. Time Delay Effect of Groundwater and Driving Factors

3.4.1. Continuous Wavelet Analysis of Groundwater and Main Driving Factors

The study of hydrological problems is not only limited to the analysis of multi-scale periodic changes of a single element, but also explores the correlation between different elements under multiple scales. In this section, continuous wavelet transform is used to analyze the periodic evolution characteristics of rainfall, surface runoff and groundwater depth in Xining from 1983 to 2020, as shown in Figure 11. The V-shaped black contour in the figure represents the 95% significance level and the black curve is the wavelet influence cone (COI). The region inside the V-shaped curve (outside the COI) is the valid spectral value, indicating that the oscillation time scale has passed the standard spectrum test of red noise at the significance level of 0.05, while the region outside the curve (inside the COI) represents the invalid spectral value [44]. According to Figure 11a–c, rainfall, surface runoff and groundwater depth within the V-shaped curve range all have a continuous main oscillation period of 12a with strong cohesiveness. It can be seen from

Table 3. Statistical information of continuous wavelet transform of rainfall, runoff and groundwater.

Sequence	Red Noise Test Period	Period Outside the COI	Main Oscillation Period/a
Precipitation	1983.2–2020.10	1984.5–2018.7	9–14
Runoff	1983.2–2020.10	1984.5–1990.10 1992.8–1996.2 1998.3–2001.1 2004.9–2018.7	9–14
Groundwater depth	1983.2–2020.10	1984.5–1996.11 1998.3–2003.10 2013.4–2018.7	9–14

that, in the time scale of 9–14a, the continuity of rainfall is the strongest among the three sets of sequences, and the distribution of surface runoff and groundwater depth is relatively dispersed. Among them, the whole band of rainfall passes the 95% red noise test, and the distribution outside the COI period is concentrated and has good continuity. The 95% red noise test of surface runoff in 1984–1990, 1992–1996, 1998–2001 and 2004–2018 has strong continuity, and the other periods are scattered. The groundwater depth in the three periods (1984–1996, 1998–2003 and 2013–2018) has strong continuity. Therefore, the groundwater depth in the same section is closely related to the hydraulic power of rainfall and surface runoff.

3.4.2. Cross-Correlation Analysis between Groundwater and Driving Factors

The groundwater depth of Xining was set as the input signal Xt, and the rainfall and runoff were set as the output signal Yt. Cross-correlation analysis was conducted on the rainfall, surface runoff and groundwater depth from January 1980 to December 2020, as shown in Figure 12. Formula (1) was used to calculate the time delay and correlation degree of groundwater depth to rainfall and runoff, with the maximum time delay set as 20 months. As can be seen from Figure 12a, the cross-correlation coefficient between rainfall and groundwater depth presents a 12-month periodic variation rule, indicating that the periodic variation of rainfall in a hydrological year also determines the variation characteristics of groundwater depth. According to the time delay that passes the black dotted 95% confidence interval test in Figure 12a, there is a time lag of 1–4 months between rainfall and groundwater depth, and the cross-correlation coefficient (0.297) reaches the maximum at the third month. It can be seen from Figure 12b that the cross-correlation coefficient between runoff and groundwater depth also presents a 12-month cyclical change, but the volatility and correlation are weaker than that of rainfall. The results show that the periodic variation of surface runoff is less than that of rainfall and the intensity of influence on groundwater depth is weaker than that of rainfall. According to the time delay of the 95% confidence interval test illustrated in Figure 12b, there is a time lag of 0–4 months between runoff and groundwater depth, and the cross-correlation coefficient (0.223) is the largest in the second month. Therefore, the lag time of groundwater depth to rainfall is slightly larger than that of surface runoff.

3.4.3. Continuous Wavelet Analysis of Groundwater and Driving Factors

The cross-wavelet transform was used to further verify the multi-scale correlation and time delay effect of rainfall, surface runoff and groundwater depth. The cross-wavelet transform nephogram and wavelet coherence nephogram of rainfall, surface runoff and groundwater depth from 1980 to 2020 are calculated using Formulas (4) and (6), as shown in Figure 13 and Figure 14. The arrow in the figure represents the phase relationship between the two groups of sequences, and the arrow to the left indicates that the two groups of signals are opposite, indicating that there is a negative correlation between the sequences. The right arrow indicates the same signal, indicating that there is a positive correlation between the sequences; the up arrow indicates that the preceding sequence is ahead of the following sequence; and the down arrow indicates that the preceding sequence lags behind the following sequence.

As can be seen from Figure 13a, in the frequency domain, in the high energy spectrum region, the cross-wavelet signal intensity at the time scale of 9–14a is the largest, the resonance relationship is significant and the arrow is to the left, indicating that there is a strong negative correlation between rainfall and groundwater depth, and it has passed the 0.05 significance level test. A complete period of the cross-phase angle (2π) represents 12a, and the phase difference arrow is about 8° to the left of the front, indicating that the rainfall is 1/45 cycle (0.27a) ahead of the groundwater depth, that is, 3.2 months. It can be seen from Figure 13b that in the low-energy spectrum region, the signal intensity at the time scale of 9–14a is the largest and has a negative correlation, with a correlation coefficient of more than 0.9. Therefore, there is a negative correlation between groundwater depth and rainfall in both high-energy and low-energy spectrum regions, and the groundwater depth lags behind the rainfall by 3 months.

As can be seen from Figure 14a, in the high energy spectrum region, the signal oscillation intensity reached its maximum and the arrow turned to the left at the time scale of 9–14a, indicating that the discharge of the Huangshui River was negatively correlated with the groundwater depth, and passed the significance level test of 0.05. The phase difference arrow is about 6° to the left, indicating that the Huangshui River discharge is 1/60 cycle (0.2a) ahead of the groundwater depth, namely, 2.4 months. As can be seen from Figure 14b, in the low energy spectrum region, the signal intensity reaches its maximum and presents a negative correlation at the time scale of 9–14a, with a correlation coefficient of more than 0.85. Therefore, there is a negative correlation between the groundwater depth and the discharge of the Huangshui River in both high-energy and low-energy spectrum regions, and the groundwater depth lags behind the surface runoff by 2 months.

It is worth noting that in the time domain, the high-energy and low-energy spectra of groundwater depth and rainfall and groundwater depth and Huangshui River discharge passed the red noise test from February 1983 to October 2020, and were continuously distributed outside of the COI period from May 1984 to July 2018. The lag time of groundwater depth on rainfall and surface runoff is basically the same, and the lag time of rainfall is slightly larger than that of runoff in the Huangshui River. The main reason is that there is little change in the hydrogeological conditions of the groundwater directly charged by river valley infiltration, which is only affected by the lithology, geological structure, water richness of the aquifer and other factors. The rainfall infiltration groundwater is affected by the lithology, geological structure, human activities and other comprehensive factors. With the continuous progress of urbanization, the hard impervious interface such as the construction area of Xining city continues to increase, and the way of rainfall infiltration recharge groundwater is lengthened and the lag time is prolonged, so the response speed of groundwater to rainfall becomes slower. Despite the intensification of human activities, rainfall and surface runoff are still the main driving forces affecting groundwater change in Xining under natural conditions.

4. Conclusions

In this study, cross-correlation analysis, wavelet analysis and cross-wavelet transform were used to explore the response mechanism and time delay effect of groundwater, rainfall and surface runoff in the northeastern margin of the Tibetan Plateau. The groundwater depth increased with the increase in extraction from 1980 to 2020, which was positively correlated with no time lag. The groundwater depth was negatively correlated with the ratio of rainfall to exploitation and the ratio of surface runoff to exploitation. With the continuous progress of urbanization, the groundwater in Xining is greatly affected by human exploitation activities, and the response degree of groundwater dynamics to rainfall and surface runoff decreases gradually with a time lag. The cyclic evolution process of driving factors was analyzed using complex Morlet wavelet transform. It was concluded that both the rainfall and the discharge of the Huangshui River had the same time scales of 9–14a and 17–25a, and the first main cycle was 12a, which showed an increasing trend for the future. The time delay effect between groundwater and the key driving factors was revealed by continuous wavelet and cross-wavelet transform. On the annual scale, the groundwater depth, rainfall and runoff of the Huangshui River have a strong condensate 12a continuous main oscillation cycle. On the monthly scale, the groundwater depth lagged behind the rainfall and the runoff of the Huangshui River for 3 months and 2 months, respectively.

The above conclusions prove that the cross-wavelet approach is suitable for applications related to groundwater in the northeast margin of the Qinghai–Tibet Plateau. However, due to the limitations of the length of the paper and the availability of monitoring data, this study only carried out groundwater research on key representative wells. In future studies, we will continue explore in detail the spatial difference in the time lag between groundwater and influencing factors.