Variation of Groundwater Flow Caused by Any Spatiotemporally Varied Recharge

: The objective of this study was to develop a complete analytical solution to determining the e ﬀ ect of any varying rainfall recharge rates on groundwater ﬂow in an unconﬁned sloping aquifer. The domain of the unconﬁned aquifer was assumed to be semi-inﬁnite with an impervious bottom base, and the initial water level was parallel to the impervious bottom of a mild slope. In the past, similar problems have been discussed mostly by considering a uniform or temporally varying recharge rate, but the current study explored the variation of groundwater ﬂow under temporally and spatially distributed recharge rates. The presented analytical solution was veriﬁed by comparing its results with those of previous research, and the practicability of the analytical solution was validated using the 2012 and 2013 data of a groundwater station in Dali District of Taichung City, Taiwan.


Introduction
Groundwater has been closely linked with human life since the beginning of civilization. Agriculture, industry, and people's livelihoods inevitably require groundwater resources. This study focused on the distribution of groundwater during replenishment. Obtaining accurate information of groundwater changes and distribution in an extensive aquifer is difficult, even with modern technical equipment and mathematical models.
Studies on groundwater flow in unconfined aquifers have mostly used the Boussinesq equation as the governing equation based on the Dupuit-Forchheimer approximation. Most researchers have derived the Boussinesq equation for a sloping aquifer by incorporating Darcy's law and the continuity equation (Bansal and Das [1]; Chapman [2]; Childs [3]; Wooding and Chapman [4]). Glover [5] and Hantush [6] were the first ones to analyze the evolution and stability of a free surface with uniform replenishment in a semi-infinite aquifer. Since then, many mathematical models have been developed to predict the dynamic behavior of groundwater in response to constant or periodic recharge in unconfined aquifers (Manglik et al. [7]; Rai and Manglik [8]).
Verhoest and Troch [9] and Pauwels et al. [10] obtained analytical solutions of a linearized Boussinesq equation by considering a uniform rainfall recharge and used the solutions to estimate changes in groundwater level and outflow discharge. Sanford et al. [11] used another linearization approach to derive analytical expressions for water table profiles and verified their solution with laboratory experiments. Su [12] derived an analytical solution of the Boussinesq equation by linearizing and transforming the equation into an ordinary heat conduction equation and then used this solution to estimate the temporally varying recharge of the groundwater table. Ramana et al. [13] presented analytical solutions of the linearized Boussinesq equation to predict the development of a phreatic surface by approximating the transient recharge as an exponentially decaying function of time.
Recently, Rahman et al. [14] developed a free surface boundary condition called Groundwater Flow Boundary (GFB) to discuss groundwater dynamics in a more computationally efficient way than a three-dimensional model, and evaluated the impact of the assumptions on model results. Sappa et al. [15] presented a rainfall-flow cross-correlation model to evaluate the spring flow patterns of the main Capodacqua di Spigno Spring in the study area. The results obtained from the model were compared with existing methods using the Standard Precipitation Index (SPI) to estimate the minimum annual spring discharge.
Bansal and Das [1] divided their study domain into a recharge area and a non-recharge area to simulate groundwater level changes in a semi-infinite aquifer. They used the Laplace transform method to develop an analytical expression of the water head distribution in an aquifer because of uniform infiltration from a single recharge basin. The present study used Darcy's law as a theoretical basis and applied the continuity equation to govern the groundwater flow; the study considered the effect of the inclination of an unconfined aquifer without separating it into two zones.

Mathematical Formulation
The schematic of an unconfined aquifer with a semi-infinite domain overlying a sloping impermeable bed is shown in Figure 1. The aquifer's recharge varies with space and time; is the range of recharge in space, and h 0 is the initial water depth. According to Bansal and Das [1], Darcy's law results in the following equation: where q is the flow rate in the x direction per unit width of the aquifer (L 2 /T), k is the hydraulic conductivity (L/T), H w is the elevation of the groundwater table measured perpendicular to the underlying impermeable bottom (L), and θ is the aquifer inclination angle.  [14] developed a free surface boundary condition called Groundwater Flow Boundary (GFB) to discuss groundwater dynamics in a more computationally efficient way than a three-dimensional model, and evaluated the impact of the assumptions on model results. Sappa et al. [15] presented a rainfall-flow cross-correlation model to evaluate the spring flow patterns of the main Capodacqua di Spigno Spring in the study area. The results obtained from the model were compared with existing methods using the Standard Precipitation Index (SPI) to estimate the minimum annual spring discharge. Bansal and Das [1] divided their study domain into a recharge area and a non-recharge area to simulate groundwater level changes in a semi-infinite aquifer. They used the Laplace transform method to develop an analytical expression of the water head distribution in an aquifer because of uniform infiltration from a single recharge basin. The present study used Darcy's law as a theoretical basis and applied the continuity equation to govern the groundwater flow; the study considered the effect of the inclination of an unconfined aquifer without separating it into two zones.

Mathematical Formulation
The schematic of an unconfined aquifer with a semi-infinite domain overlying a sloping impermeable bed is shown in Figure 1. The aquifer's recharge varies with space and time; ℓ is the range of recharge in space, and ℎ 0 is the initial water depth. According to Bansal and Das [1], Darcy's law results in the following equation: where is the flow rate in the x direction per unit width of the aquifer (L 2 /T), is the hydraulic conductivity (L/T), is the elevation of the groundwater table measured perpendicular to the underlying impermeable bottom (L), and is the aquifer inclination angle. Considering a net rainfall recharge into the ground, a continuity equation based on the law of conservation of mass can be expressed as follows: Considering a net rainfall recharge into the ground, a continuity equation based on the law of conservation of mass can be expressed as follows: To investigate groundwater flow, Equation (1) can be combined with Equation (2) and the Dupuit-Forchheimer approximation can be invoked, yielding the following Boussinesq equation for a sloping unconfined aquifer: where S is the specific yield of the aquifer (-) and r(x, t) is the spatiotemporally varying recharge (L/T). Equation (3) is a nonlinear partial differential equation and thus does not produce general solutions. According to Brutsaert [16], the nonlinear term in Equation (3) can be linearized by replacing the first H w with εD, where D is the thickness of the initially saturated aquifer and ε is a parameter given by 0 < ε < 1. The Boussinesq Equation (3) can be linearized to the following equation: with the following initial and boundary conditions: In practice, a hyetograph is usually used to record and display the rainfall distribution. Kazezyılmaz-Alhan [17] introduced the Heaviside function, u(−), to describe the hyetograph, denoting r = r(t) only. In this study, the rainfall distribution could be expressed as follows: in which r ij (L/T) (i = 1, 2, 3, . . . , n and j = 1, 2, 3, . . . , m) indicates the rainfall amount over the corresponding time and space interval. This technique was used in the work of Wu and Hsieh [18], but in this study we appended the spatial variation of recharge and estimated the recharge from rainfall. In Equation (8), when j is 1, the whole domain (0 < X < ∞) is divided into two areas-recharge and non-recharge areas, as indicated in Figure 1, but the solution can be obtained once without incorporating two solutions for the two sub-domains.

Analytical Solution
To nondimensionalize the governing equation, initial condition, and boundary conditions, the following relations can be introduced: where t d is the duration of recharge. Notably, H is a dimensionless variable denoting the increase in the groundwater level after recharge. Moreover, all the other variables on the left-hand side of the  (8) and (9) into Equations (4)-(7) yields G.E. : To comply with the generalized integral transform method (GITM), we must initially eliminate the first-order spatial derivative term by using the following transformation: Substituting Equation (14) into Equations (10)-(13) results in We define the integral transform and inversion formula for H v (X, T) in Equation (15) with respect to the space variable X in the semi-infinite region 0 ≤ X < ∞ as follows: Integral transform Inversion formula where ξ(λ, X) is the kernel function. For the present problem, both boundary conditions are of the first kind, and, according to Özisik [19], the kernel function ξ(λ, X) can be expressed as follows: where λ assumes all values from zero to infinity continuously. Taking the integral transform to Equation (15) yields the following equation: The first integral on the left-hand side of Equation (22) can be evaluated using integration by parts twice, consequently resulting in the following equation: Water 2020, 12, 287

of 17
After Equations (17) and (18) are substituted into the first term on the right-hand side of Equation (23), the term is cancelled, and the result reads as follows: Therefore, Equation (22), according to Equation (19), becomes in which, The solution of Equation (25) can be determined as follows: If the inversion formula, Equation (20), is applied to Equation (27), the solution of H v can be expressed as follows: Then, substituting Equation (28) into Equation (14) yields Equation (29) indicates that the groundwater level fluctuates with the spatiotemporally varying recharge rate. Furthermore, the actual groundwater level can be expressed as follows: Without loss of generality, the recharge in the space interval X = 0 to 1 is illustrated so that the upper and lower limits of the space integral in Equation (29) are reduced to X = 0 to 1; therefore, with the substitution of dimensional variables, we obtain the following equation: The presence of the term tanθ in Equation (32) implies that a constant amount of water always tends to flow in the aquifer because of the bottom slope. The first term with the negative sign is due to water mounding in the aquifer, which creates impedance in the groundwater flow along the semi-infinite extent. Consequently, a net inflow or outflow exists depending on the recharge rate and the slope parameter. Notably, the present solutions, Equations (31) and (32), are superior to that in previous research, for example, the work by Bansal and Das [1], because spatiotemporally varied parameters can be handled in this solution.

Results and Discussion
Bansal and Das [1] used the Laplace transform method to solve the linearized Boussinesq equation and explored the effect of the constant recharge rate into a semi-infinite region on the groundwater table. The present analytical solution was verified with the results of Bansal and Das [1]. The following parameters were used: ε = 1/3, recharge length = 100 m, soil layer thickness D = 7 m, hydraulic conductivity k = 0.001 m/s, soil water storage constant S = 0.34, initial water depth h 0 = 2.5 m, and the recharge rate r = 4 mm/h. The results are illustrated in Figure 2. As shown by the three curves, the peak water level and its location were determined to be related to the inclination angle. The aquifer with a higher slope had a lower peak water level because of a higher groundwater discharge rate. The peak water level in the steep sloping aquifer was close to the end of recharge extent, X = 1. Because of the large inclination, the peak water level was reached rapidly. Figure 2 shows that the results of this study are consistent with the results of Bansal and Das [1], especially for the cases of steeper slopes.
On the right-hand side of Equation (10), the second-order partial differential term and the first-order partial differential term are the hydraulic diffusion term and the extended term, respectively. In addition to the head difference, the main factors controlling the groundwater flow are hydraulic conductivity k and aquifer inclination angle θ. Therefore, we investigated the effects of the permeability factor α and speed factor U on the groundwater response. As presented in Figure 3, the peak was determined to occur when the recharge was completed, and the water level decreased considerably as the permeability factor increased. The decaying time at which the groundwater level returned to the initial water level changed with various permeability factors. Specifically, the dimensionless decaying time at which the water levels separately decayed to half the maximum water levels was approximately T = 10 for α = 0.0397, T = 4 for α = 0.0794, T = 3.3 for α = 0.2382, T = 2.2 for α = 0.3971, and T = 2 for α = 0.7941, signifying that the flow through the aquifer increased with the permeability, thus reducing the water level. However, as shown in Figure 4, the water levels decreased as the flow velocity factors increased. The decaying time at which the water levels separately decayed to half the maximum water levels was approximately T = 4 for U = 0, T = 2.7 for U = 0.0772, T = 2.3 for U = 0.1556, T = 1.6 for U = 0.3212, and T = 1.4 for U = 0.5094, indicating that the flow increased with the slope, thus reducing the water level. To summarize, the water level variation is more sensitive to α than to U; that is, the influence of α on water level changes is more prominent than that of U. On the right-hand side of Equation (10), the second-order partial differential term and the firstorder partial differential term are the hydraulic diffusion term and the extended term, respectively. In addition to the head difference, the main factors controlling the groundwater flow are hydraulic conductivity k and aquifer inclination angle . Therefore, we investigated the effects of the permeability factor and speed factor U on the groundwater response. As presented in Figure 3, the peak was determined to occur when the recharge was completed, and the water level decreased considerably as the permeability factor increased. The decaying time at which the groundwater level returned to the initial water level changed with various permeability factors. Specifically, the dimensionless decaying time at which the water levels separately decayed to half the maximum water levels was approximately T = 10 for = 0.0397, T = 4 for = 0.0794, T = 3.3 for = 0.2382, T = 2.2 for = 0.3971, and T = 2 for = 0.7941, signifying that the flow through the aquifer increased with the permeability, thus reducing the water level. However, as shown in Figure 4, the water levels decreased as the flow velocity factors increased. The decaying time at which the water levels separately decayed to half the maximum water levels was approximately T = 4 for U = 0, T = 2.7 for U = 0.0772, T = 2.3 for U = 0.1556, T = 1.6 for U = 0.3212, and T = 1.4 for U = 0.5094, indicating that the flow increased with the slope, thus reducing the water level. To summarize, the water level variation is more sensitive to than to U; that is, the influence of on water level changes is more prominent than that of U.    We designed a unimodal rainfall event with a duration of 11 h as a source of recharge to understand the effects of the speed factor U and the permeability factor on the groundwater translation. As illustrated in Figure 5, when the speed factor U was not considered (i.e., U = 0), the groundwater translation was only controlled by the permeability factor , therefore, the downstream translation of the groundwater level peak was slow. When the speed factor increased, as shown in Figures 6 and 7, the downstream translation of the groundwater became more pronounced and rapid. For a fixed location (X = constant), the groundwater level considerably varied with the speed factor. This observation is corroborated from Figures 6 and 7; when the speed factor was doubled, the peak We designed a unimodal rainfall event with a duration of 11 h as a source of recharge to understand the effects of the speed factor U and the permeability factor α on the groundwater translation. As illustrated in Figure 5, when the speed factor U was not considered (i.e., U = 0), the groundwater translation was only controlled by the permeability factor α, therefore, the downstream translation of the groundwater level peak was slow. When the speed factor increased, as shown in Figures 6 and 7, the downstream translation of the groundwater became more pronounced and rapid. For a fixed location (X = constant), the groundwater level considerably varied with the speed factor. This observation is corroborated from Figures 6 and 7; when the speed factor was doubled, the peak water level swiftly moved downstream. Specifically, the speed factor had a more significant influence on the groundwater translation than did the permeability factor. However, when the permeability factor α was the only factor controlling the transmission of groundwater, the downstream flow rate of the groundwater was slow, and the groundwater easily accumulated in the aquifer within the recharge area. By contrast, when we considered the speed factor U, we observed a fairly significant downstream movement of groundwater, and the groundwater was easily drained from the aquifer in the recharge area. These phenomena were graphically similar to those of water level; therefore, the corresponding figures are not presented herein. Moreover, as indicated in Figures 5-7, the groundwater level increased rapidly within the duration, and the level peaked at the end of the duration and was less affected by both the speed factor and the permeability factor. downstream movement of groundwater, and the groundwater was easily drained from the aquifer in the recharge area. These phenomena were graphically similar to those of water level; therefore, the corresponding figures are not presented herein. Moreover, as indicated in Figures 5-7, the groundwater level increased rapidly within the duration, and the level peaked at the end of the duration and was less affected by both the speed factor and the permeability factor.  To validate the presented solution, we collected observational data regarding the groundwater level and rainfall distribution at the gauge stations in the Wu River watershed in Taichung, Taiwan. Figure 8 illustrates the locations of the rainfall stations and groundwater wells within an 8-km range. In the Wu River basin, the soil in the areas from Taiping Mountain to Taichung is mostly Ultisols and Alfisols, which are categorized as red soil. The average hydraulic conductivity (k) of Ultisols and Alfisols measured by Obi and Ebo [20] in a laboratory experiment was 3 m/day.   To validate the presented solution, we collected observational data regarding the groundwater level and rainfall distribution at the gauge stations in the Wu River watershed in Taichung, Taiwan. Figure 8 illustrates the locations of the rainfall stations and groundwater wells within an 8-km range. In the Wu River basin, the soil in the areas from Taiping Mountain to Taichung is mostly Ultisols and Alfisols, which are categorized as red soil. The average hydraulic conductivity (k) of Ultisols and Alfisols measured by Obi and Ebo [20] in a laboratory experiment was 3 m/day. In reality, water table changes do not fully reflect rainfall changes, but there is a considerable correlation between the two. Clark et.al. [21] estimated the magnitude of surface infiltration from the change of the relationship between groundwater storage and outflow. Based on the local geography data of this study, we used 25% of the rainfall intensity as the recharge rate by referring to the suggestion of Clark et al. [21]. The rainfall data and local soil type are based on the public information of the Taiwan government, and then according to the indoor experiments of Obi and Ebo [20], the typical values of the corresponding parameters k and S can be obtained. After incorporating the relative parameters and recharge data into the derived solution, we observed that the simulated water level was relatively consistent with the groundwater level data, as shown in Figure 9. Figure  9a depicts a heavy recharge event with a recharge intensity of 50 mm/day at the Taichung Dali station on 10 June 2012, which resulted in an approximately 1-m increase in the groundwater level. Later, a 45-day non-recharge period caused a slight decrease of approximately 0.5 m in the groundwater table. Subsequently, a typhoon brought 125 mm/day rainfall, increasing the groundwater level by approximately 2 m. The variation between the present solution and the measured groundwater level is acceptable. Figure 9b depicts the measured groundwater level data at the Taichung Dali Station in 2013. The rainfall in 2013 was more evenly distributed, and the average rainfall intensity was lower than that in 2012. To quantitatively assess the difference between analytical and observational data, we performed an error analysis by evaluating the relative percentage difference (RPD), which is defined as follows: From the results, it is known that the error is less than 5% throughout the simulation duration, as shown in Figure 9c. The simulated curves of this analytical solution are close to the observed data; thus, the results are satisfactory. In reality, water table changes do not fully reflect rainfall changes, but there is a considerable correlation between the two. Clark et.al. [21] estimated the magnitude of surface infiltration from the change of the relationship between groundwater storage and outflow. Based on the local geography data of this study, we used 25% of the rainfall intensity as the recharge rate by referring to the suggestion of Clark et al. [21]. The rainfall data and local soil type are based on the public information of the Taiwan government, and then according to the indoor experiments of Obi and Ebo [20], the typical values of the corresponding parameters k and S can be obtained. After incorporating the relative parameters and recharge data into the derived solution, we observed that the simulated water level was relatively consistent with the groundwater level data, as shown in Figure 9. Figure 9a depicts a heavy recharge event with a recharge intensity of 50 mm/day at the Taichung Dali station on 10 June 2012, which resulted in an approximately 1-m increase in the groundwater level. Later, a 45-day non-recharge period caused a slight decrease of approximately 0.5 m in the groundwater table. Subsequently, a typhoon brought 125 mm/day rainfall, increasing the groundwater level by approximately 2 m. The variation between the present solution and the measured groundwater level is acceptable. Figure 9b depicts the measured groundwater level data at the Taichung Dali Station in 2013. The rainfall in 2013 was more evenly distributed, and the average rainfall intensity was lower than that in 2012. To quantitatively assess the difference between analytical and observational data, we performed an error analysis by evaluating the relative percentage difference (RPD), which is defined as follows:

RPD =
Hw ana − Hw obs Hw obs (33)   The variation of the flow discharge at different locations, X = 0.5, 1, and 1.5, under given rainfall events is shown in Figure 10. The outflow increased with rainfall, and further downstream, more flow discharge occurred because of the accumulation effect. It also can be seen from Figure 10 that in the case of low recharge rates, the difference in flow rate at different locations of the aquifer is not significant. From the results, it is known that the error is less than 5% throughout the simulation duration, as shown in Figure 9c. The simulated curves of this analytical solution are close to the observed data; thus, the results are satisfactory.
The variation of the flow discharge at different locations, X = 0.5, 1, and 1.5, under given rainfall events is shown in Figure 10. The outflow increased with rainfall, and further downstream, more flow discharge occurred because of the accumulation effect. It also can be seen from Figure 10 that in the case of low recharge rates, the difference in flow rate at different locations of the aquifer is not significant. observational data.
The variation of the flow discharge at different locations, X = 0.5, 1, and 1.5, under given rainfall events is shown in Figure 10. The outflow increased with rainfall, and further downstream, more flow discharge occurred because of the accumulation effect. It also can be seen from Figure 10 that in the case of low recharge rates, the difference in flow rate at different locations of the aquifer is not significant.

Conclusions
Because the observation of the groundwater level of an aquifer is not as easy as that of the surface water level, this study proposed an analytical solution instead of a numerical model for simulating or predicting the groundwater level under a tempo-spatially varying recharge rate. When the presented analytical approach is used, research funding can be considerably reduced and feasible results can be rapidly obtained. Some remarkable points are summarized as follows: (1) The present results obtained by solving the Boussinesq equation using the GITM are consistent with the results of Bansal and Das [1] obtained using the Laplace transform; this is because both analytic solutions are close to each other. However, the presented analytical approach can handle a tempo-spatially varying recharge rate, whereas the approach in the previous study can handle only a constant recharge rate. (2) In a constant bottom slope, has a considerable influence on water level changes, whereas U

Conclusions
Because the observation of the groundwater level of an aquifer is not as easy as that of the surface water level, this study proposed an analytical solution instead of a numerical model for simulating or predicting the groundwater level under a tempo-spatially varying recharge rate. When the presented analytical approach is used, research funding can be considerably reduced and feasible results can be rapidly obtained. Some remarkable points are summarized as follows: (1) The present results obtained by solving the Boussinesq equation using the GITM are consistent with the results of Bansal and Das [1] obtained using the Laplace transform; this is because both analytic solutions are close to each other. However, the presented analytical approach can handle a tempo-spatially varying recharge rate, whereas the approach in the previous study can handle only a constant recharge rate. (2) In a constant bottom slope, α has a considerable influence on water level changes, whereas U is the major influencing factor on groundwater translation for a fixed permeability. (3) To validate the practicability of the present solution, a simulation was performed under given rainfall events, and the results were satisfactory when compared with the observed data. Therefore, the presented approach can be applied to predict the groundwater level and flow discharge, and thus it would be beneficial to groundwater management.