Analytical Study of Unsteady Nested Groundwater Flow Systems

Two analytical solutions using segregation variable method to calculate the hydraulic head under steady and unsteady flow conditions based on Tóth’s classical model were developed.The impacts of anisotropy ratio, hydraulic conductivity (K), and specific yield (μ s ) on the flow patterns were analyzed. It was found that the area of the equal velocity region increases and the penetrating depth of the flow system decreases at steady state with anisotropy ratio increases, which is defined as ε = √K x /K z . In addition, stagnant zones can be found in the flow field where the streamlines have opposite directions.These stagnant zones move toward the surface as the horizontal hydraulic conductivity increases. The results of the study on transient flow indicate that a relative increase in hydraulic conductivity produces a positive impact on hydraulic head and a relative enhancement in specific yield produces a negative effect on hydraulic head at early times.


Introduction
Tóth [1,2] analytically calculated patterns of gravity-driven groundwater flow systems for small drainage basins using the hydraulic potential theory of Hubbert [3].It was found that only one large regional groundwater flow system could be developed for a homogeneous and isotropic basin with a linearly sloping water table.However, local, intermediate, and regional groundwater flow systems could develop with a water table attributable to the topographic relief of a sinusoidally undulating water table superimposed on a linearly sloping straight line.The nested groundwater flow systems are essential for understanding a variety of geologic processes, such as fluid flow in sedimentary basins [4], surface water-groundwater interaction [5], regional heat and solute transport [6], and the role of groundwater in the study of ore deposition and petroleum accumulation [7,8].Tóth [9] presented an extensive overview of the history, principles, methods, applications, and natural effects of gravity-driven groundwater flow system.
A progression of groundwater flow system investigations has been conducted based on the assumption that the medium is isotropic and homogeneous [2,25,26] or heterogeneous [14,27,28].For example, Tóth [2] investigated the impact of different depth models on the groundwater flow systems; and Jiang et al. [14] investigated the impact of hydraulic conductivity on groundwater flow systems by using the COMSOL software.
Based on a critical review of the literature, one can see that most of the previous studies focused on the steadystate flow of groundwater.However, the groundwater flow is usually transient in the fields.Models of steady-state groundwater flow can result in significant limitations in applications.Besides, the prevailing analytical solutions of unsteady groundwater flow are essentially targeted on determining well flow issues.There are few specific reports about analytical solutions for unsteady nested groundwater flow systems [29].In order to have a better understanding about the evolution process of groundwater flow systems from unsteady to steady under natural and mining conditions, it is necessary to carry out the theoretical investigation on the unsteady groundwater flow systems.This study focuses on the analytical methodology of the two-dimensional (2D) unsteady groundwater flow systems in Tóth's classical model, which will provide new insights into understanding the groundwater flow systems.

Hydraulic Head Solution under Steady-State Condition
In most of the previous studies, the media were assumed to be homogeneous and isotropic [25,26,30], and the heterogeneity has been considered in few studies [27,28].However, the porous medium is always anisotropic.The horizontal hydraulic conductivity is usually greater than the vertical one.In this paper, an analytical solution of the steadystate groundwater flow system in a small basin considering the anisotropy of hydraulic conductivity based on Tóth's classical model [2] (Figure 1) was derived.
The governing equation of the potential function for a 2D steady-state flow field considering the effect of anisotropy can be described as follows: The boundary conditions can be written as The following dimensionless variables are defined:   = /  ,   = /  , ℎ  = ℎ/  ,  = √  /  ,  =   /  , and  = /  .So the problem can be transformed to the following equations: The mathematical model can be solved by using the segregation variable method.The details can be found in Appendix A. The solution of the hydraulic head can be expressed as where If we set  = 1,   = /  ,   = /  , ℎ  = ℎ/  ,  =   /  ,  = /  ,  = 2/, and   = 4 in ( 4), (4) will be changed to Equation ( 6) is the same as the solution of Tóth [2].One can use the new solution equation ( 4) to investigate  the impact of anisotropy on groundwater flow systems.The new solution can also be used to study the unsteady flow.A MATLAB program [31] was developed to calculate (4).The parameters were given as   = 7000 m,   = 3500 m,  = 15, and tan  = 0.02.Then  = 2,  = 15/3500.According to Darcy's law, The dimensionless hydraulic head, ℎ  , and the velocity distribution were calculated with those parameters and different values of anisotropy,  = √  /  , as shown in Figure 2.
The area of equal velocity region decreases significantly with the anisotropy ratio (the contour interval is 0.003).The penetrating depth of groundwater flow systems decreases with the horizontal hydraulic conductivity, or the horizontal flow velocity increases with the horizontal hydraulic conductivity.Additionally, stagnant zones can be found from the simulating results.They are located at the places that have opposite direction for the streamlines (stagnant zone-1) and at both sides at the bottom in the field (stagnant zone-2).
The stagnant zone-1 moves toward ground surface as   increases.

Hydraulic Head Solution under
Unsteady-State Condition 3.1.Mathematical Model and Its Analytical Solutions.According to Wang and Wan [29], a 2D unsteady flow may be considered as a steady flow plus a relative head.The initial condition is ℎ(, ,  = 0) = ℎ 0 (ℎ 0 is a constant).Assuming ℎ(, , ) = ℎ  (, , ) + ℎ V (, , ), where ℎ  is the hydraulic head at the steady flow (4), ℎ V is relative head.Thus, the hydraulic head at unsteady flow can be obtained by solving ℎ V and then plus ℎ  .The mathematical model can be described as follows: The initial condition and boundary conditions can be written as where   (m −1 ) is specific yield and  (d) is the time.
The potential function can be obtained using segregation variable technique (the detailed derivation can be found in Appendix B): where Equation ( 15) is the same as (4) when the time is large enough (  > 1000, in Figure 4).
The hydraulic heads at the right side are significantly higher than those at the left side if a sinusoidally undulating water table is set at  > 0 (Figure 3).In addition, the distribution of ℎ is irregular at early times (Figure 3(a)).When time is large enough, the ultimate distribution of unsteady flow is similar to that of the steady state (Figures 3(c) and 2(a)).Figure 4 shows the ℎ-time behavior for three different positions A, B, and C (Figure 1).Each point desires a precise time to achieve a steady state indicating the flow from unsteady to steady, and the flow at the right side approaches steady state faster than that of the left side.

Sensitivity Analysis of the Parameters.
Sensitivity analysis is a tool to analyze the impact of the input parameters on the results of a model.It is essential to analyze the sensitivity of parameters in order to improve the accuracy of groundwater flow models and to reduce errors induced from the uncertainty of hydrogeological parameters [32].Studying the response of simulation results to the variations of important parameters can not only facilitate improving the reliability and accuracy of the models but also play a vital role in hydrogeological survey.
In unsteady groundwater flow systems, hydraulic head and specific yield can have great impacts on the evolution process of unsteady flow.The normalized sensitivity method [33,34] was used to analyze the sensitivity of the results to  and   .The normalized sensitivity of  (  =   ) and   to the relative modification of a given parameter can be expressed as follows: where ℎ (,) is ℎ of point (, ) changes with time before the parameters change and ℎ  (,) is ℎ of point (, ) changes with time after the parameters change.
The values of  and   were increased by 10%, 50%, and 100% to investigate their impacts on the hydraulic head ℎ at the middle of the domain (D, Figure 1) (Figure 5). Figure 5  shows that a relative increase in  produces a positive effect on hydraulic head and a relative increase in   produces a negative effect on hydraulic head at early times.Hydraulic head is sensitive to  and   at early times but insensitive to them at late times.

Summary and Conclusions
On the basis of Tóth's classical model, analytical solution of hydraulic head containing hydraulic conductivity under steady and unsteady flow conditions is obtained.From this study, the following conclusions can be drawn: (1) For the steady flow, the area of equal velocity region becomes much smaller with a larger anisotropy ratio, and the penetrating depth of groundwater flow systems becomes smaller.Stagnant zones locate at the places which have opposite directions for the streamlines and both sides at the bottom of the field.Finally, when   becomes larger, stagnant zones at the places with opposite streamlines become much closer to the surface.
(2) For the transient flow, when time is large enough, the ultimate distribution of unsteady flow is consistent with the steady groundwater flow system model.The closer it is to the right side, the faster it is for the flow to approach steady state.
(3) For the transient flow, a relative increase (e.g., 10%, 50%, and 100%) in K produces a positive impact on hydraulic head, and a relative increase in   produces a negative impact on hydraulic head at early times.(A.12) is the Fourier series; then Then one can obtain  0 ,   .

B. Analytical Solution of Potential Function for Transient Flow
Analytical solution of hydraulic head for ( 10) is presented as follows.

Figure 5 :
Figure 5: Time series of hydraulic head of the midpoint (D) of the model with different  and   .