# Analytical Study of Unsteady Nested Groundwater Flow Systems

Authors: Hong Niu; Xing Liang; Sheng-nan Ni; Zhang Wen
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

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 ε=Kx/Kz. 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.


---

## Body

## 1. 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.Furthermore, the physical mechanisms of groundwater flow system have been explained by Tóth [4] and Engelen and Jones [10]. Numerical simulations, including the concepts and methods [11–18] and afterword physical simulations [19, 20], have also been developed. The theory of groundwater flow systems has been applied to study many practical problems in science and engineering [21–24].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 steady-state 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.


## 2. 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 steady-state 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.Figure 1
Idealized cross section of a drainage basin’s valley flank of the sinusoidal water table [2].The governing equation of the potential function for a 2D steady-state flow field considering the effect of anisotropy can be described as follows:(1)∂∂xKx∂h∂x+∂∂zKz∂h∂z=0.The boundary conditions can be written as(2)Kx∂h∂xx=0,x=Lx=0,Kz∂h∂zz=0=0,hx,zz=Lz=Lz+xtan⁡α+asin⁡bx/cos⁡αcos⁡α,in which h (m) is the hydraulic head (h); Kx (m/d) and Kz (m/d) are the horizontal and vertical hydraulic conductivity, respectively; Lx (m) is the length of the basin; Lz (m) is the depth of the basin, and the datum level is chosen at the bottom of the basin; α represents the average slope of the valley flank; a and b represent the amplitude and frequency of the sinusoid, respectively; b=2π/λ, Lx=4λ, and λ represents the wavelength.The following dimensionless variables are defined:xD=x/Lx, zD=z/Lz, hD=h/Lz, ε=Kx/Kz, σ=Lx/Lz, and η=a/Lz. So the problem can be transformed to the following equations:(3)ε2σ2∂∂xD∂hD∂xD+∂∂zD∂hD∂zD=0,0≤xD≤1,0≤zD≤1,∂hD∂xDxD=0,xD=1=0,∂hD∂zDzD=0=0,hDxD,zDzD=1=1+xDσtan⁡α+ηsin⁡8πxD/cos⁡αcos⁡α.The mathematical model can be solved by using the segregation variable method. The details can be found in AppendixA. The solution of the hydraulic head can be expressed as(4)hDxD,zD=a0+∑n=1∞Dncos⁡nπxDcosh⁡zDnπεσ,where a0=1+σtanα/2+η/8π1-cos(8π/cosα), (5)Dn=2cosh⁡nπε/σσtan⁡αnπ2-1n-1+8η1--1ncos⁡8π/cos⁡απ64-n2cos2α.If we setε=1, xD=x/Lx, zD=z/Lz, hD=h/Lz, σ=Lx/Lz, η=a/Lz, b=2π/λ, and Lx=4λ in (4), (4) will be changed to(6)h=Lz+tan⁡α·Lx2+ab·Lx1-cos⁡b·Lxcos⁡α+2∑n=1∞ab1-cos⁡b·Lx/cos⁡αcos⁡nπcos2α·b/cos⁡α2-nπ/Lx2+tan⁡α·Lx2n2π2cos⁡nπ-1cos⁡nπx/Lxcosh⁡nπz/LxLx·cosh⁡nπLz/Lx.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 Lx=7000m, Lz=3500m, a=15, and tan⁡α=0.02. Then σ=2, η=15/3500. According to Darcy’s law,(7)Vx=-Kx∂hD∂xD,Vz=-Kz∂hD∂zD,V=Vx2+Vz2.The dimensionless hydraulic head,hD, and the velocity distribution were calculated with those parameters and different values of anisotropy, ε=Kx/Kz, as shown in Figure 2.Figure 2
Distribution of hydraulic head (line with dates) and velocity in basin with different anisotropy ratio ((a)ε=1, (b) ε=2, and (c) ε=5, Kz=10).
(a)
(b)
(c)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 asKx increases.


## 3. 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 h(x,z,t=0)=h0 (h0 is a constant). Assuming h(x,z,t)=hs(x,z,t)+hv(x,z,t), where hs is the hydraulic head at the steady flow (4), hv is relative head. Thus, the hydraulic head at unsteady flow can be obtained by solving hv and then plus hs. The mathematical model can be described as follows:(8)∂∂xKx∂hv∂x+∂∂zKz∂hv∂z=μs∂hv∂t.The initial condition and boundary conditions can be written as(9)hvx,z,tt=0=h0-hsx,z,Kx∂hv∂xx=0,x=Lx=0,Kz∂hv∂zz=0=0,hvx,z,tz=Lz,t>0=0,where μs (m−1) is specific yield and t (d) is the time.Similarly, some additional dimensionless variables can be defined:hDv=hv/Lz, tD=Kxt/Lx, ρ=μsLz, and hD0=h0/Lz. Thus, the problem can be transformed to the subsequent equations:(10)ε2σ2∂∂xD∂hDv∂xD+∂∂zD∂hDv∂zD=ε2ρσ∂hDv∂tD,0≤xD≤1,0≤zD≤1,(11)hDvxD,zD,tD=0=hD0-hDs,(12)∂hDv∂xDxD=0,xD=1=0,(13)∂hDv∂zDzD=0=0,(14)hDvxD,zD,tDzD=1,tD>0=0.The potential function can be obtained using segregation variable technique (the detailed derivation can be found in AppendixB):(15)hDxD,zD,tD=a0+∑n=1∞Dncos⁡nπxDcosh⁡zDnπεσ+∑n=0∞∑m=0∞An·Cmcos⁡nπxDcos⁡1+2mπ2zD·exp⁡-4nπ2ε2+1+2m2π2σ24ε2ρσtD,where(16)An·Cm=4hD0-a0·-1m1+2mπ,n=0;m=0,1,2,3,…,-4Dncoshnπε/σ1+2m-1m4πn2ε2/σ2+1+2m2π,n=1,2,3,…;m=0,1,2,….Equation (15) is the same as (4) when the time is large enough (tD>1000, in Figure 4).


### 3.2. Results and Discussions

A MATLAB program was developed to calculate (15). The parameters are given as Lx=7000 m, Lz=3500 m, a=15, tan⁡α=0.02, h0=3500 m, μs=0.3 m−1, and Kx=Kz=10 m/d. So ε=1, σ=2, ρ=1050, hD0=1, and η=15/3500. When time is large enough, the flow approaches quasi-steady state, which means that the hydraulic head does not change with time (Figure 3(c)).Figure 3
Distribution of hydraulic head (line with dates) and velocity in basin at different periods: (a)tD=143, (b) tD=286, and (c) tD=1143.
(a)
(b)
(c)Figure 4
Time series of hydraulic head in different nodes of the upper boundary.The hydraulic heads at the right side are significantly higher than those at the left side if a sinusoidally undulating water table is set att>0 (Figure 3). In addition, the distribution of h 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 h-time behavior for three different positionsA,B, andC (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.


### 3.3. 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 K and μs. The normalized sensitivity of K(Kx=Kz) and μs to the relative modification of a given parameter can be expressed as follows:(17)βK=KhDi,j′-hDi,jΔK,βμs=μshDi,j′-hDi,jΔμs,where hD(i,j) is h of point (i,j) changes with time before the parameters change and hD(i,j)′ is h of point (i,j) changes with time after the parameters change.The values ofK and μs were increased by 10%, 50%, and 100% to investigate their impacts on the hydraulic head h at the middle of the domain (D, Figure 1) (Figure 5). Figure 5 shows that a relative increase in K produces a positive effect on hydraulic head and a relative increase in μs produces a negative effect on hydraulic head at early times. Hydraulic head is sensitive to K and μs at early times but insensitive to them at late times.Figure 5
Time series of hydraulic head of the midpoint (D) of the model with different K and μs.


## 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 h(x,z,t=0)=h0 (h0 is a constant). Assuming h(x,z,t)=hs(x,z,t)+hv(x,z,t), where hs is the hydraulic head at the steady flow (4), hv is relative head. Thus, the hydraulic head at unsteady flow can be obtained by solving hv and then plus hs. The mathematical model can be described as follows:(8)∂∂xKx∂hv∂x+∂∂zKz∂hv∂z=μs∂hv∂t.The initial condition and boundary conditions can be written as(9)hvx,z,tt=0=h0-hsx,z,Kx∂hv∂xx=0,x=Lx=0,Kz∂hv∂zz=0=0,hvx,z,tz=Lz,t>0=0,where μs (m−1) is specific yield and t (d) is the time.Similarly, some additional dimensionless variables can be defined:hDv=hv/Lz, tD=Kxt/Lx, ρ=μsLz, and hD0=h0/Lz. Thus, the problem can be transformed to the subsequent equations:(10)ε2σ2∂∂xD∂hDv∂xD+∂∂zD∂hDv∂zD=ε2ρσ∂hDv∂tD,0≤xD≤1,0≤zD≤1,(11)hDvxD,zD,tD=0=hD0-hDs,(12)∂hDv∂xDxD=0,xD=1=0,(13)∂hDv∂zDzD=0=0,(14)hDvxD,zD,tDzD=1,tD>0=0.The potential function can be obtained using segregation variable technique (the detailed derivation can be found in AppendixB):(15)hDxD,zD,tD=a0+∑n=1∞Dncos⁡nπxDcosh⁡zDnπεσ+∑n=0∞∑m=0∞An·Cmcos⁡nπxDcos⁡1+2mπ2zD·exp⁡-4nπ2ε2+1+2m2π2σ24ε2ρσtD,where(16)An·Cm=4hD0-a0·-1m1+2mπ,n=0;m=0,1,2,3,…,-4Dncoshnπε/σ1+2m-1m4πn2ε2/σ2+1+2m2π,n=1,2,3,…;m=0,1,2,….Equation (15) is the same as (4) when the time is large enough (tD>1000, in Figure 4).


## 3.2. Results and Discussions

A MATLAB program was developed to calculate (15). The parameters are given as Lx=7000 m, Lz=3500 m, a=15, tan⁡α=0.02, h0=3500 m, μs=0.3 m−1, and Kx=Kz=10 m/d. So ε=1, σ=2, ρ=1050, hD0=1, and η=15/3500. When time is large enough, the flow approaches quasi-steady state, which means that the hydraulic head does not change with time (Figure 3(c)).Figure 3
Distribution of hydraulic head (line with dates) and velocity in basin at different periods: (a)tD=143, (b) tD=286, and (c) tD=1143.
(a)
(b)
(c)Figure 4
Time series of hydraulic head in different nodes of the upper boundary.The hydraulic heads at the right side are significantly higher than those at the left side if a sinusoidally undulating water table is set att>0 (Figure 3). In addition, the distribution of h 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 h-time behavior for three different positionsA,B, andC (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.


## 3.3. 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 K and μs. The normalized sensitivity of K(Kx=Kz) and μs to the relative modification of a given parameter can be expressed as follows:(17)βK=KhDi,j′-hDi,jΔK,βμs=μshDi,j′-hDi,jΔμs,where hD(i,j) is h of point (i,j) changes with time before the parameters change and hD(i,j)′ is h of point (i,j) changes with time after the parameters change.The values ofK and μs were increased by 10%, 50%, and 100% to investigate their impacts on the hydraulic head h at the middle of the domain (D, Figure 1) (Figure 5). Figure 5 shows that a relative increase in K produces a positive effect on hydraulic head and a relative increase in μs produces a negative effect on hydraulic head at early times. Hydraulic head is sensitive to K and μs at early times but insensitive to them at late times.Figure 5
Time series of hydraulic head of the midpoint (D) of the model with different K and μs.


## 4. 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, whenKx 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%) inK produces a positive impact on hydraulic head, and a relative increase in μs produces a negative impact on hydraulic head at early times.



---

*Source: 284181-2015-10-05.xml*