A Modified Water-Table Fluctuation Method to Characterize Regional Groundwater Discharge

Abstract

A modified Water-Table Fluctuation (WTF) method is developed to quantitatively characterize the regional groundwater discharge patterns in stressed aquifers caused by intensive agricultural pumping. Two new parameters are defined to express the secondary information in the observed data. One is infiltration efficiency and the other is discharge modulus (recurring head loss due to aquifer discharge). An optimization procedure is involved to estimate these parameters, based on continuous groundwater head measurements and precipitation records. Using the defined parameters and precipitation time series, water level changes are calculated for individual wells with fidelity. The estimated parameters are then used to further address the characterization of infiltration and to better quantify the discharge at the regional scale. The advantage of this method is that it considers recharge and discharge simultaneously, whereas the general WTF methods mostly focus on recharge. In the case study, the infiltration efficiency reveals that the infiltration is regionally controlled by the intrinsic characteristics of the aquifer, and locally distorted by engineered hydraulic structures that alter surface water-groundwater interactions. The seasonality of groundwater discharge is characterized by the monthly discharge modulus. These results from individual wells are clustered into groups that are consistent with the local land use pattern and cropping structures.

1. Introduction

Human activities in many regions have greatly altered groundwater systems [1,2]. In areas where groundwater is greatly exploited, such as areas with intensive agricultural pumping, many previously connected aquifer systems have evolved into multiple isolated hydrogeological units. As a result, vertical water exchanges including infiltration and pumping have become the dominant flow patterns in those highly stressed unconfined aquifers [3,4]. In such systems, groundwater recharge and discharge are critical components for understanding fluid and contaminant transport [5]. Actual rates of groundwater recharge and discharge, however, are difficult to estimate.

Recharge to unconfined aquifers is derived from the downward infiltration from the land surface of precipitation, irrigation water, and/or seepage from surface water bodies. It is a complicated process influenced by land use, land cover, vegetation type, and soil properties [5,6]. In a natural system, recharge to the shallow aquifer is largely dominated by precipitation characteristics; primarily duration, magnitude, and intensity of precipitation [7,8,9]. Whereas, in an agriculturally stressed aquifer, the irrigation return flow can contribute a large portion of total recharge [10]. This part of the recharge is governed by the pumping and irrigating related to agricultural practices [11]. Moreover, in such systems, the quantity of the discharge is significant and varies spatially and temporally based on crop types, planting schedule, and climatic factors. Since the precipitation process and the pumping and irrigating processes follow intrinsically different periods, deriving the underlying mechanisms based on the water level signals is a complicated task [10,12]. Therefore, for stressed aquifers where groundwater flow is predominately vertical, new paradigms are needed to understand and characterize such aquifers to develop integrated water resource management plans.

Scientists have long been using water-table hydrographs to estimate key aquifer parameters. The Water-Table Fluctuation (WTF) method was used as early as the 1920s [13] and since then has been applied and modified in numerous studies [14,15,16] to estimate recharge from the infiltration of precipitation. Theoretically, when the system meets certain assumptions, the WTF method can also be used to estimate evapotranspiration from fluctuations in groundwater levels [17,18].

The WTF method has several important variations, including the RISE method, the MRC (Master Recession Curve) method, and the EMR (Episodic Master Recession) method [19]. The distinctions among these variations are whether they consider the water level records in terms of a continuous process or hydrologic episodes, data requirement for the resolution, and whether the method contains correction for water level recession caused by different mechanisms.

The RISE Method is the simplest variation of WTF for the estimation of recharge. The continuous water level records are divided into equal time intervals, and the water level rise for a given interval is the amount the water table changed from the previous interval. A decline is considered as a rise of zero. The RISE method is simple and objective, but this method does not correct for the recessional processes and therefore tends to underestimate the recharge [20]. The RISE method requires high frequency water level records.

The MRC method does not identify individual recharge episodes. It was developed as an automated or semi-automated procedure to determine the rate of water-table decline (dH/dt) as a function of H, and to predict the slope of the hydrograph as a function of the current level in the absence of recharge [20,21]. Once the MRC is established, positive deviations of the hydrograph from the MRC are attributed to recharge. Like RISE, the MRC method requires continuously monitored water level data. It avoids the subjectivity through the estimation of the MRC parameters, but does not correct for hydrograph recession.

The EMR method [19] utilizes high resolution data records to identify discrete recharge episodes based on WTF rates and quantify the recharge for individual episodes, and uniquely associates each episode with a causal water-input period. A recharge episode is defined as a period during which the total recharge rate significantly exceeds its steady-state condition in response to a substantial recharge event. The EMR method eliminates the recession effects by adding two parameters to MRC, with the advantage of the high resolution data records. The EMR method, facilitated with a water level model developed by Halford [22], was applied to study the relationships between storm characteristics and episodic groundwater recharge [23].

In summary, there are two major limitations in the application of the general WTF methods. First, the WTF methods are mostly designed for natural systems and are not accounting for water level fluctuations from pumping discharge at the regional scale. Second, all WTF methods require the specific yield ($S_y$) of the aquifer to be properly determined beforehand. However, determining the specific yield for WTF analysis is a difficult task, even with carefully planned and executed aquifer tests [24].

The objective of this study is to develop a new method named as Water-Table Fluctuation Regression (WTFR) in response to the two issues with the general WTF methods mentioned above. It conceptualizes the pumping influences into recurring monthly steps, with recession corrections. As the basis of WTFR method, the long-term water-table hydrograph of a single well is defined by two parameters: infiltration efficiency $\beta$ and discharge modulus $DM$. The former is a dimensionless ratio between precipitation amount in length and groundwater head increment; and the latter is a 12-element array that depicts the seasonal pattern of net discharge. Both parameters are objective derivations from a water-table hydrograph, and for many wells in stressed aquifers these two parameters can reproduce the entire water-table hydrograph with fidelity. Thus, it is not necessary to estimate the specific yield. The authors of this paper developed a user interface to perform WTFR. A beta-version of this interface (WTFR 1.4) and related documentations are provided in the supplementary materials.

2. Study Area and Data Sets

We selected the Dagu Aquifer located in the east coastal area of China as the study area due to its highly stressed condition with groundwater overexploitation, as well as the availability of sufficient and high-resolution water table data. In addition, the Dagu Aquifer is a relative isolated hydrogeologic unit, which is facilitative for the method testing. The Dagu Aquifer (Figure 1) is located on the alluvial plain of the Dagu River watershed and encompasses 430 square kilometers. The Dagu Aquifer extends from north to south along the ancient Dagu River valley. Bedrock outcrops on the east and west are the boundaries of the alluvial terraces. The modern river channel overlies ancient alluvial materials that fill the valley and comprises the majority of the Dagu Aquifer. The bottom of the valley is Cretaceous sandstone and shale bedrock. The thickness of the aquifer ranges from 5 m to 20 m, and the depth of water table ranges from 2 m to 7 m, with an average of 5 m [25].

Under natural conditions, the Dagu River drains toward the south and eventually empties into the Jiaozhou Bay. However, a non-permeable barrier was constructed at the southern boundary of the aquifer in 1998 to prevent saltwater intrusion. This structure greatly reduced the quantity of the groundwater discharging to the sea and slowed down groundwater circulation within the aquifer [26]. Moreover, this shallow sandy unconfined aquifer provides municipal water supply and irrigation water for the thriving local agriculture industry near the Dagu River. The groundwater withdrawals for agricultural and municipal purposes have become the largest hydrologic stress on the aquifer.

The long-term dynamics of the local groundwater resources are dictated by climate cycles and human exploitation. Based on reports from the local Hydrologic Bureau, the annual net agricultural and municipal withdrawals are approximately 6.2 × 10⁷ m³, which is equivalent to 31% of the total aquifer storage. The recharge of the aquifer is dominated by rainfall infiltration and irrigation return flow, accounting for 58% and 20% of the total recharge, respectively. Given the limited upstream runoff from both surface water and groundwater, precipitation and pumping have become the dominant process of local groundwater cycling. Because of its shallow sandy properties, the Dagu Aquifer quickly responds to precipitation, as well as anthropogenic derived pumping and irrigating activities.

The average annual precipitation of 625 mm is distributed unevenly, both spatially and temporally. Precipitation decreases from the south to north and it has significant seasonal variations, with a wet season in the summer from July to September where 60% to 70% of annual precipitation occurs. The average potential water surface evapotranspiration is 984 mm, which is concentrated in the summer season between April and September. The temperature ranges from lows of −20 °C in January to highs of 37.4 °C in August, with a yearly average temperature of 12.5 °C.

As an important agricultural economic area and a source of municipal water supply, the Dagu aquifer has 115 monitoring wells instrumented to collect data for water resource management. Based on the data quality, 82 of these were selected for WTFR method application. Precipitation data are recorded at three meteorological stations located in the Dagu Aquifer, along with other meteorological parameters. As shown on Figure 1, the 82 wells are distributed relatively evenly within the aquifer. They were divided into three precipitation zones based on their distances from the three meteorological stations. The well names starting with A, B, and C represent the three meteorological zones, respectively.

3. Methodology

3.1. The Basis of WTFR Method

The general WTF methods assume that the water level fluctuation is directly proportional to the precipitation flux entering the aquifer, as expressed in Equation (1):

$$S_y\times \Delta H=P\times \alpha$$

(1)

where $S_y$ is the specific yield [Dimensionless], $\Delta H$ is the groundwater level rise, $P$ is the accumulated precipitation [L] during the period of water level rises, and $\alpha$ is the infiltration rate defined as the portion of water applied at the ground surface that reaches the saturated aquifer [Dimensionless].

The WTFR method is designed to characterize systems that are driven by both precipitation recharge and net discharge processes. The discharge processes are conceptualized into two components, as shown in Equation (2):

$$S_y\times \Delta H=P\times \alpha -ND-RV$$

(2)

where $ND$ is the net groundwater pumping discharge volume per unit area, expressed as water column height [L]; and $RV$ is the groundwater recession volume per unit area, expressed as water column height [L]. Net discharge represents the combined effect of pumping and irrigation return. $RV$ is defined as head decrease due to horizontal flow imbalances as governed by Darcy’s Law. Groundwater recession is driven by various factors. Depending on the system flow conditions, the recessional processes might be negligible or determined using case-specific information.

By rearranging Equation (2), the governing equation for the WTFR method is expressed as:

$$\Delta H_i=P_i\times \beta -D{M}_i-R_i$$

(3)

where $\Delta H_i$ is the groundwater head change during the ith day [L]; $P_i$ is precipitation during the ith day [L]; $\beta$ is the infiltration efficiency, defined as the ratio of head increase divided by corresponding precipitation height [Dimensionless]; $D{M}_i$ is the discharge modulus (defined as head decrease due to net aquifer discharge) during the ith day [L]; and $R_i$ is groundwater level recession during the ith day [L].

The parameter $\beta$ regulates precipitation; $D{M}_i$ reflects the net discharge; and $R_i$ reflects the aquifer recessional process. They are defined as:

$$\beta =\frac{\alpha }{S_y}$$

(4)

$$D{M}_i=\frac{N{D}_i}{S_y}$$

(5)

$$R_i=\frac{R{V}_i}{S_y}$$

(6)

Therefore, for a given monitoring well, a water-table hydrograph is reconstructed using the initial value of water level $H_0$, and the sequential head changes $\Delta H_i$, as shown in the following equation:

$$H_{i+1}=H_i+\Delta H_i$$

(7)

Given water-table data and precipitation records that span n days, Equation (3) is populated into n equations. Here, the daily accumulated precipitation records ($P_i$) are known data series, and groundwater level recession $R_i$ is estimated case by case, whereas, $\beta$ and $D{M}_i$ are parameters to be estimated. It is theoretically possible to have multiple combinations of $\beta$ and $D{M}_i$ to generate the same $\Delta H_i$. Therefore, extra constraints are developed for the proper solving of $\beta$ and $D{M}_i$.

3.2. Periodic Water Balance Periods

Figure 2 shows a typical water-table hydrograph. Several instances are marked in the figure where water levels rose as a result of recharge and subsequently declined to the same levels due to discharge and groundwater recession. According to the Equation (3), for any such periods, $\Delta H=0$, therefore:

$$P_p=\frac{D{M}_p+R_p}{\beta }$$

(8)

where $P_p$ is cumulative precipitation during the pth period [L], $D{M}_p$ is cumulative DM during the pth period [L], and $R_p$ is cumulative recession during the pth period [L]. For any water-table hydrograph characterized with periodic water level fluctuations, a number of such periods can be located to constrain the process of solving $\beta$ and $D{M}_i$.

With identified m periodic water balanced periods, m equations are generated based on Equation (3). Unlike the previously defined n equations, these equations are not based on head changes, but rather on periodic lumped recharge and discharge quantities, respectively. This new constraint stabilizes the parameter estimation process and greatly decreases the uncertainty of the WTFR method. Based on daily water level change equation (Equation (3)) and lumped periodic water balances (Equation (8)), we have set up the water-table hydrograph simulation model. The model parameters are $\beta$ and $D{M}_i$, and the model outputs are calculated $\Delta H_i$ in daily interval and the periodically accumulated precipitation $P_p$. Correspondingly, n observed $\Delta H_i$ and m observed $P_p$ from a given monitoring well and nearby precipitation records are employed to set the objective function for the model. The approach used to estimate parameter $\beta$ and $D{M}_i$ of the WTFR method is described in the following section.

3.3. Parameterization of DM Values

It is not practical to solve directly for $\beta$ and $D{M}_i$ (i = 1, 2, 3, …, n) as there are n + 1 unknowns. Additional assumptions are established to reduce the number of parameters. First, since regional pumping and irrigating is most likely related to seasonal agricultural activities, it is logical to reduce the daily $DM$ value into monthly steps. Additionally, when we perform a multi-year calculation, assigning a $DM$ value for every single month would still result in many arbitrary parameters. Therefore, we assume the seasonality of the discharge is consistent, i.e., the monthly discharge quantities have the same pattern each year. This assumption is reasonable as the probabilities of extreme situations, which largely alter the discharge patterns, are much lower than the general conditions in a given hydrologic system. As a result, we defined a DM(12) array, which contains 12 elements and each of these elements represents a recurring monthly $DM$ value. These $DM$ values in turn are used to populate constant daily $D{M}_i$ values for the entire month. Experiences dictate that this approach greatly strengthens the seasonal pattern and reduces the uncertainty introduced by defining and solving DM values for each individual month. However, during our application of WTFR, we discovered that sometimes reducing DM into a 12-element array over-constrains the fitting process. This occurs, particularly, when dealing with long-term data series, as the recurring 12-month pattern tends to shift in response to extreme weather conditions such as drought or storm, the development of engineered hydraulic structures, and alteration to a new farming pattern, etc. Under such circumstances, additional parameters can be introduced to provide relaxation to the fitting process.

After reducing the numerous daily DM values to a recurring monthly DM(12) array, the n + m equations are redefined and used to estimate the 13 parameters.

3.4. Solving for β and DM Array

Solving an overdetermined system with the ordinary least squares method is a standard approach that has been implemented in many existing codes. PEST, a model-independent parameter optimizer, developed by Doherty [27], was used as the solver in this study. A user interface was developed to run the PEST executable and perform the WTFR calculation. This interface breaks down a multi-well data set into individual cases, prepares these cases into input files recognizable by the PEST solver, creates a necessary run-time environment for PEST calculations, and post-processes optimization results into user-friendly spreadsheet formats. PEST uses a sub-class of the non-linear least squares method known as the Marquardt-Levenberg method.

3.5. Recession Corrections

Groundwater cycling vertically and water levels responding mainly to precipitation and pumping are two of the fundamental assumptions of the WTFR method. When these assumptions are not valid, the recession term R in Equation (3) can be activated to provide add-on packages to WTFR. Several prevalent conceptual models are discussed hereafter.

3.5.1. Natural Recession $R_N$

Natural recession is the groundwater level decline due to natural aquifer discharge. Consider the conceptual model illustrated in Figure 3. This hydrogeological unit is bounded by impermeable boundaries to the north, west, and south; and a discharge boundary to the east. This unit receives precipitation infiltration from a total area of F [L²], and discharges groundwater through the eastern boundary at a rate of Q’ [L³/T]. During the dry season, Q’ is predominantly sustained by the drainage of the aquifer. When ignoring all other sources or sinks, we have:

$$M=\frac{Q^{\prime }}{F}$$

(9)

where M is the groundwater discharge modulus, defined as the groundwater discharge rate per unit area aquifer [L/T].

For the entire hydrogeological unit, the average recession during a given time period of $\Delta t$ can be expressed as:

$$R_N=\frac{Q^{\prime }}{F\times S_y}\Delta t=\frac{M}{S_y}\Delta t$$

(10)

Similarly, if we consider only a strip within the hydrogeological unit bounded by two neighboring groundwater contours (shadowed area in Figure 3), a similar relationship can be derived:

$$R_N=\frac{Q_{dn}-Q_{up}}{A\times S_y}\Delta t=\frac{M}{S_y}\Delta t$$

(11)

where A is the area of the strip [L²]; and $Q_{up}$ and $Q_{dn}$ are the upstream and downstream groundwater discharge rates, respectively.

3.5.2. Dynamic Recession $R_D$

Dynamic recession is the groundwater level decline due to local disturbance of hydrogeological conditions. This type of groundwater recession was addressed by Simon and Rorabaugh [28]. The Rorabaugh solution was expanded to simulate bank storage of a one-dimensional aquifer [29]. Here, we take this basis of the Rorabaugh solution to consider a conceptual model (Figure 4) where an infinite, straight constant head boundary induces dynamic recession. Consider a homogenous aquifer bounded on one side by a linear constant head boundary (Figure 5). Well k is close to the boundary and Well j is far enough away not to be influenced by the boundary. Initially, the water level of the aquifer equals the head at the boundary. Suppose a precipitation event uniformly and instantaneously increases the water level in the aquifer by $\Delta H_j$. The head increase in the aquifer near the constant head boundary is illustrated by the following equations. Since Well j is too far away from the boundary to receive any of its influence, the dynamic recession asserted to Well k by the boundary is expressed by:

$$R_D=\Delta H_j-\Delta H_k$$

(12)

Suppose t = 0 when the aquifer water level rises from $H_b-\Delta H_j$ to $H_b$ during the initial precipitation event with the following initial and boundary conditions:

$$H\left(x,t=0\right)=H_b \forall x$$

(13)

$$H\left(x=0,t\right)=H_b-\Delta H_j \forall t>0$$

(14)

$$\underset{x\to \infty }{\lim }H\left(x,t\right)=H_b \forall t$$

(15)

For simplicity, suppose the aquifer has a uniform thickness of $m$ [L]. The one-dimensional groundwater flow equation is:

$$S_y\frac{\partial H}{\partial t}-Km\frac{{\partial }^{2}H}{\partial x^2}=0$$

(16)

where $S_y$ is the specific yield [Dimensionless]; H is the groundwater head [L]; K is the hydraulic conductivity [L/T]; and x is the perpendicular distance from the boundary [L].

The analytical solution to Equation (16) under definite conditions of Equations (13)–(15) is expressed as:

$$H\left(x,t\right)=\left(H_b-\Delta H_j \right)+\left[H_b-\left(H_b-\Delta H_j\right)\right]erf\left(\frac{1}{2}\sqrt{\frac{x^2{s}_y}{tmK}}\right)$$

(17)

For Well k where x = $x_k$ and $H=H_b-\Delta H_k$, we have:

$$\frac{\Delta H_j-\Delta H_k}{\Delta H_j}=erf\left(\frac{1}{2}\sqrt{\frac{x^2{s}_y}{tmK}}\right)$$

(18)

By substituting Equation (12) into Equation (18), we have:

$$R_{Dt}=erf\left(\frac{1}{2}\sqrt{\frac{x^2{s}_y}{tmK}}\right)\times \Delta H_j$$

(19)

where R_Dt is the dynamic recession at location x induced by the constant head boundary between t = 0 and t = t.

3.5.3. Leakance Recession $R_L$

The leakance recession accounts for the groundwater level fluctuations due to vertical leakance. As shown in Figure 6, the target (upper) aquifer may be leaking to the lower aquifer, or receiving leakance from the lower aquifer. According to Darcy’s law:

$$Q_{\mathrm{leak}}=\frac{A{k}^{\prime }}{m^{\prime }}\left(H_i-H_L\right)=\sigma A\left(H_i-H_L\right)$$

(20)

where $Q_{\mathrm{leak}}$ is the leakance rate [L³/T]; A is the horizontal area of leakance [L²]; $k^{\prime }$ is the vertical hydraulic conductivity [L/T] of the low permeable layer; $m^{\prime }$ is its thickness [L]; $H_i$ is the water head in the upper aquifer at moment i [L]; $H_L$ is the constant water head in lower aquifer [L]; and $\sigma =\frac{k^{\prime }}{m^{\prime }}$ is defined as the leakance rate [1/T]. Therefore, the leakance recession $R_{Li}$ during period i is expressed as:

$$R_{Li}=\frac{Q_{leak}}{A{S}_y}\times \Delta t=\frac{\sigma }{S_y}\left(H_i-H_L\right)\times \Delta t$$

(21)

4. Results

4.1. Data Processing and Parameterization

Short term water-table fluctuations can be caused by changes in temperature or atmospheric pressure [24], air trapped in the medium [30,31], and temporary horizontal flow interferences, etc. These processes are not related to the conceptual model of WTFR, and have no effect on the intrinsic aquifer property $\beta$ and long-term pumping strength indicator DM. Therefore, the original water–table hydrographs were first normalized into daily data series (from 1 August 2010 to 31 July 2012, 730 days in total) and subsequently processed with a five-day moving average procedure to remove daily spikes while retaining seasonal trends.

The WTFR method was applied to the individual well’s water-table hydrograph and the corresponding daily precipitation time series. Each pair of data generate a set of optimized parameters including $\beta$ and DM(12), the calculated water-table hydrograph, and the periodic accumulated precipitation. As an example, the detailed simulation results for one individual well are shown in Figure 7. Figure 7a shows the correspondence between the modeled and observed hydrographs, indicating that the majority of the groundwater level fluctuations are explained by the parameters defined by WTFR. Similarly, Figure 7b shows the good agreement between the modeled and observed periodic cumulative precipitation ($P_p$ in Equation (8)) during the periods bounded by the identified equal water-level pairs. Figure 7c compares the residues from two different sources: water level modeling and cumulative precipitation modeling. For this well, the objective function is dominated by water level residuals. Note that the weights assigned to the two residual sources can be adjusted should the modeler decided to emphasize either one of the sources.

In general, WTFR was effective in reproducing the observed water levels recorded at all 82 wells. The Root Mean Square (RMS) of the water level residual at all wells ranges from 0.05 m to 0.49 m, with an average of 0.22 m. The agreement between the observed and calculated water level implies that the groundwater level fluctuations in this aquifer are controlled by a linear system that allows the simple superposition of water levels. The good agreement between observed and modeled water levels indicates that precipitation and regional pumping combined explain the majority of the underlying dynamics of groundwater level fluctuations. Although precipitation is not a true linear process, it approximates one on the temporal and spatial scale of this exercise. The overall goodness of optimization is represented by the final objective function $\mathrm{\Phi }$, which is a measure for the model-to-observation discrepancies as described in the PEST [27]. A lower $\mathrm{\Phi }$ value represents better agreement between the model and real system. The averaged $\mathrm{\Phi }$ value for the 82 wells was 44.4, calculated for a range from 2.2 to 175.6.

The robustness of the optimization results was tested by applying different combinations of initial values of $\beta$ and DM(12). The optimization results were generally consistent, which further demonstrates the linear characteristics of the system. Sensitivity analysis shows that $\beta$ is the most sensitive parameter in the simulation. It is generally more sensitive than the individual DM(12) parameters by one or more orders of magnitudes. This is consistent with the fact that $\beta$ is a global parameter, while each DM(12) only controls one twelfth of the simulated period.

4.2. Seasonality of Discharge

The calculated monthly DM values vary from 0.1 mm to 238.3 mm in 82 wells, with the average of 25.1 mm and the standard deviation of 26.0 mm. This result implies that both the temporal and spatial variability of pumping is relatively high in the study area. Furthermore, we used the cluster analysis method, an un-supervised classification approach where similar data are placed into related or homogeneous categories without advanced knowledge of the categories’ definitions [32], to categorize the large number of DM(12) values into a smaller, manageable number of groundwater discharge categories. Each of these categories consists of similar discharge patterns (i.e., pumping schedules), which can be used as a quantitative indicator for further study. In this study, a combination of the Ward’s minimum variance method, the Euclidean distance measure, and internal normalization of DM(12) vectors produced the most logical clustering of discharge patterns within the Dagu Aquifer. As the result, the 82 sets of the calculated DM(12) array were grouped into three categories. The monthly DM values for each well were plotted as the histograms in Figure 8b–d according to the clustering category. Additionally, the averaged monthly DM values for each category are plotted in the figures as red boxes.

4.3. Infiltration Efficiency

The calculated β (infiltration efficiency) values for the 82 wells range from 2.3 to 9.7, with an average of 5.2. These $\beta$ values were spatially interpolated within the extent of the Dagu Aquifer using the Ordinary Kriging method, and the contours calculated with a linear variogram model are shown in Figure 9.

5. Discussions

Recharge and discharge parameters are essential inputs for regional groundwater flow models. However, the distributed values are difficult to acquire and calibrate. In the absence of more complicated analytical methods or costly field measurements, WTFR analysis can serve as a first attempt to understand the dynamics of the local groundwater cycle.

5.1. Timing of Hydrologic Signals

WTFR characterizes the recharge-discharge rivalry between precipitation and regional pumping. Accurate regressions of these two processes demand near instantaneous responses and drastically different temporal patterns. Note that the discharge modulus is defined as the net head reduction due to the collective result of pumping discharge and irrigation return. In some cases, the pumping signals and the irrigation return signals do not arrive at the monitoring well simultaneously. As a result, when these signals superpose at the monitoring well, phase shifts can be included in the final DM signal. Such discrepancies tend to reduce the quality of the signal, and in extreme cases, demand new conceptual models.

5.2. Physical Significance of $\beta$

As shown in Figure 9, the distribution of infiltration efficiency is marked by sporadic depressions along the Dagu River channel. These depressions correlate with the locations of the existing dams on the river channel. Since $\beta$ reflects the water table sensitivity to rainfall, the low values of $\beta$ near these dams indicate that the ponded water dampens the rise and fall of the water table. In addition to the sporadic depressions associated with dams, there are regional depressions away from the river channel, notably in the center of the northern section, and in the far southern sections of the aquifer. The locations of these regional depressions correspond with areas of thicker saturated aquifer materials. This correspondence is not fully understood and warrants further comparative studies.

Infiltration efficiency $\beta$ is considered to be constant during the WTFR calculating period. Although this assumption may not hold in every temporal or spatial range, the application of WTFR in the case study area suggests that $\beta$ being constant is a valid assumption for a multi-year hydrograph. Note that this assumption becomes invalid if there have been drastic water level changes, for instance, relating to man-made structures, as such head change would not only affect the infiltration process, but also potentially change the specific yield because the actual values of specific yield may vary vertically while the water table rises or falls at a large scale [1,24].

Infiltration efficiency $\beta$ is a lumped concept that is linked to the infiltration rate and specific yield (Equation (4)). The case study in this paper shows that the distribution of infiltration efficiency is regionally controlled by the intrinsic characteristics of the aquifer, and locally distorted by surface water. Insights into the control factors of infiltration efficiency need to be gathered from study sites with sufficient surface water-groundwater monitoring information and a robust estimate of specific yield. Better understanding of this parameter will contribute to the characterization of the surface water-groundwater interaction mechanisms. Note that the localized infiltration efficiency depressions near the dams on the river channel represent distortions by a nearby head boundary. Until such distortions are explicitly expressed as recessions in the WTFR governing equation, these outliers should not be used to calculate infiltration rates using Equation (4).

5.3. Spatio-Temporal Aspects of the Aquifer Discharges

The locations of each category correspond geographically with the land use types, as shown in Figure 8a. In the Dagu Aquifer, the agricultural fields are located mainly in the northwest, which corresponds to cluster 1; the green house areas are concentrated in the northeastern and middle sections, which corresponds to cluster 2; whereas the southern part is more urbanized, which corresponds to cluster 3.

The agricultural fields are generally planted with summer vegetable—winter wheat rotations. Accordingly, the discharge pattern in these areas is marked by relatively consistent water consumption throughout the growing season from April to October (Figure 8b). The winter wheat demands intensive irrigation from April to early June; and the summer vegetable crop demands irrigation from the middle of June to October. The winter wheat is planted in October and the dormancy period ends in the next March.

The areas with concentrated greenhouses are identified based on the local satellite image, as shown in Figure 8a. This DM category has a different irrigation schedule from the agricultural field category, as shown in Figure 8c. In the Dagu River region, the greenhouses are operated all year round, with the exception of between late July and late August due to the high temperature and the need to maintain the soil fertility.

The southern portion of the aquifer is largely urbanized and groundwater pumping is driven by municipal water supply. There are some municipal water supply wells located on both sides of the Dagu River downstream reach. As shown in the third clustered category (Figure 8d), the greatly increased DM in July reflects the intensive groundwater pumping that occurs in July due to the lack of water supply from surface water resources in this month. This is consistent with the local river gauge data during the studied years.

5.4. Uncertainty and Limitation

Since the WTFR method is designed for the stressed aquifer, the following assumptions are embodied in the method:

• Groundwater cycles vertically, and horizontal flow has negligible or quantifiable contributions to water-table fluctuation;
• Groundwater table fluctuations are dominated by the combined effect of rainfall infiltration and regional groundwater pumping;
• The delay time of groundwater table responding to such stimuli is on the order of several days or less.

Accordingly, the WTFR method has limitations and uncertainties regarding recession corrections, response time, and evapotranspiration, etc. As defined by the governing equations of WTFR, besides the precipitation signal as regulated by infiltration efficiency, all remaining signals are lumped into discharge modulus. For well-studied sites, some recession signals may be retrieved from the discharge modulus, and explicitly expressed in the governing equation. Such modifications are site specific and should be individually evaluated. Defining the recession items involves subjective judgment, and requires careful estimation and calibration to reduce the uncertainty.

The distance between the wells and precipitation gauges is another potential source of uncertainty in the WTFR. We mitigated this uncertainty by grouping the wells to the closest meteorological stations. Evapotranspiration, which is not considered in the WTFR method, is potentially an important factor for water-table fluctuation if the aquifer is shallow and the evapotranspiration process is intensive. If that is the case, the regional pumping process can be overestimated because the evapotranspiration contributions are lumped into the discharge modulus.

As a relatively simple and effective aquifer management instrument, WTFR has potential applications on multiple fronts: (1) generating recharge distributions for flow models; (2) locating hydrogeological anomalies; (3) understanding surface water-groundwater interactions; (4) locating hidden hydrogeological boundaries such as in an urban setting; (5) calibrating regional agricultural discharge; and (6) serving as a feedback mechanism for groundwater model assimilation.

6. Conclusions

A Water-Table Fluctuation Regression (WTFR) method was developed to simulate the groundwater head changes and characterize the regional discharge patterns in stressed aquifers. The assumptions of the WTFR method are: (1) groundwater mainly cycles in the vertical direction; (2) water table fluctuations are caused by the combined effect of rainfall infiltration and regional pumping; and (3) the delay time in the water-table responding to such stimuli is small. The two parameters in WTFR are $\beta$ (infiltration efficiency, the water table’s sensitivity to recharge) and DM (discharge modulus, the recurring head loss due to aquifer discharge). WTFR includes a fully automatable procedure to estimate the two parameters based on daily groundwater head measurements and the corresponding precipitation records. Using the estimated parameters and precipitation time series, water-table hydrographs from the case study area were reproduced with fidelity. The average RMS for the modeled and measured water levels is 0.22 m. The spatial distribution of the parameter $\beta$ is supported by real-world configuration of the aquifer, and the DM patterns are supported by the existing local land use type.