**1. Introduction**

With the increase in extreme weather events induced by climate change, understanding the response mechanism of groundwater systems to precipitation is of critical importance [1–4]. To date, two valuable correlation methods including the auto-correlation and cross-correlation function have been typically applied to assess the responses of aquifers to precipitation [5–7]. The auto-correlation analysis characterizes the degree of "inertia" [8,9], or "persistence" [3], or "memory effect" [10] of an individual time series. The cross-correlation analysis can provide the response time between rainfall and groundwater levels (GWLs). In addition, the correlation methods are usually combined with spectral analysis, such as Fourier analysis, to detect the periodicity of a signal.

**Citation:** Wang, C.; Dai, F.; Liu, Y.; Wang, Y.; Li, H.; Qu, W. Shallow Groundwater Responses to Rainfall Based on Correlation and Spectral Analyses in the Heilonggang Region, China. *Water* **2023**, *15*, 1100. https:// doi.org/10.3390/w15061100

Academic Editors: Alban Kuriqi and Luis Garrote

Received: 14 February 2023 Revised: 7 March 2023 Accepted: 9 March 2023 Published: 13 March 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Jenkins and Watts [11] and Box and Jenkins [12] explained the basic principles of correlation and spectral analyses. Mangin [13] applied these methods to three karstic systems in the Pyrenean Mountains (France–Spain), and their work indicated that these analyses provide an excellent method for the investigation of a hydrological system. After that, correlation and spectral analyses have been widely used in different aquifer systems, including karstic systems [14–16], alluvial aquifers [3,5], and coastal aquifers [17,18], etc.

Although correlation and spectral analysis have been widely used in hydrology and hydrogeology, these methods cannot describe how the frequency of a signal changes with time. Instead, the wavelet transform method is an effective tool for detecting the periodicity of a nonlinear system and can provide localized intermittent periodicities [19–21]. Based on the cross wavelet and wavelet coherence methods, Yu and Lin [22] found that the temporal lag from rainfall to groundwater was about 3.71–72.07 days for the Pingtung Plain in Taiwan. Zhang et al. [23] found that the time lag of groundwater depth to precipitation in the Yellow River Delta during 2006–2010 ranged from 35. 51 to 178. 36 days, and the relationship between groundwater depth and precipitation is largely affected by land use types, soil texture, and micro-geomorphic types. Cai et al. [24] show that the response time of groundwater levels to rainfall during 2006–2018 extended from 80 to ~190 days in Puyang, Henan, China, and that it increases with the burial depth of groundwater. Although there are many studies investigating the relationship between rainfall and aquifers, seldom does research focus on areas facing severe water shortage such as the Heilonggang region, China.

The Heilonggang region is located in the southeast of Hebei Province (Figure 1), a part of the North China Plain, and serves as an important agricultural planting area for the nation. Agricultural water accounts for 76% of water resource consumption in this area, and more than 80% of agricultural irrigation water comes from groundwater [25]. The increasing demand for groundwater resources has led to the decline in GWLs and caused a series of ecological problems such as land subsidence, etc. Previous studies of this area mainly focus on the optimization of irrigation and planting regimes [26–28], and few researchers have paid attention to the long-term groundwater dynamics and associated responses to rainfall in this region. Clarifying the relationship between groundwater and precipitation can help us better understand the local groundwater circulation, provide references for groundwater resource management, and be conducive to ecological protection [29,30]. Therefore, in this study, we aim to (1) estimate the persistence of aquifers and their responses to rainfall through the auto-correlation and cross-correlation functions, (2) identify groundwater periodicity through spectral analyses, (3) determine the groundwater response time through cross-correlation and cross spectral analyses, and (4) explore the influences of rainfall intensity, humidity index, and groundwater pumping on the response time.

**Figure 1.** Location of the study area.

#### **2. Materials and Methods**

#### *2.1. Study Site*

The Heilonggang region covers an area of 34,700 km2. The terrain in this region is gentle and inclines slightly from southwest to northeast with a topographic slope of 0.2~0.1‰. From west to east, the geomorphic type is dominated by the mountain alluvial-diluvial plain, the central alluvial-lacustrine plain and the coastal plain, respectively. Accordingly, from west to east, sediments change from gravel in front of the mountain, to medium coarse and medium fine sand in the middle, then to fine sand in the coastal area.

The study site is mainly affected by the warm temperate semi-arid and semi-humid continental monsoon with four distinct seasons. The annual average precipitation is 500~600 mm, mostly concentrated in July and August. The river system in the Heilonggang area comprises the Zhangweinan Canal system, the Heilonggang Canal system, and the Ziya River system. All rivers finally drain into the Bohai Sea (Figure 1). Reservoirs have been built at the upper reaches of almost all rivers, and more than 80% of the surface runoff is detained by the reservoirs, resulting in long-term drying of the rivers flowing through the area and greatly reducing the amount of river leakage.

Aquifers in the Heilonggang area are mainly composed of Quaternary strata. The Quaternary aquifer system can be divided into four aquifer groups from top to bottom (Figure S1). Among them, the first aquifer group is the aquifer we focused on in this study. It is composed of Holocene and Upper Pleistocene loose sediments, with a bottom depth of 20~50 m [31]. The vadose zone lithology in this area is mainly composed of silt and silty clay (Figure 2).

**Figure 2.** Vadose zone lithology and the flow direction map (2020).

In the 1950s, the degree of groundwater exploitation in this area was very low, and the whole groundwater flow system was in a natural state. Shallow groundwater generally flows from southwest to northeast (Figure S2a). Under the influence of human activities, local groundwater depression cones have formed in the north part of the region (Figure S2b). At present, the groundwater depression cones have further expanded towards the west (Figure 2), and the southwest–northeast hydraulic gradient has decreased greatly compared to the natural state.

Shallow fresh groundwater mainly occurs in the west and north of the research area, and shallow brackish water mainly occurs in the east of the region. The distribution area of the shallow fresh groundwater accounts for 51.6% of the total Heilonggang area [32]. Rainfall serves as the main source for the river and shallow groundwater.

#### *2.2. Datasets*

The groundwater data were obtained from the "Groundwater Almanacs of Geological Environment Monitoring in China". Four boreholes monitoring a substantial length of shallow groundwater were selected in this study, as shown in Figure 1. Detail characteristics of each monitoring well are presented in Table 1. The rainfall data before 2018 was extracted daily from the "China meteorological forcing dataset (1979–2018)" [33], and the data after 2018 was collected from China Meteorological Administration (http://data.cma.cn/, accessed on 1 June 2022).

**Table 1.** Detail information of the monitoring wells. The thickness of unsaturated zones was calculated as the average depth over the observation period.


#### *2.3. Mann–Kendall Trend Analysis*

The Mann–Kendall trend test is a nonparametric statistical test [34,35]. It is not necessary to assume that samples obey a certain distribution and are not disturbed by a few outliers.

The trend of time series is determined using *Z* values:

$$Z = \begin{cases} (S-1) / \sqrt{n(n-1)(2n+5)/18} & \text{ $S > 0$ } \\ 0 & \text{ $S = 0$ ,} \\ (S+1) / \sqrt{n(n-1)(2n+5)/18} & \text{ $S < 0$ } \end{cases} \tag{1}$$

where *S* is the testing statistic, and *n* is the sample size. A positive value of *S* indicates an increasing trend, and vice versa. If |*Z*| ≥ 1.96, it indicates the time series passes the significance test with 95% confidence.

#### *2.4. Auto-Correlation and Cross-Correlation Functions*

The auto-correlation functions were used to depict the persistence degree of the time series. The auto-correlation coefficient *rxx*(*k*) is expressed as [12]:

$$r\_{xx}(k) = \frac{\frac{1}{N} \sum\_{i=1}^{N-k} (\mathbf{x}\_i - \overline{\mathbf{x}})(\mathbf{x}\_{i+k} - \overline{\mathbf{x}})}{\sigma\_x^2} = \frac{\sum\_{i=1}^{N-k} (\mathbf{x}\_i - \overline{\mathbf{x}})(\mathbf{x}\_{i+k} - \overline{\mathbf{x}})}{\sum\_{i=1}^{N} (\mathbf{x}\_i - \overline{\mathbf{x}})^2},\tag{2}$$

where *N* is the length of the series; *k* is the time lag; and *x* is the arithmetic mean of the series. The cross-correlation coefficient between series *x* and *y* is defined as [12]:

$$r\_{xy}(k) = \frac{\frac{1}{N} \sum\_{i=1}^{N-k} (\mathbf{x}\_i - \overline{\mathbf{x}})(y\_{i+k} - \overline{\mathbf{y}})}{\sigma\_x \sigma\_y} = \frac{\sum\_{i=1}^{N-k} (\mathbf{x}\_i - \overline{\mathbf{x}})(y\_{i+k} - \overline{\mathbf{y}})}{\sqrt{\sum\_{i=1}^{N} (\mathbf{x}\_i - \overline{\mathbf{x}})^2 \sum\_{i=1}^{N} (y\_i - \overline{\mathbf{y}})^2}},\tag{3}$$

where *N* is the length of the series; *k* is the time lag; and *σ<sup>x</sup>* and *σ<sup>y</sup>* are the standard deviations.

The cross-correlation functions characterize the relationships between the input and output signals. Here, we take the rainfall as input signals and the GWLs as output signals. If the rainfall can be considered random, the cross-correlation functions give the impulse response of the aquifer.

The sliding-window cross-correlation method followed by Delbart et al. [2] was also carried out to investigate the influences of rainfall intensity on the response time. The sliding-window cross-correlation method consists of slicing the *x* (precipitation) and *y* (GWLs) series with partially superposed windows. For each window, the *rxy* values between rainfall and GWLs are computed, and the corresponding response time is identified. Then, a time series of response times is obtained for different windows. In this study, the window length is 6 months, and the sliding interval is 1.5 months. Note that the correlation coefficients should be not lower than the standard error of ~2/√*N*, where *<sup>N</sup>* is the sample size, and "2" is the critical value for the 0.95 probability of the normal distribution. That is, values for which the *rxy* peak was not significant at a 95% confidence level were left out.

#### *2.5. Wavelet Transform*

Hydrologic time series are usually nonstationary with temporal variations in both frequency and amplitude. Through the continuous wavelet transform, a complete timescale representation of localized and transient phenomena occurring at different time scales can be obtained [36]. The wavelet spectrum is defined as the modulus of wavelet coefficients. This wavelet spectrum can also be averaged in time, known as the global averaged wavelet spectrum [20], which helps to identify the characteristic periods within a single time series. Here, the Morlet wavelet serves as the wavelet mother function [37].

After the continuous wavelet transform, the cross wavelet transform (XWT) is used to determine the cross wavelet spectrum of two time series, and to examine relationships between the two series. The XWT is defined as:

$$\mathcal{W}\_t^{XY}(\mathbf{s}) = \mathcal{W}\_t^X(\mathbf{s})\mathcal{W}\_t^{Y\*}(\mathbf{s}),\tag{4}$$

where *W<sup>X</sup> <sup>t</sup>* (*s*) is the wavelet transform of time series *xt* (rainfall) at frequency scale *s*; and *WY*<sup>∗</sup> *<sup>t</sup>* (*s*) is the complex conjugate of wavelet transform *W<sup>Y</sup> <sup>t</sup>* (*s*) for *yt* (GWLs). The XWT can be represented using polar coordinates:

$$\mathcal{W}\_t^{XY}(s) = \left| \mathcal{W}\_t^{XY}(s) \right| e^{t\rho\_l(s)},\tag{5}$$

where # #*WXY <sup>t</sup>* (*s*) # # is the power of the cross wavelet; and *<sup>ϕ</sup>t*(*s*) is the phase angle, which denotes the delay between the two series at time *t* and scale *s*.

The time lag Δ*T* at a scale *s* between the two signals is calculated by:

$$
\Delta T = T(s) \times \frac{q\_l(s)}{2\pi},\tag{6}
$$

where *T*(*s*) is the period relative to the scale *s*.

The distribution of the cross wavelet power is:

$$P(\frac{|\mathcal{W}\_t^{XY}(s)|}{\sigma\_X \sigma\_Y} < p) = \frac{Z\_\mathcal{V}(p)}{\nu} \sqrt{P\_k^X P\_{k\ \prime}^Y} \tag{7}$$

where *P<sup>X</sup> <sup>k</sup>* and *<sup>P</sup><sup>Y</sup> <sup>k</sup>* are Fourier background spectra of the two series *xt* and *yt*; and *Zν*(*p*) is the confidence level associated with the probability *p*, and *p* = \$ *<sup>Z</sup>ν*(*p*) <sup>0</sup> *fν*(*z*)*dz*. The 5% significance level is determined using *Z*2(95%) = 3.999.

This study performed XWT on the original time series of rainfall and GWLs to understand their relationships. To exclude the edge effects, the cone of influence is introduced for all wavelet transforms. Codes used for wavelet analyses were based on those originally written by Grinsted et al. [19] and finished by Matlab 2014b (MathWorks, Natick, MA, USA).

#### **3. Results**

#### *3.1. Observed Time Series and Trend Analysis*

Piezometric levels of all boreholes shown in Figure 3 characterize the behavior of unconfined aquifers. The average GWLs decrease from upstream to downstream, i.e., from 28.77 m at W1, to 21.90 m at W2, to 6.42 m at W3, and to 1.09 m at W4. The standard deviation of the GWLs also declines from 1.25 m at W1 in the west, to 0.70 m at W2 in the middle, and to ~0.5 m for the other two wells in the east of the region. According to the Pearson Type-III distribution curve, the wet year and dry year are identified and denoted by the blue and orange bands, respectively. For example, in Figure 3d, the years 2010 and 2014 are identified as the wet year and dry year, respectively.

**Figure 3.** Observed time series of rainfall and GWLs at (**a**) W1, (**b**) W2, (**c**) W3, and (**d**) W4. The blue and orange bands indicate the wet year and dry year, respectively.

According to the Mann–Kendall statistics (Table 2), the GWLs at W1 and W2 present a downward trend while the GWLs at W3 and W4 show an upward trend (all at the 95% confidence level). For the annual averaged GWLs, piezometric levels of W1 and W2 decrease at a rate of 0.11 m/y and 0.04 m/y, respectively. Annual averaged GWLs of W3 and W4 increase at a rate of 0.03 m/y and 0.01 m/y, respectively. Both W3 and W4 are located in Cangzhou City. In 2008, Cangzhou City was listed as a national pilot area for groundwater protection, and the amount of exploitation decreased rapidly. Yan et al. [38] also reported that after 2008, the amount of artificial mining decreased and the depth of groundwater decreased. Our study result is in line with them.


**Table 2.** Test results of the Mann–Kendall trend analysis.
