Characterizing Groundwater Level Response to Precipitation at Multiple Timescales in the Lubei Plain Region Using Transfer Function Analysis

Abstract

Groundwater is essential for ecosystem stability and climate adaptation, with precipitation variations directly affecting groundwater levels (GWLs). Human activities, particularly groundwater exploitation, disrupt the recharge mechanism and the regional water cycle. In this study, we propose a new research framework: On the basis of analyzing the spatiotemporal variability characteristics of precipitation and shallow GWL, we used transfer function analysis (TFA) to quantify the multi-timescale characteristics of precipitation–GWL response under the effects of climate change and human activities. In addition, we evaluated the GWL seasonality and seasonal response while also considering apportionment entropy. We applied this framework to the Lubei Plain (LBP), and the findings indicated the following: (1) Annual precipitation in the LBP decreased from southeast to northwest, with July and August contributing 51.5% of total rainfall; spatial autocorrelation of GWL was high and was influenced by geological conditions and cropland irrigation. (2) The coherence between GWL and precipitation was 0.96 in the high-precipitation areas but was only 0.6 in overexploited areas, and sandy soils enhanced the effective groundwater recharge, with a gain of 1.65 and a lag time of 2.1 months. (3) Over interannual scales, GWL response was driven by precipitation distribution and aquifer characteristics, while shorter timescales (4 months) were significantly affected by human activities, with a longer lag time in overexploited areas, which was nearly 60% longer than areas that were not overexploited. (4) Groundwater exploitation reduced the seasonality of GWL, and irrigation reduced the coherence between GWL and precipitation (0.5), with a gain of approximately 0.5, while a coherence of 0.8 and a gain of 3.5 were observed in the non-irrigation period. This study clarified the multi-timescale characteristics of the precipitation–GWL response, provided a new perspective for regional research on groundwater response issues, and proposed an important basis for the short-term regulation and sustainable development of water resources.

1. Introduction

Groundwater, the largest form of freshwater storage, supplies drinking water to nearly 50% of the global population and is an important guarantee for the sustainable development of ecosystems and human response to climate change [1,2]. Precipitation infiltration is a primary method for recharging groundwater [3], and factors such as precipitation amount and intensity significantly influence the response mechanisms of groundwater systems [4,5,6]. However, in recent decades, climate change and human activities have posed serious threats to groundwater resources [7]. Therefore, the ability to clarify the groundwater level (GWL) response of precipitation under human impact is vital for our understanding of regional water-cycle processes and to effectively manage groundwater resources [8,9].

Precipitation and its spatial distribution directly affect groundwater recharge, which in turn constrains groundwater-dependent ecosystems and human activities [10]. When precipitation occurs, a portion of the precipitation infiltrates into the aquifer. The process may be influenced by factors such as precipitation [11], aquifer characteristics [12], and human activities [13]. With climate change increasing the complexity of precipitation patterns [14], continuous groundwater exploitation and irrigation [6,15], as well as inter-basin water transfers [16], have further complicated the regional water cycle. Therefore, identifying the precipitation–GWL response process in a complex and changing environment has become particularly important. Because the GWL fluctuations are influenced by the combined effects of multiple factors, the GWL response to precipitation shows significant temporal variability, ranging from within a day to years [17]. Therefore, the exploration of the multi-timescale response provides important insight into the mechanisms driving the groundwater system response and adaptations to climate change [12,18,19,20].

At present, two main methods are used to study the response of groundwater to precipitation changes. The first method involves using numerical simulation models of groundwater flow, such as MODFLOW [21,22] and HYDRUS-1D [23], to simulate the response of GWL to precipitation. This method often requires a large amount of input data (e.g., groundwater recharge, soil properties, hydrogeologic conditions) to support it, making it difficult to implement for data-poor areas. The second method employs mathematical and statistical methods to analyze the response characteristics relation to long-term precipitation and GWL. For example, Asoka [24] investigated the correlation between precipitation and groundwater recharge across various regions in India using multiple linear regression and empirical orthogonal functions, Tashie [3] assessed the response of groundwater recharge to precipitation by calculating the ratio of groundwater recharge to precipitation, and Qiu [25] used cross-correlation function and cross-wavelet analysis to explore the lag effect of GWL on precipitation. These studies quantified the existence, strength, and lag time of the precipitation–groundwater response characteristics from various perspectives. These methods, however, analyze the response of the groundwater system only from a single perspective, making it difficult to capture the complex characteristics of multi-timescale processes. Meanwhile, human activities, such as irrigation [26,27,28] and groundwater exploitation [29], which influence groundwater recharge and its response to precipitation significantly, cannot be ignored.

Transfer function analysis (TFA) is a time-series analysis method used to understand the relationship between inputs and outputs of a system by simplifying complex physical processes [17]. It effectively characterizes the coherence, gain, and phase within the response system [30]. These three characteristic parameters combine the correlation of the linear regression, the response strength indicated by the ratio of groundwater recharge to precipitation, and the lag time calculated by the cross-correlation function. TFA addresses the issue of singularity present in the aforementioned approaches, allowing for a more comprehensive description of GWL response to precipitation. Additionally, TFA can decompose component processes in different timescales, enabling an analysis of the characteristics of the response of GWL to precipitation in multiple timescales. Unlike the standard groundwater index [12,31] and wavelet transform [32,33,34], TFA divides multiple timescales from the frequency domain perspective, clarifying inherent response characteristics throughout the entire time series [35]. This method can effectively capture the multi-timescale response of the groundwater system to precipitation and provide a novel perspective on the complexity of GWL fluctuations.

The Lubei Plain (LBP), situated in the southeastern North China Plain, is predominantly an agricultural area. Influenced by natural conditions, the total water resources are insufficient. Water resource utilization primarily relies on groundwater and the Yellow River. More than 40 years of continuous exploitation has significantly lowered the GWL, creating regional groundwater depression cones centered on Dezhou City, Binzhou City, and Dongying City, and leading to geological issues such as land subsidence [36,37,38,39,40,41] and seawater intrusion [42]. Following the groundwater restriction policy and the South-to-North Water Diversion Project, the pressure on groundwater resources has been alleviated, and the level has begun to rebound [43]. Consequently, the regional water cycle and the groundwater response to precipitation have exhibited new characteristics. Current research on the LBP, however, has featured two main aspects. First, scholars believe that the continuous exploitation of groundwater is the main reason for the decline of the groundwater level in the LBP, which has triggered regional land subsidence [37,38,40,41]. Secondly, the BP neural network [44] and control evaluation [45] methods were used to simulate and predict the groundwater management schemes of different levels in the LBP. These studies provide a lot of information for understanding the dynamic characteristics of the groundwater level and its evolution in the LBP. These studies have provided a great deal of information to support the understanding of the dynamic characteristics of the groundwater level and its evolution in the LBP, but they have not yet investigated the characteristics of the groundwater response to precipitation and the impact of groundwater exploitation and irrigation in the region on this response process. Given the complexity of the water cycle in the LBP, examining groundwater response to precipitation and the effects of exploitation and irrigation is crucial for the rational development and sustainable management of groundwater resources.

Therefore, focusing on the question of how shallow groundwater levels respond to precipitation under the dual influence of climate change and human activities, this paper proposes a new research framework: On the basis of an analysis of the spatial and temporal variability characteristics of regional precipitation and shallow GWL by using the Global Precipitation Measurement (GPM) dataset and shallow GWL monitoring data, the multi-timescale characteristics of the response of shallow groundwater to precipitation changes are quantified using the TFA method to divide different timescales in the frequency domain. Next, apportionment entropy (AE) and Moran’s Index are introduced to characterize the seasonal fluctuation in the shallow GWL and its spatial correlation. On this basis, a transfer function model is established for the irrigated and non-irrigated periods to investigate the impact of cropland irrigation on the seasonal response characteristics of groundwater. This research enhances our understanding of how GWLs respond to precipitation changes influenced by climate change and human activities, provides a new perspective for regional research on groundwater response issues, and proposes an important basis for the short-term regulation and sustainable development of water resources. The rest of the paper is organized as follows: Section 2 describes the study area, dataset, and the methods of semivariogram function model, TRA, apportionment entropy seasonality index, and Moran’s I index. Section 3 presents a detailed analysis of the results, while Section 4 presents a discussion of the results. Finally, the conclusions are summarized in Section 5.

2. Materials and Methods

2.1. Study Area

The Lubei Plain is located in the northwestern of Shandong Province, China, and the southeast section of the North China Plain. It is bordered by the Taishan Mountains to the south, lying between 115°20′E and 119°45′E, 35°47′N and 38°16′N, encompassing a total area of 43,400 km². The region includes four major cities: Dezhou City, Binzhou City, Dongying City, and Liaocheng City, as shown in Figure 1.

2.1.1. Climate

The region has a temperate monsoon climate with four distinct seasons: dry and windy in spring, hot and rainy in summer, irregular droughts and floods in autumn, and cold with minimal snowfall in winter. The average annual precipitation is 667.2 mm, which is characterized by significant variability throughout the year, with rainfall primarily concentrated in July and August, accounting for 52.6% of the annual total [46]. The average annual evaporation in the LBP ranges from 1000 to 1300 mm, peaking from March to June, when it accounts for more than 50% of the annual evaporation [46]. As a result of the uneven distribution of precipitation and evaporation, the region is prone to climatic characteristics, such as spring drought, summer flooding, and late autumn drought, leading to frequent occurrences of drought and flood disasters.

2.1.2. Hydrogeological Condition

The landform type of the LBP is primarily characterized by the Yellow River alluvial plain (Figure 2b). The north of the Yellow River is divided into three sections from west to east: the Yellow River alluvial plain, the Yellow River delta alluvial–marine plain, and the coastal plain. Conversely, the south is divided into the Piedmont alluvial–diluvial inclined plain, the Yellow River alluvial plain, the Yellow River delta alluvial–marine plain, and the coastal plain from southwest to northeast [47].

Groundwater is mainly stored in the loose sedimentary pore fissures of the Quaternary and Neoproterozoic Minghuazhen Formation. An aquifer group with a buried depth of less than 60 m is defined as a shallow aquifer. The thickness of the water-bearing sand layer is about 5 to 25 m, which is one of the main exploitation layers. A buried depth of 60 to 200 m is defined as a medium aquifer, and a depth of more than 200 m is defined as a deep aquifer [46]. The shallow aquifer is composed predominantly of fine sand and silt that has been deposited by the Yellow River since the Late Pleistocene, with additional silt formed during the lake and marine phases. This aquifer is recharged primarily by atmospheric precipitation, irrigation return flows from the Yellow River, and the river and lake leakage, serving as the main water source for agricultural irrigation [48]. Furthermore, the water cycle exhibits strong fluctuations and a rapid renewal rate, with a water level typically buried at depths of 3 to 20 m, which is favorable for groundwater recharge and exploitation [49].

2.1.3. Water Resources

The surface water system in the region is a part of the basins of the Haihe River, the Yellow River, and the coastal rivers. All of the rivers in this area are rain-fed, exhibiting a dramatic increase in flow during the rainy season and low flow in the dry season. Because of the natural conditions, the total water resources are insufficient, with a total of 4.219 billion m³ of water resources, including 3.786 billion m³ of groundwater resources and 1.872 billion m³ of surface water resources [50], and the per capita water resources possession only reaching 195 m³·a⁻¹. Water resource utilization primarily relies on the Yellow River and groundwater.

Since the 1970s, large-scale exploitation of groundwater has been implemented in the LBP, and groundwater exploitation was about 0.3 billion m³ in the 1970s. Rapid industrial and agricultural development has significantly increased the demand for water resources. In DeCheng district alone, the number of deep groundwater exploitation (DGE) wells rose sharply from 33 in 1973 to 278 in 2000. In 2009, the DGE increased to 3.41 billion m³. Long-term groundwater exploitation has caused the GWL of aquifers to drop sharply, forming a regional groundwater depression cone centered on Dezhou City, Binzhou City, and Dongying City. From 2000 to 2015, the depth of the middle and deep GWL increased by nearly 30 m. Following the implementation of comprehensive groundwater management policies by the Shandong Provincial Government [43], GWL began to rebound after 2015, leading to new features in the regional water cycle.

2.2. Dataset

2.2.1. Precipitation Datasets

Because of the limited number of continuously recorded meteorological stations, existing data are insufficient to adequately describe the spatial variability of precipitation in the LBP. Therefore, remote sensing precipitation datasets are introduced in this paper. Four commonly used precipitation datasets—Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS), Multi-Source Weighted-Ensemble Precipitation (MSWEP), Climate Prediction Center MORPHing technique (CMORPH) and Global Precipitation Measurement (GPM))—were selected calibrated using data from meteorological stations (Figure 3c). On the basis of the verification results, we ultimately selected the GPM dataset, which demonstrated a low root mean square error (RMSE), good performance during extreme precipitation years, and wide applicability, as the primary precipitation dataset for this study.

GPM is an international satellite network initiated by NASA and the Japan Aerospace Exploration Agency (JAXA) to measure precipitation. As a follow-on to the Tropical Rainfall Measuring Mission (TRMM), GPM integrates core observatory satellites equipped with advanced radar and radiometer systems to provide global precipitation data (https://gpm.nasa.gov/ (accessed on 31 January 2024)). GPM-IMERG is a combined product derived from mutual calibration, fusion, and interpolation of precipitation measurements from all microwave satellites [51]. In this study, we utilized GPM-IMERG monthly precipitation data from June 2000 to December 2020, with a spatial resolution of 0.1° × 0.1° (approximately 10 km × 10 km).

2.2.2. Groundwater Level

The GWL data were obtained from the geological survey department of Shandong Province, China. This study focuses on dynamic observation data from shallow groundwater monitoring wells, specifically from June 2000 to December 2020. Only data with continuous GWL monitoring for more than 10 years and missing data accounting for less than 10% of the time series were included. After quality checking and preprocessing, 46 shallow groundwater monitoring wells were selected for analysis. The distribution of those wells is shown in Figure 2a, and the details of 13 typical groundwater level monitoring wells are shown in

Table 1. The details of the typical groundwater level monitoring wells.

Well ID	Location	Monitoring Layer	Time Range
Well-01	Lingcheng, Dezhou City	Shallow	Jan. 2000–Dec. 2020
Well-02	Decheng, Dezhou City
Well-03	Decheng, Dezhou City
Well-04	Decheng, Dezhou City
Well-05	Lingcheng, Dezhou City
Well-06	Wucheng, Dezhou City
Well-07	Xiajin, Dezhou City
Well-08	Guangrao, Dongying City	Shallow	Jan. 2000–Dec. 2020
Well-09	Guangrao, Dongying City
Well-10	Guangrao, Dongying City
Well-11	Guangrao, Dongying City
Well-12	Guangrao, Dongying City
Well-13	Kenli, Dongying City

. Monthly precipitation data and monthly GWL monitoring data were used to analyze the relationship between precipitation and GWL in the LBP. To minimize errors caused by the distance between the meteorological stations and the GWL monitoring wells, we assigned meteorological data for each monitoring well based on the nearest GPM precipitation data grid unit [24].

Additionally, to illustrate the potential impact of human activities on GWL response to precipitation, this paper incorporates land use type data, groundwater overexploitation areas, annual shallow groundwater exploitation (SGE), and deep groundwater exploitation (DGE) statistics for each district and county in the LBP, the details are shown in Figure 4. The land use information is extracted from the China Land Cover Dataset (CLCD) [52], the groundwater overexploitation distribution data are from the notice published by the People’s Government of Shandong Province in 2012, and the groundwater exploitation volume data are from the water resources bulletins of various districts and counties.

2.3. Method

In this paper, a new research framework is established and applied to the LBP to analyze the multi-timescale characteristics of the response of the shallow GWL to precipitation changes under the influence of human activities, as shown in Figure 5.

First, the spatiotemporal variation characteristics of precipitation and shallow GWL in the LBP region are analyzed. Second, the multi-timescale characteristics of the response of the shallow GWL to precipitation under the influence of human activities are quantified by the transfer function analysis method, which divides different timescales in the frequency domain. Third, apportionment entropy and the Moran Index are introduced to characterize the seasonal fluctuation in the shallow GWL and its spatial correlation. On this basis, a transfer function model is established for the irrigated and non-irrigated periods to investigate the impact of cropland irrigation on the seasonal response characteristics.

2.3.1. Semivariogram Function

The semivariogram function reflects the spatial relation between the sampling point and its adjacent sampling point, and the spatial variability of each GWL monitoring point in the study area can be evaluated based on the semivariogram parameters [53]. The generalized formula of the semivariogram function is as follows:

$$r\left(h\right)={\displaystyle \frac{1}{2N\left(h\right)}}{\sum }_{i=1}^{N\left(h\right)}{\left[z\left(x_i\right)-z\left(x_i+h\right)\right]}^2,$$

(1)

where $r\left(h\right)$ is the semivariogram function, $z\left(x\right)$ is a regionalized random variable satisfying the second-order stationary assumption, $h$ is a step size (i.e., the distance between ordered data), $z\left(x_i\right)$ is the sample value at the spatial point $x_i$, and $N\left(h\right)$ is the total number of sample point pairs at the separation distance $h$.

The most commonly used semivariogram models are spherical, exponential, Gaussian, and pure nugget models [54]. By comparison, the semivariogram function of GWL in the LBP can be described by the following Gaussian model.

The Gaussian model is described using normal probability distribution curves and is suitable for situations where phenomena are at close distances; the Gaussian model is formulated as follows:

$$r\left(h\right)=\left\{\begin{array}{c}0\\ C_0+C\left[1-exp{\left(-h/a\right)}^2\right]\\ C_0+C\end{array}\right. \begin{array}{c}h=0;\\ 0<h\le \sqrt{3}a\\ h>\sqrt{3}a\end{array},$$

(2)

where $C_0$ represents the nugget, which indicates the value of the semivariogram function due to measurement errors and spatial variation when two sampling points are very close together. It reflects random changes in regional variables over small distances, influenced by factors such as local human activities, such as regional groundwater exploitation, and meteorological conditions. $C_0+C$ is the sill, a constant value at which the semivariogram stabilizes as the distance between sampling points increases by $h$.

This parameter captures broader variations caused by structural factors, including geological conditions as well as human influences such as extensive cropland irrigation. $a$ is the range, which refers to the distance between the sampling point when the semivariogram reaches the sill from the nugget, reflecting the range of spatial correlation [53]. $C_0/(C_0+C)$ refers to the spatial variability ratio, reflecting the degree of spatial correlation of GWL in the region. Higher values suggest greater spatial heterogeneity due to random components, while lower values indicate strong spatial autocorrelation of the variable. From the perspective of structural factors, a ratio of less than 25% indicates strong spatial correlation, a ratio between 25% and 75% indicates moderate spatial correlation, and a ratio greater than 75% suggests weak spatial correlation of the variables [55]. The calculation of the semivariogram model is implemented in GS+ 9.0.

2.3.2. Transfer Function Analysis (TFA)

The transfer function analysis describes the relationship between the input and output signals of a system in the frequency domain, aiming to estimate the parameters of GWL response to precipitation. In general, the transfer function decomposes the stationary input and output signals into a sum of sines and cosines at multiple frequencies using the Fourier transform. Assuming linearity, the sine wave at the input will be converted into a sine wave with the same frequency at the output. In this paper, the transfer function analysis method is used to analyze the multi-timescale spatiotemporal characteristics of the response of the shallow GWL to precipitation in the LBP.

In the time domain, the spectral estimation in the frequency domain is obtained by Fast Fourier Transform, and the main period of the time series is extracted. The transfer function analysis method is introduced in the frequency domain, used to determine the coherence, gain, and phase of the GWL response to the precipitation signal, in order to quantify the existence, strength, and lag time, which are the properties of the response system. The conceptual illustration of transfer function analysis method is shown in Figure 6. Before the TFA, the time series should be detrended to remove the climatological annual cycle and linear trend.

The transfer function is specified as follows:

Setting shallow groundwater level as a response signal $X$, with precipitation set to forced signal $Y$, the transfer function between the two signals is $H\left(f\right)$ [30]:

$$H\left(f\right)=S_{xy}\left(f\right)/S_{xx}\left(f\right),$$

(3)

The $S_{xx}\left(f\right)$ is the spectrum of the $X$ signal, $S_{xy}\left(f\right)$ is the cross-spectrum of the two signals. The transfer function magnitude (gain) $\left|H\left(f\right)\right|$ and phase spectrum $\varphi \left(f\right)$ are derived from the real part $H_r\left(f\right)$ and the imaginary part $H_i\left(f\right)$ of the complex transfer function as

$$\left|H\left(f\right)\right|={\left\{{\left[H_r\left(f\right)\right]}^2+{\left[H_i\left(f\right)\right]}^2\right\}}^{{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}},$$

(4)

$$\varphi \left(f\right)=arctan\left[H_i\left(f\right)/H_r\left(f\right)\right],$$

(5)

The existence of the relationship between $X$ and $Y$ is evaluated by the coherence function $\theta \left(f\right)$, which is defined as

$$\theta (f)={\left|S_{xy}\left(f\right)\right|}^2/\left[S_{xx}\left(f\right)\ast S_{xy}\left(f\right)\right],$$

(6)

Under the linear assumption, each term in the forcing series is transformed to have the same frequency as the response term, but a different amplitude (i.e., gain), usually shifted in time.

Using the coherence function aids in identifying unreliable cases of gain and phase estimation. Coherence ranges from 0 to 1 at each frequency, ensuring the consistency and coordination of input and output signals. In a linear system characterized by a high signal-to-noise ratio (SNR) and a single-variable input–output relationship, coherence approaches 1. Conversely, when SNR is low, or when the system exhibits high nonlinearity or is influenced by additional variables, coherence tends toward 0. Thus, a higher coherence (closer to 1) indicates a stronger correlation between GWL and precipitation.

The gain reflects the ability of a response system to amplify or attenuate an input signal at different frequencies. In the frequency domain, gain represents the ratio of the amplitude of the response signal amplitude to the forcing signal, characterizing the efficiency of signal transfer. If the gain remains constant, the system exhibits uniform transfer effects for input signals across all frequencies. Conversely, if gain varies with frequency, the system’s response to signals of different frequencies will differ. A higher gain signifies a more robust response of the response signal to the forcing signal; specifically, when precipitation changes by one unit, the corresponding fluctuation in GWL is more pronounced.

The phase represents the delay effect of the output signal relative to the input signal within a system. In the frequency domain, as signals of varying frequencies traverse the system, the output signal’s phase shifts concerning the input signal due to the system’s internal factors. This shift is the phase. The phase captured by transfer function analysis reflects the response and feedback process between the input and the output signals. However, due to the negligible effect of the groundwater on precipitation in the LBP [56], this study primarily considers how precipitation influences GWL fluctuations. Positive phase values may indicate long-term response processes that current data cannot adequately capture [32]. When combined with frequency information, the phase lag can be converted into time lag, quantifying the temporal lag of GWL in response to precipitation.

Given the local heterogeneity of factors affecting GWL fluctuations, the groundwater response to precipitation exhibits significant temporal variability, ranging from rapid response within a day to delays spanning months and years [17]. To elucidate the characteristics of the precipitation–GWL response across different time scales, four distinct time scales are defined based on frequency ranges: annual, 6-month, 4-month, and less than 4-month time scales. The maximum coherence for each frequency range is sequentially calculated as the coherence parameter for the corresponding time scale. The gain and time lag for frequencies exhibiting significant coherence are averaged within each frequency range, and these averages are designated as the gain and time lag for that time scale. The TFA method was implemented with the “mtspec” package in Python.

2.3.3. Apportionment Entropy (AE) Seasonality Index

Apportionment entropy is an indicator that describes the randomness of data distribution [57]. We introduced it to assess the uniformity of the GWL distribution in the LBP throughout the year, thereby quantifying the seasonal variation in shallow GWL.

In order to evaluate the AE value of GWL change in $k$ years, calculate the sum of the groundwater levels $x_m\left(m=1, 2,\cdots ,12\right)$ in each month of $k$ years firstly, denoted as $X_k$.

$$X_k={\sum }_{m=1}^{12}{x}_{m,k},$$

(7)

Then, the AE value of the year $k$ can be calculated as

$${AE}_k=-{\sum }_{m=1}^{12}{\displaystyle \frac{x_{m,k}}{X_k}}{\mathit{log}}_2\left({\displaystyle \frac{x_{m,k}}{X_k}}\right),$$

(8)

From the formula, when the GWL is evenly distributed in each month of the year (the groundwater level is equal in each month), the AE value reaches the maximum value of ${\log }_{2}12$. In contrast, when GWL peaks in one month throughout the year and is very low in the remaining months, the AE value is close to 0.

The AE values characterize the seasonal variation in GWL. In the study, a high AE means that GWL is evenly distributed during the year with less seasonal variation, while a low AE implies that GWL fluctuates more dramatically with stronger seasonality.

2.3.4. Moran’s Index

The local Moran’s I can determine the spatial autocorrelation of each spatial unit with neighboring spatial units, as well as quantitatively identify the differences in spatial autocorrelation in different regions [58]. At spatial location $i,$ the local Moran’s I_i is calculated as follows:

$$I_i=Z_i{\sum }_{j=1}^n{W}_{ij}{Z}_j,$$

(9)

where $Z_i$ and $Z_j$ are the standardized deviations of the observed values from the mean of units $i$ and $j$, and $W_{ij}$ is the spatial weight matrix. When calculating the model, the spatial weight matrix needs to be row-standardized.

The statistic score for the local Moran’s $I_i$ significance test is $Z\left(I_i\right)$. The formula is as follows:

$$Z\left(I_i\right)={\displaystyle \frac{I_i-E\left(I_i\right)}{\sqrt{VAR\left(I_i\right)}}},$$

(10)

where $I_i$ means the local Moran’s I, $E\left(I_i\right)$ is the hope of $I_i$, and $VAR\left(I_i\right)$ means the variance of $I_i$.

In this study, the local Moran’s I was employed to assess spatial autocorrelation in seasonal fluctuations in GWLs across various monitoring wells in the LBP. When the $Z\left(I_i\right)$ score is significant, if $I_i$ is positive, it indicates that spatial unit $i$ and its neighbors share similar attribute values (both high or low). A positive $Z\left(I_i\right)$ suggests that the local Moran’s $I_i$ for spatial unit $i$ exceeds the expected value, designating it as a high–high agglomeration area. Conversely, a positive $I_i$ paired with a negative $Z\left(I_i\right)$ indicates that the local Moran’s $I_i$ is below the expected value, characterizing it as a low–low agglomeration area.

3. Results

3.1. Spatiotemporal Characteristics of Precipitation and Shallow Groundwater Level

3.1.1. Spatiotemporal Characteristics of Precipitation in the LBP

The statistical analysis of the GPM precipitation dataset revealed significant spatial inhomogeneity in LBP precipitation, exhibiting a discernible pattern. The northeastern coastal region received the highest average annual precipitation, totaling 1047.1 mm. Precipitation decreased from southeast to northwest, reaching a minimum of 577.3 mm in Guanxian County, Liaocheng City. The average annual precipitation in the LBP from 2000 to 2020 was 667.2 mm, with 2002 recorded as the driest year (383.41 mm) and 2020 as the rainiest year (944.44 mm). Overall, we observed a slightly upward trend in average annual precipitation, with a growth rate of 8.06 ± 5.10 mm/yr. Additionally, the precipitation pattern in the LBP also had obvious seasonal differences, with 51.5% of the annual precipitation occurring in July and August, whereas December and January saw significantly less precipitation, contributing only about 1%. This annual precipitation distribution was markedly uneven.

3.1.2. Characteristics of Spatiotemporal Variability of Shallow Groundwater Level

We effectively modeled the semivariogram function of the shallow GWL in the LBP using the Gaussian model, as illustrated in Figure 7. The parameters derived from the semivariogram function facilitated further analysis of the spatial and temporal variations in the autocorrelation of the shallow GWL. In this study, the nugget $C_0$, reflects the influence of small-scale human activities, such as regional groundwater exploitation, on the spatial distribution of GWL, and the sill $C_0+C$, captures the effects of structural factors, including geological conditions and topography, as well as anthropogenic influences like extensive crop irrigation and large-scale water diversion projects.

As shown in Figure 8, the nugget $C_0$ fluctuates around 0.5, while the sill $C_0+C$ is significantly higher, approximately 2.5. The ratio of the two, $C_0/\left(C_0+C\right)$, remains stable below 25%, indicating strong spatial autocorrelation of the GWL in the LBP, where change patterns in adjacent areas exhibit notable similarity. In addition to geological influences, large-scale crop irrigation contributes to a consistent change trend in GWL over extensive areas, further enhancing spatial autocorrelation. From 2000 to 2004, $a$ showed a significant downward trend, decreasing from 406 km to 350 km in five years. $C_0$ decreased from 0.26 to 0.16 by 2004 (Figure 8), suggesting a reduction in the impact of regional variables and an increase in the spatial autocorrelation of GWL. The GWL in the LBP showed greater stability, particularly in the center of the groundwater depression cone along the southeastern Boxing–Guangrao County, which stabilized at −15 ± 3 m. The minimal influence of small-scale human activities resulted in high spatial autocorrelation of GWL. From 2008 to 2018, $C_0+C$ fluctuates steadily around 2.4 ± 0.1, remaining higher than $C_0$, indicating that the changes in GWL during this period were predominantly driven by large-scale structural factors. Throughout this time, the groundwater depression cone in the Boxing–Guangrao area fluctuated around −20 m. In 2018, the $C_0$ value shifted from a gradual decline to a rapid increase, rising from 0.19 to 0.26 within a year, suggesting that small-scale human activities, such as regional groundwater exploitation, became the primary influence during this period. Figure 9e–f show that the groundwater level along the western Ningjin region stabilized within the range of 2 ± 2 m. Meanwhile, the GWL in the southeastern Boxing–Guangrao area increased from −20.22 m in 2018 to −15.8 m in 2020. In contrast, continuous exploitation of shallow groundwater in Yangxin County led to a decline from 4.25 m to 2.15 m (Figure 9b); the GWL in Yangxin County dropped from 4.25 m to 2.15 m, causing an imbalance in GWL changes across the LBP and diminishing spatial autocorrelation. By this time, the entire LBP region had evolved into a spatial distribution characterized by multiple groundwater depression cones in Ningjin, Boxing–Guangrao, and Wudi–Yangxin.

3.2. Spatiotemporal Characteristics of Precipitation–GWL Response

3.2.1. Spatial Variability of the Precipitation–Groundwater Level Response in the LBP

On the basis of the analysis of spatial and temporal variability characteristics of precipitation and shallow GWL, in this study, we further explored the response relationship between them in the LBP region from 2000 to 2020, utilizing the TFA method. The spatial distribution is shown in Figure 10. The coherence spatial distribution indicates significant coherence between GWL and precipitation, with the highest coherence value of 0.96 observed in the Zouping area south of the LBP, decreasing to 0.76 in the Ningjin area in the northwest direction. This trend aligned with the spatial variation in annual precipitation in the LBP (Figure 3a). As regional precipitation increased, the volume of water that effectively infiltrated into the aquifer also rose. Consequently, in areas with heavy precipitation, GWL responded significantly to changes in precipitation, resulting in a higher coherence between the two values. In contrast, we observed anomalously low coherence in Dezhou City and Guangrao County, Dongying City, with coherence values of 0.62 and 0.60, respectively. We attributed this result to both regions being classified as groundwater-overexploited areas (Figure 2c). Long-term groundwater exploitation has disrupted the recharge process of precipitation to groundwater, thereby weakening the correlation between GWL and precipitation in these areas.

In the groundwater-overexploited areas, we identified two typical zones, A and B, which exhibited spatial differences in precipitation–GWL response characteristics. In zone A, the lag time of GWL response to precipitation was longer in area a, averaging 5.3 months, whereas in area b, the lag time was significantly shorter, approximately 2.1 months (Figure 10b). The gain values of GWL response to precipitation within these areas displayed similar spatial distribution characteristics, with a lower gain value of approximately 0.63 in area a and a higher value of about 1.65 in area b (Figure 10c and Figure 11a,b). This indicates that GWL fluctuates more significantly in area b during unit precipitation events, reflecting a higher effective recharge rate of precipitation to groundwater [11]. The typical zone A is located in the alluvial plain (Figure 2b), where the primary lithology of the aquifer consists of fine sand and silt, with medium sand present along the Guanxian–Ningjin area. This aquifer has a high sandy soil content (Figure 2d). Due to the low water-holding capacity and high permeability of sandy soil, precipitation can effectively recharge the aquifer. Consequently, the lag time of GWL fluctuations in response to precipitation is shorter, and there are more pronounced fluctuations. Thus, the gain value of GWL response to precipitation in area b is high, with a short lag time. In typical area B, the lag time of GWL response to precipitation in area c shows a low-value cluster of about 2.7 months, while the lag time in area d near Guangrao County is longer, reaching up to 9.7 months. Conversely, the gain value of GWL response to precipitation in area c (0.9) is slightly higher than in area d (0.52) (Figure 11c,d). Typical area B is situated at the junction of the Piedmont alluvial plain, alluvial plain, and alluvial–marine plain, with area c located in the Piedmont alluvial plain and alluvial plain. The high permeability of the aquifer results in a shorter lag time for GWL to respond to precipitation. In contrast, area d, located in the alluvial–marine plain, is primarily composed of silty clay soil (Figure 2c). The low permeability of the clay renders the water level in the aquifer less sensitive to changes in precipitation, leading to a low gain value for GWL response to precipitation and a longer lag time. Lorenzo-Lacruz’s research [12] supports the findings of this study, indicating that clay soil has a significant impact on the precipitation–GWL response process.

In the long-term time-series analysis, precipitation served as the primary source of groundwater recharge, with GWL fluctuations being primarily influenced by climate change. We observed a significant correlation between shallow GWL fluctuations and precipitation across the entire region. The response characteristics exhibited regional differences, primarily influenced by climatic factors, aquifer characteristics, and human activities [18]. In regions with substantial precipitation, the correlation between GWL and precipitation was stronger. Conversely, in areas with intense groundwater exploitation, the relationship between GWL and precipitation was disrupted, resulting in a weaker correlation. Additionally, aquifer characteristics largely determine the strength and lag time of the GWL’s response to precipitation.

3.2.2. Multi-Timescale Response of Groundwater Level to Precipitation

The response of GWL to precipitation is a complex hydrological process influenced not only by regional geological conditions but also by human activities, such as groundwater exploitation. In conjunction with the TFA results, we further analyzed this multi-timescale response across several intervals: annual, 6 months, 4 months, and less than 4 months.

On the annual timescale, the coherence between GWL and precipitation was generally weak in most regions, with a mean coherence value of 0.66, except for the southern Zouping region, which had a value of 0.94. On the 6-month timescale, the correlation improved, reaching approximately 0.73 (Figure 12a,d). Regarding the gain value and lag time of GWL response to precipitation, the annual mean gain value was relatively low at 0.57. In contrast, on the 6-month timescale, the gain value was significantly higher in the west (approximately 1.2) but remained lower in the east (Figure 12b,e), showing strong spatial variability consistent with hydrogeologic unit zoning (Figure 2b). As discussed in Section 3.2.1, the western alluvial plain of the LBP features an aquifer primarily composed of fine sand and silt, exhibiting high permeability. Consequently, the recharge rate of precipitation to groundwater was higher, as reflected in the elevated gain value of groundwater response. Additionally, the high permeability resulted in a shorter lag time of about 3.5 months. Conversely, in the eastern alluvial–marine plains, the low permeability of clay limited the infiltration of precipitation into the aquifer, leading to a weaker GWL response and a longer lag time of approximately 4.2 months (Figure 12c,f).

On the 4-month and less than 4-month timescales, the response of GWL to precipitation showed more significant coherence, with mean values reaching up to 0.78 (Figure 12g,j). Except for the main groundwater-overexploited areas of Dezhou City and Dongying City, coherence in other regions predominantly exceeded 0.8. Correspondingly, the gain value of GWL response precipitation significantly improved compared with the annual scale (0.57) and the 6-month scale (0.92) (Figure 12b,e,h,f). On the less than 4-month timescale, the gain value peaked at 3.49 in the Lingcheng County of Dezhou City. Thus, this response was likely influenced by the seasonal characteristics of precipitation on this timescale (Figure 3b), demonstrating a more pronounced correlation between groundwater response and precipitation. This finding aligned with the results of the study by Shamsudduha et al. [59]. From a spatial perspective, the distribution characteristics of gain values were generally consistent with those observed on the annual timescale. The difference in coherence between GWL and precipitation, however, was more pronounced in overexploited groundwater areas compared with areas that were not overexploited. Coherence reached 0.96 in the Zouping area but only 0.55 in the southern Guangrao (Figure 12j). Additionally, the lag time for GWL response to precipitation in areas that were not overexploited was approximately 0.6 months, whereas it was about 1.17 months in groundwater-overexploited areas (Figure 12l), nearly doubling the response time. In groundwater-overexploited areas, the impact of exploitation on GWL was immediate and pronounced. Therefore, on shorter timescales (e.g., less than 4 months), the coherence effect of groundwater exploitation on precipitation–GWL response was more evident, thus highlighting spatial differences in coherence. Furthermore, groundwater exploitation could lead to a dropdown of the GWL, with an increasing unsaturated zone thickness, resulting in longer infiltration times for precipitation. This extended response time manifested as a prolonged lag time for GWL to react to precipitation changes. [29,60,61]

On shorter timescales, however, the response of shallow GWL to precipitation exhibited significant immediate characteristics, indicating a more rapid response. For instance, during a brief period of heavy precipitation, GWL may rise sharply within 5 days [62]. Nonetheless, due to the limitations imposed by the time resolution of the time series and the Nyquist sampling theorem, this study did not analyze the immediate response of GWL to precipitation.

In keeping with the results noted in Section 3.2.1, we found that GWL fluctuations were influenced primarily by precipitation fluctuations on longer timescales (i.e., interannual scales). The characteristics of the response of GWL to precipitation were mainly affected by factors such as the spatial distribution of precipitation and aquifer characteristics. Conversely, on shorter timescales (e.g., 4 months or less), this response relationship was more readily influenced by human activities, such as groundwater exploitation. This conclusion aligned with findings from previous studies [32,63].

4. Discussion

4.1. Seasonal Fluctuations in the Shallow Groundwater Level

The shallow GWL in the LBP exhibited significant seasonal characteristics, as shown in Figure 13. From January to March, the GWL remained relatively stable, fluctuating around 9.1 m. Starting in April, the GWL gradually decreased, reaching its lowest point of 8.4 m in July. As precipitation increased in July, however, the groundwater was recharged and rebounded, reaching a peak of 9.2 m in September. As a result of regional crop irrigation practices, the winter wheat irrigation period extends from November to April of the following year. During this time, precipitation was minimal, and irrigation relied primarily on groundwater and diversions from the Yellow River. Consequently, the GWL was influenced by exploitation and the diversion of water, resulting in a tendency for this level to stabilize or decline. Research by Yan et al. [64] confirmed the impact of groundwater exploitation and cropland irrigation on GWL fluctuations. During the non-irrigation season (May to October), fluctuations in GWL are influenced primarily by precipitation recharge, leading to a tendency for recovery. Therefore, this intrayear fluctuation in GWL reflected the effects of irrigation-induced groundwater exploitation and seasonal precipitation recharge.

To further analyze the seasonal fluctuations in GWL in the LBP, we introduce the apportionment entropy seasonality index [57] to quantify the uniformity of GWL distribution over a year. Figure 14a illustrates the spatial autocorrelation results of GWL fluctuations in the LBP. Red indicates areas with high AE values, signifying weak seasonal fluctuations, and blue indicates areas with low AE values, signifying strong seasonal fluctuations. Strong seasonal fluctuations were concentrated in the marine plain in the northeast of the LBP, whereas areas with weak seasonal fluctuations were found primarily in Dezhou City, and Boxing–Guangrao in the east of the LBP. As discussed in Section 3.2., groundwater exploitation in the Xiajin and Boxing–Guangrao areas has weakened the aquifer’s response capacity to precipitation recharge, resulting in GWLs that do not fully recover after summer precipitation. Consequently, the amplitude of seasonal GWL fluctuations diminishes, leading to weaker seasonal fluctuations. Therefore, high AE values cluster in the Dezhou and Boxing–Guangrao areas, indicating that the seasonal GWL fluctuations in these regions were weak.

4.2. Seasonal Characteristics of Precipitation–GWL Response Affected by Irrigation Exploitation

Six irrigation points were selected based on land use types, as shown in Figure 14b. Considering the characteristics of regional crop irrigation, we classified the period from November to April of the following year as the irrigation period and classified May to October as the non-irrigation period. The precipitation–GWL response characteristics of the six irrigation points during these periods are illustrated in Figure 15.

All of the six irrigation points exhibit similar characteristics, with the coherence of the GWL response to precipitation being higher during the non-irrigation period than during the irrigation period. The coherence value at Well-07 during the non-irrigation period reached 0.8, which was significantly higher than the 0.3 observed during the irrigation period on the same timescale (Figure 15c). Additionally, the gain value of the GWL response to precipitation during the non-irrigation period was markedly higher than during the irrigation period, with this difference becoming more pronounced as the timescale decreased. For instance, the gain value at Well-06 remained around 0.5 during the irrigation period, whereas during the non-irrigation period, it increased from 0.5 to 3.5 in the 4-month timescale (Figure 15b). This result indicated that GWL had a significant response to precipitation during the non-irrigation period, allowing the aquifer to receive more effective precipitation. Since winter wheat primarily relies on groundwater for irrigation during the irrigation period [63], irrigation water also penetrates downward, recharging groundwater [65]. Consequently, during the irrigation period, the GWL was influenced by exploitation and other factors, reducing the dominance of precipitation, which resulted in a low coherence between GWL and precipitation. In contrast, during the non-irrigation period, GWL fluctuations were influenced predominantly by precipitation recharge. Heavy summer rainfall effectively recharged groundwater in a short time [24]. Therefore, coherence between GWL and precipitation was higher in the non-irrigation season, reflecting a stronger strength of the response of GWL to precipitation. In the NCP, 70–80% of groundwater is used for irrigation [63]. Since the 1970s, the Yellow River has been utilized for irrigation in the LBP. Yellow River water effectively replenishes groundwater in the irrigation area, with infiltration recharge from Yellow River irrigation accounting for approximately 20–40% of total shallow groundwater resources [48]. Unlike the water-saving methods that restrict irrigation in other areas of the NCP [66], the Yellow River irrigation method in the LBP has alleviated the pressure placed on farmland irrigation caused by water shortages [48,67] and, to some extent, has replenished groundwater, improved the regional water cycle, and promoted the recovery of groundwater resources.

4.3. Limitations and Constraints

In this study, we combined the GPM precipitation dataset and long-term GWL dynamic monitoring data and used the TFA method to explore the multi-timescale characteristics of the response of GWL to precipitation changes, combined with human activities, such as groundwater exploitation and crop irrigation. The limitations of this study include the lack of GWL data in the study area, the complexity of the water-cycle process, and hydrogeological conditions.

First, the lack of GWL data and the uneven distribution of groundwater monitoring wells made it impossible to fully and accurately reflect the groundwater dynamics in the study area. Therefore, areas lacking monitoring wells may have had deviations or blind spots in the characteristics of the response of GWL to precipitation changes. Second, the hydrogeological conditions of the LBP were complex, and the distribution and characteristics of multiple aquifers and impermeable layers affected the regional groundwater flow and water-cycle process. At the same time, strong and complex human activities (e.g., groundwater exploitation, agricultural irrigation, inter-basin water transfer) were intertwined and together had an impact on the regional groundwater recharge–discharge process, further increasing the complexity of the precipitation–GWL response system. These complexities increased the uncertainty of this study’s findings.

In this study, we effectively captured the multi-timescale response of GWL to precipitation, providing a new perspective for understanding the response of groundwater systems. Given the limitations of this study, however, in the future, we will try to collect GWL data from a wider area to supplement deficiencies in existing data. Additionally, in view of the water-cycle process and complex aquifer characteristics, we will consider combining multiple models, such as machine learning and hydrological models, to provide a more comprehensive and accurate understanding of the physical process of the GWL response to precipitation.

5. Conclusions

The response of shallow groundwater to precipitation in the LBP exhibits new characteristics under the dual influences of climate change and human activities, such as groundwater exploitation. A new research framework is established to analyze the multi-timescale characteristics of the response of shallow groundwater to precipitation under the influence of human activities and applied to the LBP. The main conclusions are as follows:

(1) The distribution of precipitation and shallow GWL in the LBP shows significant spatial variability. The average annual precipitation in the LBP is 667.2 mm, decreasing spatially from southeast to northwest. Precipitation is unevenly distributed throughout the year, with July and August accounting for 51.5% of total annual precipitation. The semivariogram of shallow GWL can be well fitted with the Gaussian model, with the $C_0/(C_0+C)$ value stable below 25%, indicating high spatial autocorrelation. This spatial variability is influenced by large-scale structural factors, such as geological conditions and extensive cropland irrigation.

(2) Shallow GWLs generally exhibit a significant response to precipitation in the LBP. The coherence between GWL and precipitation reaches a maximum value of 0.96 in the Zouping areas with high annual precipitation. However, due to the long-term effects of groundwater exploitation, the coherence in Dezhou City and Dongying City is relatively low (approximately 0.6). Additionally, regional differences in the response of GWL to precipitation are influenced by hydrogeological conditions and aquifer characteristics. Areas with higher sandy soil content in the aquifer show greater effective recharge from precipitation, with a gain value of approximately 1.65. Correspondingly, the lag time in GWL response to precipitation is shorter (2.1 months).

(3) GWL fluctuations are primarily influenced by precipitation on the longer timescale (interannual). The characteristics of the response of GWL to precipitation are mainly affected by factors such as the spatial distribution of precipitation and aquifer characteristics. Conversely, on shorter timescales, the spatial characteristics of the precipitation–groundwater response are more constrained by human activities, such as groundwater exploitation. On timescales of less than 4 months, the lag time for GWL response in overexploitation areas (0.98 months) is nearly 60% longer than in areas without overexploitation (0.6 months).

(4) Groundwater exploitation weakens the seasonality of regional groundwater fluctuations. During the irrigation period (November to April), the GWL is influenced by various factors, including groundwater exploitation and irrigation infiltration, resulting in low coherence with precipitation (the lowest being 0.5). In contrast, during the non-irrigation period (May to October), GWL fluctuations are predominantly driven by precipitation recharge. The correlation between GWL and precipitation is high (0.8), and the response to precipitation is strong, with a gain value reaching up to 3.5.

This study focuses on the precipitation–GWL response process, incorporating the impacts of human activities, such as groundwater exploitation and cropland irrigation, to provide a deeper understanding of the characteristics of the response of shallow GWL to precipitation changes and their spatial and temporal variability. Those findings are significant for the optimization and rational allocation of water resources in the LBP. The framework provides a new perspective for other regions to study groundwater response issues and a theoretical basis for the sustainable development of water resources in both the long and short term.