*4.3. Humidity Index*

The potential evapotranspiration (PET) data were obtained from the "1 km monthly potential evapotranspiration dataset in China (1990–2021)" [42], based on the Hargreaves method [43]. The annual average PET at the four stations did not vary greatly, and ranges from 1190 mm at W1, to 1214 mm at W2, to 1187 mm at W3, and 1069 mm at W4. Seasonal variations in PET have also been observed (Figure S3). Generally, the variation pattern of PET is consistent with that of precipitation, though sometimes the peak of PET arrives ~1 month before the precipitation.

The humidity index (the aridity index in UNEP, [44]) can be used to assess the surfacewater stress or the arid degree at a given location [45]. It is defined as the ratio of precipitation P to the PET:

$$\text{HI} = \frac{\text{P}}{\text{PET}} \tag{8}$$

The results from the cross wavelet spectra between HI and GWLs are given in Table 6 and Figure S4. Compared to the time lags between P and GWLs as shown in Table 4, the lags between HI and GWLs are shortened by 14% at W1, 36% at W2, 28.5% at W3, and 19.6% at W4, respectively. This can be understood since not all rainfall can recharge aquifers. Only the effective rainfall can produce surface runoff or infiltrate into the subsurface. Compared to the precipitation, the humidity index reflects comprehensive influences of both P and PET, and to some extent the effective rainfall. Therefore, the GWLs response more quickly to the HI index than to the rainfall.

**Table 6.** Time lags from the cross wavelet spectra between HI and GWLs (period of 365 days band).


#### *4.4. Comprehensive Analysis*

To further investigate the controlling factors of the time lags between aquifers and rainfall identified using the XWT method, changes in lags induced by different factors are given in Table 7. It can be seen that variations in the time lags at W2 and W4 are most sensitive to pumping. In fact, Figure 3 and Table 4 indicate that compared to W1 and W3, W2 and W4 are in areas less affected by human activities with few pumping years. Therefore, they are very sensitive to the pumping under the regional background of groundwater level recession. The large percentage change in time lags induced by pumping at these areas highlights the fragility of local aquifer systems.

On the other hand, variations in the time lags at W1 and W3 are less affected by pumping, but most sensitive to HI. They are not sensitive to pumping because things could hardly become worse at these areas that have suffered from pumping to a certain extent. These unhealthy aquifers that have been affected by pumping, or these somewhat thirsty aquifers, will try to recover as long as they are replenished by rainfall. Therefore, they are most sensitive to HI. This highlights the resilience of local aquifers.


**Table 7.** Variation in time lags due to different factors.

Assuming that all rises in water level are due to recharge from precipitation, the water table fluctuation (WTF) method [46,47] tell us the recharge ratio (*α*) can be estimated as:

$$\alpha = S\_{\mathcal{Y}} \frac{\sum h}{\sum P} \tag{9}$$

where *Sy* is the specific yield; Σ*h* is cumulative rise in water-level; and Σ*P* is the total precipitation in the period corresponding to the water level rise.

Here, the hydrologic year during which the water table rose most significantly is chosen, and the maximum empirical values of *Sy* provided by [48] are also used. As a result, the recharge ratio we estimated as shown in Table 8 is larger than or close to the maximum value of the empirical values. As a whole, the recharge ratio decreases from upstream to downstream. Specially, it decreases from W2 to W4, which is in line with the inertia ranks as shown in Figure 4. Meanwhile, considering the smallest time lags indicated by the cross wavelet spectra, the site of W2 has a good potentiality for groundwater recharge.

**Table 8.** Recharge ratio calculated using the WTF method.


#### *4.5. Limitations*

There are many factors affecting the groundwater response time. Here, we only discussed the influences of rainfall intensity, evaporation, and groundwater pumping. Other factors, such as river water, were not taken into account. In fact, these factors are not isolated, but interact in various geological and geographical contexts. It is reported that the partial wavelet coherency method can detect localized and scale-specific bivariate relationships between predictor and response variables after removing the impact of other variables [30,49,50]. Therefore, in future research, the partial wavelet coherence method can be further implemented to distinguish the impacts of climate change and human activities on the GWLs, for better understanding their impacts on the groundwater flow system.

#### **5. Conclusions**

This research provides analyses of the long-term dynamics, the persistence, and the periodicity of shallow groundwater located in the Heilonggang region, China. The interrelation between precipitation and groundwater levels is also examined based on correlation and spectral analyses. The results of this research provide a more complete understanding of the local groundwater circulation system. The major conclusions are as follows.

Firstly, trend analysis shows that the GWLs at W1 and W2 present a downward trend while the GWLs at W3 and W4 show an upward trend over the observation period. Therefore, more attention should be paid to the upstream of the aquifer system.

Secondly, auto-correlation analysis indicates a rising trend in the memory time for aquifers from upstream to downstream. The cross-correlation analysis stresses that the shorter the memory time, the greater the correlation coefficient between rainfall and GWLs. The continuous wavelet spectra shows that the dominant period increases gradually from

2.1 years at W1 to 3.7 years at W4, further demonstrating longer regulation times from upstream to downstream.

Thirdly, both the cross wavelet spectra and the sliding-window cross-correlation method display that wells two and four respond quickly while wells one and three respond slowly to the local rainfall. The time lags between aquifers and rainfall at W1, W2, W3, and W4 are 139.14 ± 59.76 days (2008–2020), 23.27 ± 12.03 days (2005–2014), 145.01 ± 68.00 days (2007–2020), and 59.22 ± 26.14 days (2005–2019), respectively. The temporal lags of groundwater to precipitation are shortened by 10~13 days during the wet year conditions, and the lags during the pumping years are 1.3~2.0 times those during the years without pumping. The time lags between HI and GWLs are reduced by 14~36%, compared to those between rainfall and GWLs.

Further analysis shows that variations in the time lags at W2 and W4 are most sensitive to pumping, while variations in the time lags at W1 and W3 are less affected by pumping, but most sensitive to HI. The overestimated recharge ratio decreases from 0.32 at W2 to 0.17 at W4, suggesting that the site of W2 has a good potentiality for groundwater recharge.

Finally, although this research provides a new insight into relations between precipitation and groundwater in the study area, there are still some other factors, such as river water, which were not considered. Future researchers could possibly use the partial wavelet coherence method to distinguish the effects of human activities and climate change on the GWLs.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/w15061100/s1, Figure S1. Aquifer groups in the Heilonggang region. Figure S2. Shallow groundwater flow field in (a) 1959 and (b) 2005. Figure S3. Monthly precipitation (P) and monthly potential evapotranspiration (PET) at (a) W1, (b) W2, (c) W3 and (d) W4. Figure S4. Cross wavelet spectra between HI and GWLs at (a) W1, (b) W2 (c) W3 and (d) W4.

**Author Contributions:** Methodology, C.W.; formal analysis, C.W.; investigation, H.L.; data curation, Y.L.; writing—original draft preparation, C.W. and W.Q.; Writing—Review & Editing, W.Q.; funding acquisition, F.D.; visualization, Y.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by the S&T Program of Hebei (No: 22373601D), the Water Conservancy Science and Technology Project of Hebei Province, China (2021-39), the Open Fund for Hebei Province Collaborative Innovation Center for Sustainable Utilization of Water Resources and Optimization of Industrial Structure (XTZX202111), the Scientific Research Projects of the Higher University in Hebei (BJK2022007), the Natural Science funds Project in Hebei Province (D2020403022, E2021403001), National Pre-research Funds of Hebei GEO University in 2023 (KY202314), the Youth Project of Hebei GEO University (QN202118 and QN202141).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data are contained within the article.

**Acknowledgments:** We are very grateful to the three anonymous reviewers for their constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
