Analysis of three-dimensional unsaturated–saturated flow induced by localized recharge in unconfined aquifers

In the process of groundwater recharge, surface water usually enters an aquifer by passing an overlying unsaturated zone. Little attention has been given to the development of analytical solutions to a coupled unsaturated–saturated flow model due to localized recharge up to now. This paper develops a mathematical model to depict three-dimensional transient unsaturated–saturated flow in an unconfined aquifer with localized recharge on the ground surface. The model contains Richards' equation for unsaturated flow, a flow equation for saturated formation, and the Gardner constitutive model describing the behavior of unsaturated soil properties. Both flow equations are coupled through the continuity conditions of the head and flux at the water table. The semi-analytical solution to the coupled flow model is derived by the methods of Laplace transform and Fourier cosine transform. A sensitivity analysis is performed to explore the head response to the change in each of the aquifer parameters. A quantitative tool is presented to assess the recharge efficiency signifying the percentage of the water from the recharge to the aquifer. We found that the effect of unsaturated flow on the saturated hydraulic head is negligible if two criteria associated with the unsaturated soil properties and initial aquifer thickness are satisfied. The head distributions predicted from the present solution match well with those from finite-difference simulations. The predictions of the present solution also agree well with the observed data from a field experiment at an artificial recharge pond in Fresno County, California.

1 Introduction

Understanding the effect of water flow due to recharge from a surface water body such as precipitation, lake, or artificial pond on the groundwater flow system is important in water resource planning and management (e.g., Wang et al., 2010; Siltecho et al., 2015; Yang et al., 2015; Scudeler et al., 2016; Shi et al., 2016). The subsurface soil formation may be divided into unsaturated and saturated zones depending on the water saturation in void spaces of the soils. In the recharge process, the surface water may infiltrate and flow through the unsaturated zone and then arrives at the water table of the saturated zone (i.e., aquifers). Chang et al. (2016) reviewed analytical solutions describing the spatiotemporal distributions of groundwater mounds caused by localized recharge on the ground surface. They classified 17 solutions in a tabular form with flow dimensions as well as six headings of references, aquifer domain, aquifer boundary conditions, recharge region, recharge rate, and remarks. However, those solutions they reviewed all neglect the process of infiltration in the unsaturated zone and assume that the surface water directly recharges the saturated zone.

Solving Richards' equation (Richards, 1931) analytically for unsaturated flow is tricky owing to its nonlinearity. Gardner (1958) presented a model to express the relative hydraulic conductivity as an exponential function of the pressure head in unsaturated soils. Analytical methods for developing the solution to Richards' equation mostly rely on the use of linearization based on Gardner's model. Many articles used such an approach to study flow in an unsaturated zone with infiltration from a variety of surface water bodies (see, e.g., Huang and Wu, 2012; Wu et al., 2013; Wang and Li, 2015). Those articles neglected the presence of an underlying aquifer and treated its water table as the lower boundary with a condition of constant pressure head (e.g., Huang and Wu, 2012) or water content (e.g., Chen et al., 2001b). For 1-D downward flow, Srivastava and Yeh (1991) discussed distributions of the pressure head and water content in two distinct unsaturated soil layers with a constant surface flux. Chen et al. (2001a) examined the water content in an unsaturated medium with an arbitrary time-varying surface flux, by extending Warrick's (1975) solution for a flux consisting of step functions of time. Later, Wu et al. (2012) did similar work to Srivastava and Yeh (1991), but additionally considered the deformation of the two-layer soils caused by the change in the porewater pressure in the soils due to the surface flux. For 2-D flow in a vertical plane, Batu (1980) analyzed steady-state flow net affected by an array of strip surface sources with two different infiltration rates. Protopapas and Bras (1991) focused on the transient pressure head due to a uniform strip source with a finite width and infinite length on the ground surface. For 3-D flow, Chen et al. (2001b) investigated the water content induced by a surface source with an arbitrary spatiotemporal infiltration rate. Tracy (2007) studied the pressure head distribution in a cuboid soil sample with localized recharge over a rectangular area on the top. The sides of the sample are under either the Dirichlet or no-flow boundary condition.

Abovementioned solutions are applicable to either the case of saturated flow in aquifers recharged directly by surface water or the case of unsaturated flow due to surface water infiltration. So far little has been known about the combination of saturated and unsaturated flows that represents a typical process of recharge to the aquifer. This paper aims at developing a mathematical model for describing 3-D transient unsaturated–saturated flow in an unconfined aquifer with localized recharge. Richards' equation along with Gardner's model is adopted to delineate unsaturated flow between the ground surface and the water table. The 3-D groundwater flow equation is employed to depict saturated flow in the aquifer. Richards' equation is coupled with the saturated flow equation via the continuity conditions of the head and flux at the water table. Such a coupled flow model has been proposed by several articles to investigate pumping drawdown problems (e.g., Mathias and Butler, 2006; Tartakovsky and Neuman, 2007; Mishra and Neuman, 2010, 2011). They treated an extraction well as a line sink in the aquifer, while we consider the localized recharge as a plane source to the aquifer. The coupled flow model in their studies is 2-D written in cylindrical coordinates, while that in ours is 3-D expressed in Cartesian coordinates. In addition, their solutions are obtained by the Hankel transform, but ours is based on the Fourier cosine transform. The present work aims to investigate the spatiotemporal distribution of the hydraulic head due to localized recharge from the ground surface. The semi-analytical solution for the hydraulic head is obtained by the Laplace transform and the Fourier cosine transform. A finite-difference solution is built to check the correctness of the present solution. The effect of the unsaturated zone on the head in the saturated aquifer is explored by the present solution. The water quantity from the localized recharge to the aquifer is analyzed. The sensitivity analysis is executed to examine the head response to the variation in each of the aquifer parameters. Application of the present solution to a field experiment of artificial recharge is also provided.

2 Methodology

2.1 Mathematical model

Consider an unconfined aquifer system with localized recharge over a rectangular area on the ground surface of the system. The origin of the Cartesian coordinate system is located at the center of the recharge area as illustrated in Fig. 1a. The area has a size of 2l by 2w on the x−y plane. The shortest distance between an observation point (x, y) and a point (x_e, y_e) on the edge of the area is defined as $d=min\left(\sqrt{(x-x_{\mathrm{e}}{)}^{\mathrm{2}}+(y-y_{\mathrm{e}}{)}^{\mathrm{2}}}\right)$. The initial water table separates the unsaturated and saturated zones as shown in Fig. 1b and is chosen as the reference datum of the coordinate system. The initial thicknesses of the unsaturated and saturated zones prior to the recharge are denoted as b and B, respectively.

Figure 1 Schematic diagram of unsaturated–saturated flow in an unconfined aquifer system with localized recharge: (a) top view and (b) cross-sectional view.

The mathematical model for the aquifer system comprises two simultaneous equations for unsaturated and saturated flows. The equation for saturated flow in homogeneous and anisotropic aquifers is expressed as

$$\begin{array}{}\text{(1)}& {\displaystyle }{K}_x{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial x^{\mathrm{2}}}}+K_y{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial y^{\mathrm{2}}}}+K_z{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial z^{\mathrm{2}}}}=S_{\mathrm{s}}{\displaystyle \frac{\partial h}{\partial t}}\text{for}-B\le z\le \mathrm{0},\end{array}$$

where h(x, y, z, t) is the hydraulic head in the saturated zone; t is elapsed time since recharge began; K_x, K_y, and K_z are, respectively, the saturated hydraulic conductivities in the x, y, and z directions; S_s is the specific storage. Richards' equation for unsaturated flow is expressed as (Richards, 1931)

$$\begin{array}{ll}{\displaystyle }{K}_x& {\displaystyle \frac{\partial }{\partial x}}\left[k_{\mathrm{r}}\left(\mathit{\varphi }\right){\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}\right]+K_y{\displaystyle \frac{\partial }{\partial y}}\left[k_{\mathrm{r}}\left(\mathit{\varphi }\right){\displaystyle \frac{\partial \mathit{\varphi }}{\partial y}}\right]\\ \text{(2)}& {\displaystyle }& {\displaystyle }+K_z{\displaystyle \frac{\partial }{\partial z}}\left[k_{\mathrm{r}}\left(\mathit{\varphi }\right){\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\right]=C\left(\mathit{\varphi }\right){\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}},\end{array}$$

where ϕ(x, y, z, t) is the hydraulic head in the unsaturated zone. The relative hydraulic conductivity k_r(ϕ) and specific moisture capacity C(ϕ) are defined by the Gardner constitutive model (Gardner, 1958) as

$$\begin{array}{}\text{(3)}& {\displaystyle }{k}_{\mathrm{r}}\left(\mathit{\varphi }\right)=e^{a(\mathit{\varphi }-z)}\end{array}$$

and

$$\begin{array}{}\text{(4)}& {\displaystyle }C\left(\mathit{\varphi }\right)=a{S}_y{e}^{a(\mathit{\varphi }-z)},\end{array}$$

where S_y is the specific yield and a is the unsaturated exponent related to the pore-size distribution of a medium ranging from 0.2 to 5 m⁻¹ (Philip, 1969). Substituting Eqs. (3) and (4) into Eq. (2) leads to

$$\begin{array}{ll}{\displaystyle }{K}_x& {\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x^{\mathrm{2}}}}+K_y{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial y^{\mathrm{2}}}}+K_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial z^{\mathrm{2}}}}-a{K}_z{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\\ {\displaystyle }& {\displaystyle }+a\left[K_x{\left({\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}\right)}^{\mathrm{2}}+K_y{\left({\displaystyle \frac{\partial \mathit{\varphi }}{\partial y}}\right)}^{\mathrm{2}}+K_z{\left({\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\right)}^{\mathrm{2}}\right]\\ \text{(5)}& {\displaystyle }& {\displaystyle }=a{S}_y{\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}.\end{array}$$

It is essentially nonlinear and solved with difficulty by analytical methods. Kroszynski and Dagan (1975) employed the approach of perturbation expansion to simplify Richards' equation as a first-order linearized equation and developed an approximate solution for unsaturated–saturated flow induced by well pumping. The approach is extensively used in many studies on unsaturated–saturated flow (e.g., Mathias and Butler, 2006; Tartakovsky and Neuman, 2007; Mishra et al., 2012; Liang et al., 2017a). The linearized version of Richards' equation is written as

$$\begin{array}{ll}{\displaystyle }{K}_x& {\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x^{\mathrm{2}}}}+K_y{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial y^{\mathrm{2}}}}+K_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial z^{\mathrm{2}}}}-a{K}_z{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\\ \text{(6)}& {\displaystyle }& {\displaystyle }=a{S}_y{\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}\text{for}\mathrm{0}\le z\le b.\end{array}$$

The initial conditions for those two zones are

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathit{\varphi }=h=\mathrm{0}\text{at}t=\mathrm{0}.\end{array}$$

Because of symmetry of the recharge area along the x and y axes, the first quadrant (i.e., x≥0 and y≥0) of the flow domain is considered. Thus, all the horizontal outer boundaries are specified as the no-flow condition expressed as

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathrm{u}}}={\displaystyle \frac{\partial h}{\partial u}}=\mathrm{0}\text{at}u=\mathrm{0},\text{(9)}& {\displaystyle }& {\displaystyle }{lim}_{u\to \mathrm{\infty }}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial u}}={lim}_{u\to \mathrm{\infty }}{\displaystyle \frac{\partial h}{\partial u}}=\mathrm{0},\end{array}$$

where $u\in (x,y)$. The top boundary condition for the recharge area is denoted as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{K}_z{e}^{-az}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}=I[H\left(x\right)-H\left(x-l\right)][H\left(y\right)-H\left(y-w\right)]\\ \text{(10)}& {\displaystyle }& {\displaystyle }\text{at}z=b,\end{array}$$

where I is a constant recharge rate and H() is the Heaviside step function. Note that Eq. (10) can be written as $K_z\exp (-az)\partial \mathit{\varphi }/\partial z=I$ inside the recharge area $\mathrm{0}\le x\le l$ and $\mathrm{0}\le y\le w$ and denoted as $\partial \mathit{\varphi }/\partial z=\mathrm{0}$ outside that area. The impermeable boundary condition at the aquifer bottom is written as

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial h}{\partial z}}=\mathrm{0}\text{at}z=-B.\end{array}$$

The two continuity requirements of the hydraulic head and flux at the water table are expressed, respectively, as

$$\begin{array}{}\text{(12)}& {\displaystyle }\mathit{\varphi }=h\text{at}z=\mathrm{0}\end{array}$$

and

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}={\displaystyle \frac{\partial h}{\partial z}}\text{at}z=\mathrm{0}.\end{array}$$

The continuity conditions are valid when the water table change is less than 50 % of the initial saturated aquifer thickness, which is certified by a Hele-Shaw experiment (Marino, 1967).

Define the dimensionless variables and parameters as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{h}={\displaystyle \frac{h}{B}},\stackrel{\mathrm{‾}}{\mathit{\varphi }}={\displaystyle \frac{\mathit{\varphi }}{B}},\stackrel{\mathrm{‾}}{x}={\displaystyle \frac{x}{d}},\stackrel{\mathrm{‾}}{y}={\displaystyle \frac{y}{d}},\stackrel{\mathrm{‾}}{z}={\displaystyle \frac{z}{B}},\stackrel{\mathrm{‾}}{t}={\displaystyle \frac{K_{x}t}{S_{\mathrm{s}}{d}^{\mathrm{2}}}},\\ {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{l}={\displaystyle \frac{l}{d}},\stackrel{\mathrm{‾}}{w}={\displaystyle \frac{w}{d}},\stackrel{\mathrm{‾}}{b}={\displaystyle \frac{b}{B}},{\mathit{\kappa }}_z={\displaystyle \frac{K_z{d}^{\mathrm{2}}}{K_x{B}^{\mathrm{2}}}},{\mathit{\kappa }}_y={\displaystyle \frac{K_y}{K_x}},\\ \text{(14)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }=aB,\mathit{\sigma }={\displaystyle \frac{S_y}{S_{\mathrm{s}}B}},\mathit{\xi }={\displaystyle \frac{I}{K_z}},\end{array}$$

where the overbar represents a dimensionless variable or parameter. According to Eq. (14), the unsaturated–saturated flow model is rewritten as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial {\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}}}+{\mathit{\kappa }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial {\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}}}+{\mathit{\kappa }}_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial {\stackrel{\mathrm{‾}}{z}}^{\mathrm{2}}}}-\mathit{\alpha }{\mathit{\kappa }}_z{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial \stackrel{\mathrm{‾}}{z}}}\\ \text{(15)}& {\displaystyle }& {\displaystyle }=\mathit{\alpha }\mathit{\sigma }{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial \stackrel{\mathrm{‾}}{t}}}\text{for}\mathrm{0}\le \stackrel{\mathrm{‾}}{z}\le \stackrel{\mathrm{‾}}{b}\text{(16)}& {\displaystyle }& {\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{‾}}{h}}{\partial {\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}}}+{\mathit{\kappa }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{‾}}{h}}{\partial {\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}}}+{\mathit{\kappa }}_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{‾}}{h}}{\partial {\stackrel{\mathrm{‾}}{z}}^{\mathrm{2}}}}={\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{h}}{\partial \stackrel{\mathrm{‾}}{t}}}\text{for}-\mathrm{1}\le \stackrel{\mathrm{‾}}{z}\le \mathrm{0}\text{(17)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{\mathit{\varphi }}=\stackrel{\mathrm{‾}}{h}=\mathrm{0}\text{at}\stackrel{\mathrm{‾}}{t}=\mathrm{0},\text{(18)}& {\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial \stackrel{\mathrm{‾}}{u}}}={\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{h}}{\partial \stackrel{\mathrm{‾}}{u}}}=\mathrm{0}\text{at}\stackrel{\mathrm{‾}}{u}=\mathrm{0},\text{(19)}& {\displaystyle }& {\displaystyle }{lim}_{\stackrel{\mathrm{‾}}{u}\to \mathrm{\infty }}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial \stackrel{\mathrm{‾}}{u}}}={lim}_{\stackrel{\mathrm{‾}}{u}\to \mathrm{\infty }}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{h}}{\partial \stackrel{\mathrm{‾}}{u}}}=\mathrm{0},{\displaystyle }& {\displaystyle }\exp (-\mathit{\alpha }\stackrel{\mathrm{‾}}{z}){\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial \stackrel{\mathrm{‾}}{z}}}=\mathit{\xi }[H\left(\stackrel{\mathrm{‾}}{x}\right)-H\left(\stackrel{\mathrm{‾}}{x}-\stackrel{\mathrm{‾}}{l}\right)][H\left(\stackrel{\mathrm{‾}}{y}\right)-H\left(\stackrel{\mathrm{‾}}{y}-\stackrel{\mathrm{‾}}{w}\right)]\\ \text{(20)}& {\displaystyle }& {\displaystyle }\text{at}\stackrel{\mathrm{‾}}{z}=\stackrel{\mathrm{‾}}{b},\text{(21)}& {\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{h}}{\partial \stackrel{\mathrm{‾}}{z}}}=\mathrm{0}\text{at}\stackrel{\mathrm{‾}}{z}=-\mathrm{1},\text{(22)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{\mathit{\varphi }}=\stackrel{\mathrm{‾}}{h}\text{at}\stackrel{\mathrm{‾}}{z}=\mathrm{0},\text{(23)}& {\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial \stackrel{\mathrm{‾}}{z}}}={\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{h}}{\partial \stackrel{\mathrm{‾}}{z}}}\text{at}\stackrel{\mathrm{‾}}{z}=\mathrm{0},\end{array}$$

where $\stackrel{\mathrm{‾}}{u}\in (\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y})$.

2.2 Laplace domain solution

The unsaturated–saturated flow model composed of Eqs. (15)–(23) is solved by the methods of Laplace and Fourier cosine transforms. The Laplace transform is defined as

$$\begin{array}{}\text{(24)}& {\displaystyle }\stackrel{\mathrm{̃}}{f}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\stackrel{\mathrm{‾}}{f}\exp (-p\stackrel{\mathrm{‾}}{t})\mathrm{d}\stackrel{\mathrm{‾}}{t}\end{array}$$

with the property that

$$\begin{array}{}\text{(25)}& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{f}}{\partial \stackrel{\mathrm{‾}}{t}}}\exp (-p\stackrel{\mathrm{‾}}{t})\mathrm{d}\stackrel{\mathrm{‾}}{t}=p\stackrel{\mathrm{̃}}{f}-\stackrel{\mathrm{‾}}{f}{\mathrm{|}}_{\stackrel{\mathrm{‾}}{t}=\mathrm{0}},\end{array}$$

where $\stackrel{\mathrm{̃}}{f}\in (\stackrel{\mathrm{̃}}{\mathit{\varphi }},\stackrel{\mathrm{̃}}{h})$ represents the dimensionless hydraulic head in the Laplace domain, p is the Laplace transform parameter, $\stackrel{\mathrm{‾}}{f}\in (\stackrel{\mathrm{‾}}{\mathit{\varphi }},\stackrel{\mathrm{‾}}{h})$, and $\stackrel{\mathrm{‾}}{f}{\mathrm{|}}_{\stackrel{\mathrm{‾}}{t}=\mathrm{0}}=\mathrm{0}$ is from Eq. (17). Using Eqs. (24) and (25) converts $\stackrel{\mathrm{‾}}{\mathit{\varphi }}(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},\stackrel{\mathrm{‾}}{z},\stackrel{\mathrm{‾}}{t})$ into $\stackrel{\mathrm{̃}}{\mathit{\varphi }}(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},\stackrel{\mathrm{‾}}{z},p)$, $\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}/\partial \stackrel{\mathrm{‾}}{t}$ into $p\stackrel{\mathrm{̃}}{\mathit{\varphi }}$, $\stackrel{\mathrm{‾}}{h}(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},\stackrel{\mathrm{‾}}{z},\stackrel{\mathrm{‾}}{t})$ into $\stackrel{\mathrm{̃}}{h}(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},\stackrel{\mathrm{‾}}{z},p)$, $\partial \stackrel{\mathrm{‾}}{h}/\partial \stackrel{\mathrm{‾}}{t}$ into $p\stackrel{\mathrm{̃}}{h}$, and ξ into ξ∕p. The model then becomes

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\varphi }}}{\partial {\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}}}+{\mathit{\kappa }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\varphi }}}{\partial {\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}}}+{\mathit{\kappa }}_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\varphi }}}{\partial {\stackrel{\mathrm{‾}}{z}}^{\mathrm{2}}}}-\mathit{\alpha }{\mathit{\kappa }}_z{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\varphi }}}{\partial \stackrel{\mathrm{‾}}{z}}}=\mathit{\alpha }\mathit{\sigma }p\stackrel{\mathrm{̃}}{\mathit{\varphi }}\end{array}$$

and

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{h}}{\partial {\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}}}+{\mathit{\kappa }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{h}}{\partial {\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}}}+{\mathit{\kappa }}_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{h}}{\partial {\stackrel{\mathrm{‾}}{z}}^{\mathrm{2}}}}=p\stackrel{\mathrm{̃}}{h}\end{array}$$

subject to the transformed boundary conditions of $\underset{\stackrel{\mathrm{‾}}{u}\to \mathrm{0},\mathrm{\infty }}{lim}\partial \stackrel{\mathrm{̃}}{\mathit{\varphi }}/\partial \stackrel{\mathrm{‾}}{u}=\underset{\stackrel{\mathrm{‾}}{u}\to \mathrm{0},\mathrm{\infty }}{lim}\partial \stackrel{\mathrm{̃}}{h}/\partial \stackrel{\mathrm{‾}}{u}=\mathrm{0}$ with $\stackrel{\mathrm{‾}}{u}\in (\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y})$, $\exp \left(-\mathit{\alpha }\stackrel{\mathrm{‾}}{z}\right)\partial \stackrel{\mathrm{̃}}{\mathit{\varphi }}/\partial \stackrel{\mathrm{‾}}{z}=(\mathit{\xi }/p)[H\left(\stackrel{\mathrm{‾}}{x}\right)-H\left(\stackrel{\mathrm{‾}}{x}-\stackrel{\mathrm{‾}}{l}\right)][H\left(\stackrel{\mathrm{‾}}{y}\right)-H\left(\stackrel{\mathrm{‾}}{y}-\stackrel{\mathrm{‾}}{w}\right)]$ at $\stackrel{\mathrm{‾}}{z}=\stackrel{\mathrm{‾}}{b}$, and $\partial \stackrel{\mathrm{̃}}{h}/\partial \stackrel{\mathrm{‾}}{z}=\mathrm{0}$ at $\stackrel{\mathrm{‾}}{z}=-\mathrm{1}$. Moreover, the transformed continuity conditions are $\stackrel{\mathrm{̃}}{\mathit{\varphi }}=\stackrel{\mathrm{̃}}{h}$ and $\partial \stackrel{\mathrm{̃}}{\mathit{\varphi }}/\partial \stackrel{\mathrm{‾}}{z}=\partial \stackrel{\mathrm{̃}}{h}/\partial \stackrel{\mathrm{‾}}{z}$ at $\stackrel{\mathrm{‾}}{z}=\mathrm{0}$.

Afterward, one may take the double Fourier cosine transform that provides

$$\begin{array}{}\text{(28)}& {\displaystyle }\widehat{f}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\stackrel{\mathrm{̃}}{f}\cos \left({\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{‾}}{x}\right)\cos \left({\mathit{\omega }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{y}\right)\mathrm{d}\stackrel{\mathrm{‾}}{x}\mathrm{d}\stackrel{\mathrm{‾}}{y}\end{array}$$

and

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{f}}{\partial {\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}}}+{\mathit{\kappa }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{f}}{\partial {\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}}}\right)\cos \left({\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{‾}}{x}\right)\cos \left({\mathit{\omega }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{y}\right)\mathrm{d}\stackrel{\mathrm{‾}}{x}\mathrm{d}\stackrel{\mathrm{‾}}{y}\\ \text{(29)}& {\displaystyle }& {\displaystyle }=-({\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\kappa }}_y{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}})\widehat{f},\end{array}$$

where $\widehat{f}\in (\widehat{\mathit{\varphi }},\widehat{h})$ represents the dimensionless hydraulic head in the Fourier domain; ω₁ and ω₂ are the Fourier cosine transform parameters. The transform converts $\stackrel{\mathrm{̃}}{\mathit{\varphi }}(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},\stackrel{\mathrm{‾}}{z},p)$ into $\widehat{\mathit{\varphi }}({\mathit{\omega }}_{\mathrm{1}},{\mathit{\omega }}_{\mathrm{2}},\stackrel{\mathrm{‾}}{z},p)$, $\stackrel{\mathrm{̃}}{h}(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},\stackrel{\mathrm{‾}}{z},p)$ into $\widehat{h}({\mathit{\omega }}_{\mathrm{1}},{\mathit{\omega }}_{\mathrm{2}},\stackrel{\mathrm{‾}}{z},p)$,

${\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\varphi }}/\partial {\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}+{\mathit{\kappa }}_y({\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\varphi }}/\partial {\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}})$ into $-({\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\kappa }}_y{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}})\widehat{\mathit{\varphi }}$, ${\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{h}/\partial {\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}+{\mathit{\kappa }}_y({\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{h}/\partial {\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}})$ into $-({\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\kappa }}_y{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}})\widehat{h}$, and $(\mathit{\xi }/p)[H\left(\stackrel{\mathrm{‾}}{x}\right)-H\left(\stackrel{\mathrm{‾}}{x}-\stackrel{\mathrm{‾}}{l}\right)][H\left(\stackrel{\mathrm{‾}}{y}\right)-H\left(\stackrel{\mathrm{‾}}{y}-\stackrel{\mathrm{‾}}{w}\right)]$ into $\mathit{\xi }\sin \left({\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{‾}}{l}\right)\sin \left({\mathit{\omega }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{w}\right)/\left(p{\mathit{\omega }}_{\mathrm{1}}{\mathit{\omega }}_{\mathrm{2}}\right)$. Equations (26) and (27) hence become ordinary differential equations in terms of $\stackrel{\mathrm{‾}}{z}$ denoted, respectively, as

$$\begin{array}{}\text{(30)}& {\displaystyle }{\mathit{\kappa }}_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\widehat{\mathit{\varphi }}}{\partial {\stackrel{\mathrm{‾}}{z}}^{\mathrm{2}}}}-\mathit{\alpha }{\mathit{\kappa }}_z{\displaystyle \frac{\partial \widehat{\mathit{\varphi }}}{\partial \stackrel{\mathrm{‾}}{z}}}-\left(\mathit{\alpha }\mathit{\sigma }p+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\kappa }}_y{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}\right)\widehat{\mathit{\varphi }}=\mathrm{0}\end{array}$$

and

$$\begin{array}{}\text{(31)}& {\displaystyle }{\mathit{\kappa }}_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\widehat{h}}{\partial {\stackrel{\mathrm{‾}}{z}}^{\mathrm{2}}}}-\left(p+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\kappa }}_y{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}\right)\widehat{h}=\mathrm{0}.\end{array}$$

Similarly, the transformed boundary conditions are expressed as

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{\partial \widehat{h}}{\partial \stackrel{\mathrm{‾}}{z}}}=\mathrm{0}\text{at}\stackrel{\mathrm{‾}}{z}=-\mathrm{1}\end{array}$$

and

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\exp (-\mathit{\alpha }\stackrel{\mathrm{‾}}{z}){\displaystyle \frac{\partial \widehat{\mathit{\varphi }}}{\partial \stackrel{\mathrm{‾}}{z}}}=\mathit{\xi }\sin \left({\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{‾}}{l}\right)\sin \left({\mathit{\omega }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{w}\right)/\left(p{\mathit{\omega }}_{\mathrm{1}}{\mathit{\omega }}_{\mathrm{2}}\right)\\ \text{(33)}& {\displaystyle }& {\displaystyle }\text{at}\stackrel{\mathrm{‾}}{z}=\stackrel{\mathrm{‾}}{b}.\end{array}$$

The transformed continuity conditions are written as

$$\begin{array}{}\text{(34)}& {\displaystyle }\widehat{\mathit{\varphi }}=\widehat{h}\text{at}\stackrel{\mathrm{‾}}{z}=\mathrm{0}\end{array}$$

and

$$\begin{array}{}\text{(35)}& {\displaystyle \frac{\partial \widehat{\mathit{\varphi }}}{\partial \stackrel{\mathrm{‾}}{z}}}={\displaystyle \frac{\partial \widehat{h}}{\partial \stackrel{\mathrm{‾}}{z}}}\text{at}\stackrel{\mathrm{‾}}{z}=\mathrm{0}.\end{array}$$

Solving Eqs. (30) and (31) subject to Eqs. (32)–(35) and then taking the inverse Fourier cosine transform leads to the solutions in the Laplace domain written as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\varphi }}\left(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},\stackrel{\mathrm{‾}}{z},p\right)={\displaystyle \frac{\mathrm{4}}{{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{R}_{\mathrm{e}}(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},{\mathit{\omega }}_{\mathrm{1}},{\mathit{\omega }}_{\mathrm{2}}){\mathrm{\Omega }}_{\mathit{\varphi }}({\mathit{\omega }}_{\mathrm{1}},{\mathit{\omega }}_{\mathrm{2}},\stackrel{\mathrm{‾}}{z},p)\\ \text{(36a)}& {\displaystyle }& {\displaystyle }\mathrm{d}{\mathit{\omega }}_{\mathrm{1}}\mathrm{d}{\mathit{\omega }}_{\mathrm{2}}\text{for}\mathrm{0}\le \stackrel{\mathrm{‾}}{z}\le \stackrel{\mathrm{‾}}{b}\end{array}$$

and

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{h}\left(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},\stackrel{\mathrm{‾}}{z},p\right)={\displaystyle \frac{\mathrm{4}}{{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{R}_{\mathrm{e}}(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y,}{\mathit{\omega }}_{\mathrm{1}},{\mathit{\omega }}_{\mathrm{2}}){\mathrm{\Omega }}_h({\mathit{\omega }}_{\mathrm{1}},{\mathit{\omega }}_{\mathrm{2}},\stackrel{\mathrm{‾}}{z},p)\\ \text{(36b)}& {\displaystyle }& {\displaystyle }\mathrm{d}{\mathit{\omega }}_{\mathrm{1}}\mathrm{d}{\mathit{\omega }}_{\mathrm{2}}\text{for}-\mathrm{1}\le \stackrel{\mathrm{‾}}{z}\le \mathrm{0}\end{array}$$

with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{\Omega }}_{\mathit{\varphi }}={\displaystyle \frac{\mathrm{4}}{p({\mathit{\mu }}_{\mathrm{1}}+{\mathit{\mu }}_{\mathrm{2}})}}\exp \left[{\displaystyle \frac{\mathit{\alpha }\stackrel{\mathrm{‾}}{z}+\left(\mathit{\alpha }+{\mathit{\lambda }}_{\mathrm{2}}\right)\stackrel{\mathrm{‾}}{b}}{\mathrm{2}}}\right]\\ {\displaystyle }& {\displaystyle }[{\mathit{\lambda }}_{\mathrm{2}}\cosh {\mathit{\lambda }}_{\mathrm{1}}\cosh \left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{z}}{\mathrm{2}}}\right)+\left(\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\sinh {\mathit{\lambda }}_{\mathrm{1}}-\mathit{\alpha }\cosh {\mathit{\lambda }}_{\mathrm{1}}\right)\\ \text{(36c)}& {\displaystyle }& {\displaystyle }\sinh \left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{z}}{\mathrm{2}}}\right)],\text{(36d)}& {\displaystyle }& {\displaystyle }{\mathrm{\Omega }}_h={\displaystyle \frac{\mathrm{4}}{p({\mathit{\mu }}_{\mathrm{1}}+{\mathit{\mu }}_{\mathrm{2}})}}{\mathit{\lambda }}_{\mathrm{2}}\exp \left[{\displaystyle \frac{(\mathit{\alpha }+{\mathit{\lambda }}_{\mathrm{2}})\stackrel{\mathrm{‾}}{b}}{\mathrm{2}}}\right]\cosh \left[\left(\mathrm{1}+\stackrel{\mathrm{‾}}{z}\right){\mathit{\lambda }}_{\mathrm{1}}\right],\text{(36e)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{e}}=\left\{\begin{array}{l}\mathit{\xi }\sin \left({\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{‾}}{l}\right)\sin \left({\mathit{\omega }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{w}\right)\cos \left({\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{‾}}{x}\right)\cos \left({\mathit{\omega }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{y}\right)/\left({\mathit{\omega }}_{\mathrm{1}}{\mathit{\omega }}_{\mathrm{2}}\right)\\ \mathrm{for}{\mathit{\omega }}_{\mathrm{1}}\ne \mathrm{0}\mathrm{and}{\mathit{\omega }}_{\mathrm{2}}\ne \mathrm{0},\\ \mathit{\xi }\stackrel{\mathrm{‾}}{w}\sin \left({\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{‾}}{l}\right)\cos \left({\mathit{\omega }}_{\mathrm{1}}\stackrel{\mathrm{‾}}{x}\right)/{\mathit{\omega }}_{\mathrm{1}}\mathrm{for}{\mathit{\omega }}_{\mathrm{1}}\ne \mathrm{0}\mathrm{and}{\mathit{\omega }}_{\mathrm{2}}=\mathrm{0},\\ \mathit{\xi }\stackrel{\mathrm{‾}}{l}\sin \left({\mathit{\omega }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{w}\right)\cos \left({\mathit{\omega }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{y}\right)/{\mathit{\omega }}_{\mathrm{2}}\mathrm{for}{\mathit{\omega }}_{\mathrm{1}}=\mathrm{0}\mathrm{and}{\mathit{\omega }}_{\mathrm{2}}\ne \mathrm{0},\\ \mathit{\xi }\stackrel{\mathrm{‾}}{w}\stackrel{\mathrm{‾}}{l}\mathrm{for}{\mathit{\omega }}_{\mathrm{1}}=\mathrm{0}\mathrm{and}{\mathit{\omega }}_{\mathrm{2}}=\mathrm{0},\end{array}\right.\text{(36f)}& {\displaystyle }& {\displaystyle }{\mathit{\mu }}_{\mathrm{1}}=[\exp \left(\stackrel{\mathrm{‾}}{b}{\mathit{\lambda }}_{\mathrm{2}}\right)-\mathrm{1}]\left({\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}-{\mathit{\alpha }}^{\mathrm{2}}\right)\cosh {\mathit{\lambda }}_{\mathrm{1}},\text{(36g)}& {\displaystyle }& {\displaystyle }{\mathit{\mu }}_{\mathrm{2}}=\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}[\left({\mathit{\lambda }}_{\mathrm{2}}+\mathit{\alpha }\right)\exp \left(\stackrel{\mathrm{‾}}{b}{\mathit{\lambda }}_{\mathrm{2}}\right)+{\mathit{\lambda }}_{\mathrm{2}}-\mathit{\alpha }]\sinh {\mathit{\lambda }}_{\mathrm{1}},{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{1}}=\sqrt{(p+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\kappa }}_y{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}})/{\mathit{\kappa }}_z},\\ \text{(36h)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{2}}=\sqrt{{\mathit{\alpha }}^{\mathrm{2}}+\mathrm{4}(\mathit{\alpha }\mathit{\sigma }p+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\kappa }}_y{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}})/{\mathit{\kappa }}_z}.\end{array}$$

Note that Eq. (36a) is the solution for unsaturated flow, while Eq. (36b) is that for saturated flow. The inverse Laplace transform to both solutions may not be tractable. The numerical inversion of Laplace transform proposed by Stehfest (1970) is therefore used to obtain time-domain results of the solutions. The double integrals in the solutions can be evaluated numerically by the Gaussian quadrature (e.g., Gerald and Wheatley, 2004) using the dblquad Matlab built-in function (Gilat and Subramaniam, 2007) or the NIntegrate Mathematica built-in function (Wolfram, 1996).

2.3 Solution for the transient recharge rate

The present solution can be applied to the problem of time-varying recharge rates based on Duhamel's integral (Bear, 1979, p. 158). The dimensionless transient head solution ${\stackrel{\mathrm{‾}}{g}}_{\mathrm{t}}$ subject to the dimensionless time-varying recharge rate ${\mathit{\xi }}_{\mathrm{t}}\left(\stackrel{\mathrm{‾}}{t}\right)$ can be expressed as

$$\begin{array}{}\text{(37)}& {\displaystyle }{\stackrel{\mathrm{‾}}{g}}_{\mathrm{t}}={\stackrel{\mathrm{‾}}{g}}_{\mathrm{0}}+\underset{\mathrm{0}}{\overset{\stackrel{\mathrm{‾}}{t}}{\int }}{\displaystyle \frac{\partial {\mathit{\xi }}_{\mathrm{t}}\left(\mathit{\tau }\right)}{\partial \mathit{\tau }}}\stackrel{\mathrm{‾}}{g}\left(\stackrel{\mathrm{‾}}{t}-\mathit{\tau }\right)\mathrm{d}\mathit{\tau },\end{array}$$

where τ is a dummy variable, ${\stackrel{\mathrm{‾}}{g}}_{\mathrm{0}}$ denotes $\stackrel{\mathrm{‾}}{\mathit{\varphi }}$ or $\stackrel{\mathrm{‾}}{h}$ for the initial dimensionless recharge rate ${\mathit{\xi }}_{\mathrm{t}}\left(\stackrel{\mathrm{‾}}{t}=\mathrm{0}\right)$, and $\stackrel{\mathrm{‾}}{g}(\stackrel{\mathrm{‾}}{t}-\mathit{\tau })$ represents $\stackrel{\mathrm{‾}}{\mathit{\varphi }}$ or $\stackrel{\mathrm{‾}}{h}$ with $\stackrel{\mathrm{‾}}{t}$ replaced by $\stackrel{\mathrm{‾}}{t}-\mathit{\tau }$. If Eq. (37) is not an integrable function, we can evaluate numerically through the discretization method that (Singh, 2005)

$$\begin{array}{}\text{(38a)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{g}}_N={\stackrel{\mathrm{‾}}{g}}_{\mathrm{0}}+\sum _{i=\mathrm{1}}^N{\displaystyle \frac{\mathrm{\Delta }{\mathit{\xi }}_i}{\mathrm{\Delta }\stackrel{\mathrm{‾}}{t}}}G(N-i+\mathrm{1}),\text{(38b)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{\mathit{\xi }}_i={\mathit{\xi }}_i-{\mathit{\xi }}_{i-\mathrm{1}},\text{(38c)}& {\displaystyle }& {\displaystyle }G\left(M\right)=\underset{\mathrm{0}}{\overset{\mathrm{\Delta }\stackrel{\mathrm{‾}}{t}}{\int }}\stackrel{\mathrm{‾}}{g}(M\mathrm{\Delta }\stackrel{\mathrm{‾}}{t}-\mathit{\tau })\mathrm{d}\mathit{\tau },\end{array}$$

where ${\stackrel{\mathrm{‾}}{g}}_N$ signifies the dimensionless head solution at $\stackrel{\mathrm{‾}}{t}=\mathrm{\Delta }\stackrel{\mathrm{‾}}{t}\times N$; $\mathrm{\Delta }\stackrel{\mathrm{‾}}{t}$ is a dimensionless time step; G(M) is called the ramp kernel; ξ_i and ξ_i−1 are, respectively, dimensionless recharge rates at $\stackrel{\mathrm{‾}}{t}=\mathrm{\Delta }\stackrel{\mathrm{‾}}{t}\times i$ and $\stackrel{\mathrm{‾}}{t}=\mathrm{\Delta }\stackrel{\mathrm{‾}}{t}\times (i-\mathrm{1})$.

Figure 2 Schematic diagram of finite-difference grids: (a) top view and (b) cross-sectional view.

2.4 Recharge efficiency

The percentage of the water from the localized recharge reaching the water table is defined as recharge efficiency (RE) (Munevar and Marino, 1999) written as

$$\begin{array}{}\text{(39)}& {\displaystyle }\mathrm{RE}\left(t\right)=K_z\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\partial h}{\partial z}}\mathrm{d}x\mathrm{d}y/\left(Ilw\right)\text{at}z=\mathrm{0},\end{array}$$

where the denominator $I\times l\times w$ is the volumetric rate of the water entering the aquifer system from the recharge, and the double integral is the sum of the infiltration flux at the water table. According to the dimensionless quantities defined in Eq. (14), Eq. (39) becomes

$$\begin{array}{}\text{(40)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathrm{RE}}\left(p\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{h}}{\partial \stackrel{\mathrm{‾}}{z}}}\mathrm{d}\stackrel{\mathrm{‾}}{x}\mathrm{d}\stackrel{\mathrm{‾}}{y}/\left(\mathit{\xi }\stackrel{\mathrm{‾}}{l}\stackrel{\mathrm{‾}}{w}\right)\text{at}\stackrel{\mathrm{‾}}{z}=\mathrm{0},\end{array}$$

where $\stackrel{\mathrm{̃}}{\mathrm{RE}}$ represents RE in the Laplace domain and $\stackrel{\mathrm{̃}}{h}$ is defined in Eq. (36b). The RE increases from zero to a value equal to or below unity. The infiltration process does not affect the water table when RE =0. On the other hand, the water from the surface recharge totally arrives at the aquifer when RE =1.

2.5 Sensitivity analysis

The sensitivity analysis is commonly used to assess the change in the hydraulic head in response to a small change in a hydraulic parameter. The normalized sensitivity coefficient based on the present solution is defined as

$$\begin{array}{}\text{(41)}& {\displaystyle }{S}_{i,t}={\displaystyle \frac{\partial O}{\partial P_i/P_i}},\end{array}$$

where O represents the present solution for the unsaturated or saturated flow and P_i is the ith parameter. Equation (41) can be approximated as

$$\begin{array}{}\text{(42)}& {\displaystyle }{S}_{i,t}={\displaystyle \frac{O\left(P_i+\mathrm{\Delta }{P}_i\right)-O\left(P_i\right)}{\mathrm{\Delta }{P}_i/P_i}},\end{array}$$

where ΔP_i is an increment set to 10⁻³P_i (Yeh et al., 2008). Note that a large value of $\left|S_{i,t}\right|$ indicates that the head is sensitive to the change in the target parameter.

2.6 Finite-difference solution

An iterative algorithm based on an implicit finite-difference approximation to Eq. (5) is developed to solve the nonlinear unsaturated–saturated flow model. Figure 2 shows the finite-difference grids in the simulation domains of $\mathrm{0}\le x\le \mathrm{500}$ m, $\mathrm{0}\le y\le \mathrm{500}$ m, and −20 m $\le z\le \mathrm{10}$ m discretized by a non-uniform grid with small grid sizes near the recharge area of $\mathrm{0}\le x\le \mathrm{50}$ m and $\mathrm{0}\le y\le \mathrm{50}$ m and large grid sizes away from that area. The domain falls in the first quadrant due to symmetrical flow to the x-axis and y-axis. The saturated thickness is 20 m and the unsaturated thickness is 10 m. All the boundaries except the recharge region are therefore under the no-flow condition. Equation (5) is approximated as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{K}_x\left({\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{w}}}{\mathrm{\Delta }X}}{\mathit{\varphi }}_{i+\mathrm{1},j,k}^{m+\mathrm{1}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{\Delta }{x}_{\mathrm{w}}\mathrm{\Delta }{x}_{\mathrm{e}}}}{\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1}}+{\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{e}}}{\mathrm{\Delta }X}}{\mathit{\varphi }}_{i-\mathrm{1},j,k}^{m+\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }+K_y\left({\displaystyle \frac{\mathrm{\Delta }{y}_{\mathrm{s}}}{\mathrm{\Delta }Y}}{\mathit{\varphi }}_{i,j+\mathrm{1},k}^{m+\mathrm{1}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{\Delta }{y}_{\mathrm{s}}\mathrm{\Delta }{y}_{\mathrm{n}}}}{\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1}}+{\displaystyle \frac{\mathrm{\Delta }{y}_{\mathrm{n}}}{\mathrm{\Delta }Y}}{\mathit{\varphi }}_{i,j-\mathrm{1},k}^{m+\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{K_z}{\mathrm{\Delta }{z}^{\mathrm{2}}}}\left({\mathit{\varphi }}_{i,j,k+\mathrm{1}}^{m+\mathrm{1}}-\mathrm{2}{\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1}}+{\mathit{\varphi }}_{i,j,k-\mathrm{1}}^{m+\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{a{K}_z}{\mathrm{\Delta }z}}\left({\mathit{\varphi }}_{i,j,k+\mathrm{1}}^{m+\mathrm{1}}-{\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }+a\left[K_x\left(\mathrm{\Delta }{x}_{\mathrm{w}}{\mathit{\varphi }}_{i+\mathrm{1},j,k}^{m+\mathrm{1}}-(\mathrm{\Delta }{x}_{\mathrm{w}}+\mathrm{\Delta }{x}_{\mathrm{e}}){\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1}}+\right.\right.\\ {\displaystyle }& {\displaystyle }{\left.\mathrm{\Delta }{x}_{\mathrm{e}}{\mathit{\varphi }}_{i-\mathrm{1},j,k}^{m+\mathrm{1}}\right)}^{\mathrm{2}}+K_y\left(\mathrm{\Delta }{y}_{\mathrm{s}}{\mathit{\varphi }}_{i,j+\mathrm{1},k}^{m+\mathrm{1}}-(\mathrm{\Delta }{y}_{\mathrm{s}}+\mathrm{\Delta }{y}_{\mathrm{n}}){\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1}}\right.\\ {\displaystyle }& {\displaystyle }\left.{\left.+\mathrm{\Delta }{y}_{\mathrm{n}}{\mathit{\varphi }}_{i,j-\mathrm{1},k}^{m+\mathrm{1}}\right)}^{\mathrm{2}}+{\displaystyle \frac{K_z}{\mathrm{\Delta }{z}^{\mathrm{2}}}}{\left({\mathit{\varphi }}_{i,j,k+\mathrm{1}}^{m+\mathrm{1}}-{\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1}}\right)}^{\mathrm{2}}\right]\\ \text{(43a)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{a{S}_y}{\mathrm{\Delta }t}}\left({\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1}}-{\mathit{\varphi }}_{i,j,k}^m\right),\text{(43b)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }X=\mathrm{\Delta }{x}_{\mathrm{w}}\mathrm{\Delta }{x}_{\mathrm{e}}(\mathrm{\Delta }{x}_{\mathrm{w}}+\mathrm{\Delta }{x}_{\mathrm{e}})/\mathrm{2},\text{(43c)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }Y=\mathrm{\Delta }{y}_{\mathrm{s}}\mathrm{\Delta }{y}_{\mathrm{n}}(\mathrm{\Delta }{y}_{\mathrm{s}}+\mathrm{\Delta }{y}_{\mathrm{n}})/\mathrm{2},\end{array}$$

where ${\mathit{\varphi }}_{i,j,k}^m$ is the hydraulic head in the unsaturated zone at a nodal point (i, j, k); superscript m represents one time step earlier than the present time denoted as superscript m+1; Δx_w, Δx_e, Δy_n and Δy_s are grid sizes beside a nodal point (i, j, k) in the west, east, north and south, respectively; Δz is the grid size on the z-axis; Δt is the time step. Note that Eq. (43) reduces to the discretized expression of Eq. (6) when the quadratic terms are neglected. Similarly, Eq. (1) is approximated as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{K}_x\left({\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{w}}}{\mathrm{\Delta }X}}{h}_{i+\mathrm{1},j,k}^{m+\mathrm{1}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{\Delta }{x}_{\mathrm{w}}\mathrm{\Delta }{x}_{\mathrm{e}}}}{h}_{i,j,k}^{m+\mathrm{1}}+{\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{e}}}{\mathrm{\Delta }X}}{h}_{i-\mathrm{1},j,k}^{m+\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }+K_y\left({\displaystyle \frac{\mathrm{\Delta }{y}_{\mathrm{s}}}{\mathrm{\Delta }Y}}{h}_{i,j+\mathrm{1},k}^{m+\mathrm{1}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{\Delta }{y}_{\mathrm{s}}\mathrm{\Delta }{y}_{\mathrm{n}}}}{h}_{i,j,k}^{m+\mathrm{1}}+{\displaystyle \frac{\mathrm{\Delta }{y}_{\mathrm{n}}}{\mathrm{\Delta }Y}}{h}_{i,j-\mathrm{1},k}^{m+\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{K_z}{\mathrm{\Delta }{z}^{\mathrm{2}}}}\left(h_{i,j,k+\mathrm{1}}^{m+\mathrm{1}}-\mathrm{2}{h}_{i,j,k}^{m+\mathrm{1}}+h_{i,j,k-\mathrm{1}}^{m+\mathrm{1}}\right)\\ \text{(44)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{S_{\mathrm{s}}}{\mathrm{\Delta }t}}\left(h_{i,j,k}^{m+\mathrm{1}}-h_{i,j,k}^m\right),\end{array}$$

where $h_{i,j,k}^m$ is the hydraulic head in the saturated zone at a nodal point (i, j, k). The initial condition for each nodal point is expressed as

$$\begin{array}{}\text{(45)}& {\displaystyle }{\mathit{\varphi }}_{i,j,k}^{\mathrm{1}}=h_{i,j,k}^{\mathrm{1}}=\mathrm{0}\text{ at each}(i,j,k).\end{array}$$

The no-flow condition specified at the outer boundaries shown in Fig. 2a and the bottom can be written as

$$\begin{array}{}\text{(46)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{i-\mathrm{1},j,k}^{m+\mathrm{1}}={\mathit{\varphi }}_{i+\mathrm{1},j,k}^{m+\mathrm{1}}\text{and}{h}_{i-\mathrm{1},j,k}^{m+\mathrm{1}}=h_{i+\mathrm{1},j,k}^{m+\mathrm{1}}\text{at}i=\mathrm{1},n_x,\text{(47)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{i,j-\mathrm{1},k}^{m+\mathrm{1}}={\mathit{\varphi }}_{i,j+\mathrm{1},k}^{m+\mathrm{1}}\text{and}{h}_{i,j-\mathrm{1},k}^{m+\mathrm{1}}=h_{i,j+\mathrm{1},k}^{m+\mathrm{1}}\text{at}j=\mathrm{1},n_y,\text{(48)}& {\displaystyle }& {\displaystyle }{h}_{i,j,k-\mathrm{1}}^{m+\mathrm{1}}=h_{i,j,k+\mathrm{1}}^{m+\mathrm{1}}\text{at}k=\mathrm{1},\end{array}$$

where n_x and n_y are the total number of grids on the x- and y-axes, respectively. The top boundary condition is approximated as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left\{\begin{array}{l}{\mathit{\varphi }}_{i,j,k-\mathrm{1}}^{m+\mathrm{1}}={\mathit{\varphi }}_{i,j,k+\mathrm{1}}^{m+\mathrm{1}}\text{outside the recharge area},\\ {\displaystyle \frac{K_z{e}^{-ab}}{\mathrm{\Delta }z}}\left({\mathit{\varphi }}_{i,j,k+\mathrm{1}}^{m+\mathrm{1}}-{\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1}}\right)=I\text{inside the recharge area}\end{array}\right.\\ \text{(49)}& {\displaystyle }& {\displaystyle }\text{at}k=n_z,\end{array}$$

where n_z is the total number of grids on the z-axis. The grid sizes Δx_w, Δx_e, Δy_s, and Δy_n are all 5 m inside the recharge area, while outside the area they gradually increase according to the formulas Δx_e=1.2Δx_w and Δy_n=1.2Δy_s starting from $\mathrm{\Delta }{x}_{\mathrm{w}}=\mathrm{\Delta }{y}_{\mathrm{s}}=\mathrm{5}$ m and $\mathrm{\Delta }{x}_{\mathrm{e}}=\mathrm{\Delta }{y}_{\mathrm{n}}=\mathrm{1.2}\times \mathrm{5}$ m =6 m. Note that the largest grid size is set equal to 25 m for good accuracy in solution prediction. The grid size Δz is set to 0.1 m and the time step Δt is chosen as 0.1 days for the period of 0–2.5 days and 0.25 days for 2.5–100 days. The total number of nodal points is 327 789. The values of the hydraulic parameters are shown in Table 1.

Table 1 Default values of variables and hydraulic parameters.

The head solution to the nonlinear model of Eqs. (43)–(49) is obtained by an iteration method. Initially, the quadratic terms in Eq. (43) are assumed as $K_x{g}_x^{(n-\mathrm{1})}{G}_x^{\left(n\right)}+K_y{g}_y^{(n-\mathrm{1})}{G}_y^{\left(n\right)}+K_z{g}_z^{(n-\mathrm{1})}{G}_z^{\left(n\right)}$ with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{G}_x^{\left(n\right)}=\mathrm{\Delta }{x}_{\mathrm{w}}{\mathit{\varphi }}_{i+\mathrm{1},j,k}^{m+\mathrm{1},\left(n\right)}-\left(\mathrm{\Delta }{x}_{\mathrm{w}}+\mathrm{\Delta }{x}_{\mathrm{e}}\right){\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1},\left(n\right)}\\ \text{(50a)}& {\displaystyle }& {\displaystyle }+\mathrm{\Delta }{x}_{\mathrm{e}}{\mathit{\varphi }}_{i-\mathrm{1},j,k}^{m+\mathrm{1},\left(n\right)},{\displaystyle }& {\displaystyle }{G}_y^{\left(n\right)}=\mathrm{\Delta }{y}_{\mathrm{s}}{\mathit{\varphi }}_{i,j+\mathrm{1},k}^{m+\mathrm{1},\left(n\right)}-(\mathrm{\Delta }{y}_{\mathrm{s}}+\mathrm{\Delta }{y}_{\mathrm{n}}){\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1},\left(n\right)}\\ \text{(50b)}& {\displaystyle }& {\displaystyle }+\mathrm{\Delta }{y}_{\mathrm{n}}{\mathit{\varphi }}_{i,j-\mathrm{1},k}^{m+\mathrm{1},\left(n\right)},\text{(50c)}& {\displaystyle }& {\displaystyle }{G}_z^{\left(n\right)}=({\mathit{\varphi }}_{i,j,k+\mathrm{1}}^{m+\mathrm{1},\left(n\right)}-{\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1},\left(n\right)})/\mathrm{\Delta }z,{\displaystyle }& {\displaystyle }{g}_x^{(n-\mathrm{1})}=\mathrm{\Delta }{x}_{\mathrm{w}}{\mathit{\varphi }}_{i+\mathrm{1},j,k}^{m+\mathrm{1},(n-\mathrm{1})}-\left(\mathrm{\Delta }{x}_{\mathrm{w}}+\mathrm{\Delta }{x}_{\mathrm{e}}\right){\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1},(n-\mathrm{1})}\\ \text{(50d)}& {\displaystyle }& {\displaystyle }+\mathrm{\Delta }{x}_{\mathrm{e}}{\mathit{\varphi }}_{i-\mathrm{1},j,k}^{m+\mathrm{1},(n-\mathrm{1})},{\displaystyle }& {\displaystyle }{g}_y^{(n-\mathrm{1})}=\mathrm{\Delta }{y}_{\mathrm{s}}{\mathit{\varphi }}_{i,j+\mathrm{1},k}^{m+\mathrm{1},(n-\mathrm{1})}-(\mathrm{\Delta }{y}_{\mathrm{s}}+\mathrm{\Delta }{y}_{\mathrm{n}}){\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1},(n-\mathrm{1})}\\ \text{(50e)}& {\displaystyle }& {\displaystyle }+\mathrm{\Delta }{y}_{\mathrm{n}}{\mathit{\varphi }}_{i,j-\mathrm{1},k}^{m+\mathrm{1},(n-\mathrm{1})},\text{(50f)}& {\displaystyle }& {\displaystyle }{g}_z^{(n-\mathrm{1})}=({\mathit{\varphi }}_{i,j,k+\mathrm{1}}^{m+\mathrm{1},(n-\mathrm{1})}-{\mathit{\varphi }}_{i,j,k}^{m+\mathrm{1},(n-\mathrm{1})})/\mathrm{\Delta }z,\end{array}$$

where superscript (n) represents the nth iteration and gradients $g_x^{(n-\mathrm{1})}$, $g_y^{(n-\mathrm{1})}$, and $g_z^{(n-\mathrm{1})}$ cause a linearized Eq. (43) because they are known head values from the previous iteration. At the first time step (i.e., t=Δt, m=2), the first iteration solves a system of Eqs. (44)–(49) and the linearized Eq. (43) with $g_x^{\left(\mathrm{0}\right)}=g_y^{\left(\mathrm{0}\right)}=g_z^{\left(\mathrm{0}\right)}=\mathrm{1}$ and obtains the numerical solution of ${\mathit{\varphi }}_{i,j,k}^{\mathrm{2},\left(\mathrm{1}\right)}$ at each nodal point. The second iteration obtains ${\mathit{\varphi }}_{i,j,k}^{\mathrm{2},\left(\mathrm{2}\right)}$ with updated values of $g_x^{\left(\mathrm{1}\right)}$, $g_y^{\left(\mathrm{1}\right)}$, and $g_z^{\left(\mathrm{1}\right)}$ from the previous result of ${\mathit{\varphi }}_{i,j,k}^{\mathrm{2},\left(\mathrm{1}\right)}$. Repeat this iteration process for n≥3 until the convergence condition of $\left|{\mathit{\varphi }}_{i,j,k}^{\mathrm{2},\left(n\right)}-{\mathit{\varphi }}_{i,j,k}^{\mathrm{2},(n-\mathrm{1})}\right|<{\mathrm{10}}^{-\mathrm{4}}$ at each nodal point in the unsaturated zone is satisfied. The last result of ${\mathit{\varphi }}_{i,j,k}^{\mathrm{2},\left(n\right)}$ is therefore the head solution to the nonlinear model. Similarly, the iteration process is applied to obtain ${\mathit{\varphi }}_{i,j,k}^{m,\left(n\right)}$ for m≥3 at the other time steps (i.e., t=2Δt, 3Δt,…) with the convergence condition $\left|{\mathit{\varphi }}_{i,j,k}^{m,\left(n\right)}-{\mathit{\varphi }}_{i,j,k}^{m,(n-\mathrm{1})}\right|<{\mathrm{10}}^{-\mathrm{4}}$. Note that the first iteration at each time step calculates $g_x^{\left(\mathrm{0}\right)}$, $g_y^{\left(\mathrm{0}\right)}$, and $g_z^{\left(\mathrm{0}\right)}$ using ${\mathit{\varphi }}_{i,j,k}^{m,\left(n\right)}$ obtained at the previous time step.

3 Results and discussion

The default values of the parameters and variables used in the calculation of the present solution are listed in Table 1. In Sect. 3.1, the error arising from neglecting the process of infiltration in the unsaturated zone is examined. In Sect. 3.2, the recharge efficiency associated with the properties of the unsaturated zone is investigated. In Sect. 3.3, the sensitivity analysis of the hydraulic head in the unsaturated zone in regard to various hydraulic parameters is discussed. In Sect. 3.4, the present solution is compared with a finite-difference solution. In Sect. 3.5, the present solution is applied to a field problem of artificial recharge.

3.1 Effect of unsaturated flow on head distributions in aquifers

Figure 3 Temporal distributions of the dimensionless head in the saturated zone predicted by the present solution and Chang et al.'s (2016) solution for different pairs of (α, $\stackrel{\mathrm{‾}}{b}$) representing the effect of unsaturated flow.

Here we investigate the difference between the present solution and Chang et al.'s (2016) analytical solution to explore the effect of unsaturated flow on the head distributions in the aquifer. Chang et al.'s (2016) solution considers 3-D saturated flow in an unconfined aquifer with localized recharge but neglects the effect of unsaturated flow. One might expect that the difference will mainly be dominated by the magnitudes of parameters α (dimensionless unsaturated exponent) and $\stackrel{\mathrm{‾}}{b}$ (dimensionless unsaturated thickness). Figure 3 displays the predicted temporal head distributions at ($\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},\stackrel{\mathrm{‾}}{z})=(\mathrm{2},\mathrm{0},-\mathrm{0.5})$ by their solution and the present solution, Eq. (36b), for different pairs of (α, $\stackrel{\mathrm{‾}}{b}$) with α=10² or 10³ and $\stackrel{\mathrm{‾}}{b}$ from 10⁻³ to one. Significant difference in $\stackrel{\mathrm{‾}}{h}$ predicted by both solutions can be seen except for the cases (α, $\stackrel{\mathrm{‾}}{b})=({\mathrm{10}}^{\mathrm{2}},{\mathrm{10}}^{-\mathrm{1}})$, (10³, 10⁻¹), and (10³, 10⁻²) shown in the figure. The result indicates that the thickness of the unsaturated zone is less than 10 % of the saturated aquifer thickness (i.e., $\stackrel{\mathrm{‾}}{b}\le \mathrm{0.1}$) for obtaining close predictions from both solutions. When $\stackrel{\mathrm{‾}}{b}={\mathrm{10}}^{-\mathrm{2}}$, both solutions disagree if α=10² and agree well if α=10³, indicating that the magnitude of the product $\mathit{\alpha }\stackrel{\mathrm{‾}}{b}$ (=ab) should at least be 10 (i.e., $\mathit{\alpha }\stackrel{\mathrm{‾}}{b}\ge \mathrm{10}$) for good agreement of both solutions. It is worth noting that both solutions disagree even for a much thinner unsaturated zone as compared with the aquifer (i.e., $\stackrel{\mathrm{‾}}{b}={\mathrm{10}}^{-\mathrm{3}}$) because of $\mathit{\alpha }\stackrel{\mathrm{‾}}{b}<\mathrm{10}$. Additionally, the curves for α=10 all disagree with Chang et al.'s (2016) solution, whereas the curves for α=10⁴ match with this solution except for two cases where (α, $\stackrel{\mathrm{‾}}{b})=({\mathrm{10}}^{\mathrm{4}},\mathrm{0.5})$ and (10⁴, 1) (not shown in the figure). Judging from the above, one can recognize that the effect of unsaturated flow on the predicted head in saturated aquifers is negligible when $\mathit{\alpha }\stackrel{\mathrm{‾}}{b}\ge \mathrm{10}$ and $\stackrel{\mathrm{‾}}{b}\le \mathrm{0.1}$. A great number of existing analytical solutions ignoring unsaturated flow give accurate predictions only when those two inequalities are satisfied (e.g., Chang and Yeh, 2007; Illas et al., 2008; Bansal and Teloglou, 2013). Otherwise, significant deviations will happen in their predictions.

3.2 Effect of unsaturated flow on recharge efficiency

Figure 4 Temporal distribution curves of the recharge efficiency (RE) for different pairs of (α, $\stackrel{\mathrm{‾}}{b}$) representing the effect of unsaturated flow.

The effect of the unsaturated zone on the RE is explored based on the curves of the RE versus $\stackrel{\mathrm{‾}}{t}$ shown in Fig. 4 plotted using Eq. (40) for different pairs of $(\mathit{\alpha },\stackrel{\mathrm{‾}}{b})=(\mathrm{100},\mathrm{0.5})$, (10, 0.5), (1, 0.01), (1, 0.5), and (1, 1). For a given $\stackrel{\mathrm{‾}}{t}$, the RE increases with α for a fixed $\stackrel{\mathrm{‾}}{b}$ and decreasing $\stackrel{\mathrm{‾}}{b}$ for a fixed α. After $\stackrel{\mathrm{‾}}{t}={\mathrm{10}}^{\mathrm{6}}$, the RE approaches an ultimate value equalling unity when $(\mathit{\alpha },\stackrel{\mathrm{‾}}{b})=(\mathrm{100},\mathrm{0.5})$ and (1, 0.01), 0.9 when $(\mathit{\alpha },\stackrel{\mathrm{‾}}{b})=(\mathrm{10},\mathrm{0.5})$, 0.7 when $(\mathit{\alpha },\stackrel{\mathrm{‾}}{b})=(\mathrm{1},\mathrm{0.5})$, and 0.6 when $(\mathit{\alpha },\stackrel{\mathrm{‾}}{b})=(\mathrm{1},\mathrm{1})$. Those results imply that the ultimate recharge efficiency (URE) depends on the magnitudes of both α and $\stackrel{\mathrm{‾}}{b}$. Figure 5 illustrates contours of the URE at $\stackrel{\mathrm{‾}}{t}={\mathrm{10}}^{\mathrm{7}}$ in the ranges of $\mathrm{0.01}\le \stackrel{\mathrm{‾}}{b}\le \mathrm{1}$ and $\mathrm{1}\le \mathit{\alpha }\le \mathrm{100}$. The URE approaches unity when $\stackrel{\mathrm{‾}}{b}<\mathrm{0.05}$ or α>20. In contrast, it is below 0.9 and related to a given pair (α, $\stackrel{\mathrm{‾}}{b}$) when $\stackrel{\mathrm{‾}}{b}>\mathrm{0.1}$ and α<10. It is clearly seen that the RE is great for a large α and/or a small $\stackrel{\mathrm{‾}}{b}$. Those results provide useful information in the estimation of the amount of water from the recharge entering the aquifer. Notice that the case of URE <1 may be due to the problem that unsaturated flow is influenced by the water retention capacity and diffusivity in the horizontal direction.

3.3 Sensitivity analysis for flow in unsaturated zone

Chang et al. (2016) performed the sensitivity analysis to investigate the sensitivity of the hydraulic head in saturated aquifers to the change in each of the aquifer parameters. This section focuses on the sensitivity analysis of the head in the unsaturated zone. Consider the recharge area of $\mathrm{0}\le x\le \mathrm{50}$ m and $\mathrm{0}\le y\le \mathrm{50}$ m and the observation points A at (0, 0, 5 m) under the area and B at (100 m, 0, 5 m) beside the area. Other values of the parameters are given in Table 1. The temporal distribution curves of the normalized sensitivity coefficient S_i,t predicted by Eq. (42) to each of the parameters a, l, w, S_s, S_y, K_x, K_y, and K_z are exhibited in Fig. 6a for point A and Fig. 6b for point B. At a given time, a positive S_i,t means that the change in the specific parameter causes the increase in the head. In contrast, a negative S_i,t signifies that the change leads to the head decrease. The magnitude of the head remains unchanged when $S_{i,t}=\mathrm{0}$. Obviously, the parameters l, w, S_y, K_x, and K_y are important factors affecting the predicted head observed at points A and B, revealing that those parameters should be included in the flow model. The head at point A is sensitive to the changes in a and K_z, but that at point B is insensitive. The result implies that unsaturated flow prevails under the recharge area but does not away from the area. In addition, the coefficient S_i,t to S_s almost equals zero over the entire recharge period, indicating that the change in S_s does not affect the predicted head in the unsaturated zone.

Figure 5 Contours of the ultimate recharge efficiency (URE) plotted at $\stackrel{\mathrm{‾}}{t}={\mathrm{10}}^{\mathrm{7}}$ for various values of α and $\stackrel{\mathrm{‾}}{b}$.

Figure 6 Temporal distribution curves of the normalized sensitivity coefficients for the hydraulic head in the unsaturated zone in response to the change in each of the parameters a, l, w, S_s, S_y, K_x, K_y, and K_z observed at (a) (0, 0, 5 m) under the recharge area and (b) (100 m, 0, 5 m) beside the area.

3.4 Validation of the present solution

The finite-difference solution to the unsaturated–saturated flow model based on the nonlinear and linearized versions of Richards' equation, Eqs. (5) and (6), has been developed and described in Sect. 2.6. It is used to validate the present solution. Figure 7 demonstrates temporal head distributions observed at (25 m, 0, 5 m) and (25 m, 0, −10 m) under the recharge area and at (94.65 m, 0, 5 m) and (94.65 m, 0, −10 m) beside the area. The figure displays good agreement on the predicted head distributions from both the solutions. It is noteworthy that numerous attempts had been made by scholars to examine the accuracy of the linearized version of Richards' equation (e.g., Kroszynski and Dagan, 1975; Mishra and Neuman, 2010; Liang et al., 2017b). They also revealed that the linearized equation causes insignificant deviation in model predictions. We therefore conclude that the present solution is correctly developed and fairly predicts the hydraulic head for the unsaturated–saturated flow induced by localized recharge.

3.5 Application of the present solution to the field experiment

Figure 7 Temporal distributions of the hydraulic head predicted by the present solution and the finite-difference solution based on Richards' equation, Eq. (5), and its linearized version, Eq. (6), observed at (25 m, 0, 5 m) and (25 m, 0, −10 m) under the recharge area and at (94.65 m, 0, 5 m) and (94.65 m, 0, −10 m) beside the area.

Bianchi and Haskell (1968) executed a field experiment of artificial recharge from two ponds on an alluvial fan in Fresno County, California. Pond No. 2 was within a square of 90 m × 90 m on the ground. The average recharge rate of the pond was 0.107 m d⁻¹. The initial water table was 6.4 m below the ground and 24.384 m over the impervious aquifer bottom. The entire recharge period was 10.92 days on record. The values of the aquifer parameters obtained from well test data were K_x=7.925 m d⁻¹ and S_y=0.022. There are 19 observation data of the water table rise beneath the center of the pond versus time shown in Fig. 8. We apply the least square method using the Mathematica FindRoot built-in function (Wolfram, 1996) to estimate five parameters a, K_x=K_y, K_z, S_s, and S_y based on the data and the present solution. The estimated values are a=0.388 m⁻¹, K_x=5.642 m d⁻¹, K_z=1.573 m d⁻¹, $S_{\mathrm{s}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ m⁻¹, and S_y=0.102, which are all in the reasonable ranges of their parameter values; they are $\mathrm{0.2}\le a\le \mathrm{5}$ m⁻¹ (Philip, 1969), $\mathrm{8.64}\times {\mathrm{10}}^{-\mathrm{2}}\le K_x\le \mathrm{864}$ m d⁻¹, $\mathrm{0.1}{K}_x\le K_z\le \mathrm{0.33}{K}_x$, ${\mathrm{10}}^{-\mathrm{5}}\le S_{\mathrm{s}}B\le {\mathrm{10}}^{-\mathrm{3}}$, and $\mathrm{0.01}\le S_y\le \mathrm{0}$.3 for sandy aquifers (Freeze and Cherry, 1979, p. 604). Figure 8 demonstrates 19 observed data of the water table rise, the predictions from Glover's (1960) solution with K_x=7.925 m d⁻¹ and S_y=0.022 provided in Bianchi and Haskell (1968), and the present solution with the five estimated parameters. The Glover solution was developed by assuming that the flow is radially outward from a circular recharge pond with an equivalent area to the square of 90 m × 90 m and the unsaturated flow above water table is neglected. The predictions from the present solution agree well with the observed data, but those from Glover's solution do not, indicating that the effect of unsaturated flow had better be considered because $\mathit{\alpha }\stackrel{\mathrm{‾}}{b}=\mathrm{2.48}$ and $\stackrel{\mathrm{‾}}{b}=\mathrm{0.26}$ in this case do not satisfy the condition of $\mathit{\alpha }\stackrel{\mathrm{‾}}{b}\ge \mathrm{10}$ and $\stackrel{\mathrm{‾}}{b}\le \mathrm{0.1}$ concluded in Sect. 3.1. From those discussed above, the present solution has been shown to be applicable to a real-world problem for unsaturated–saturated flow due to a recharge pond.

Figure 8 Comparison of the water table rise predicted by the present solution and Glover's (1960) solution with field-observed data given in Bianchi and Haskell (1968).

4 Concluding remarks

This study develops a novel mathematical model depicting 3-D unsaturated–saturated flow for the process that surface water recharge passes through an unsaturated zone and flows down to an unconfined aquifer. The Richards equation is considered to delineate unsaturated flow induced by infiltration due to recharge from the ground surface. The Gardner model is used to describe the unsaturated soil characteristics. The transient groundwater flow equation is then employed to describe the rise of the hydraulic head in the aquifer in response to the water flow from the unsaturated zone. Both equations are coupled by the continuity equations of the head and flux at the water table. The head solution to the model is derived by means of the Laplace transform and Fourier cosine transform. The recharge efficiency defined as the percentage of the water from the recharge down to the aquifer is clearly discussed. The sensitivity analysis is performed to investigate the head response to the change in each of the hydraulic parameters in the unsaturated zone. The present solution agrees well with the finite-difference solution in predicting the time-varying head for the unsaturated–saturated flow model. In addition, the present solution is applied to study the observed data from a field experiment conducted by Bianchi and Haskell (1968). On the basis of the studies obtained from the present solution, the following conclusions can be drawn.

The effect of unsaturated flow on the hydraulic head in the aquifer is ignorable when the product of the unsaturated exponent (a) and initial unsaturated thickness (b) is greater than 10 (i.e., ab≥10) and the unsaturated thickness is less than 10 % of the initial aquifer thickness (B) (i.e., $b/B\le \mathrm{0.1}$). Otherwise, the effect should be considered to avoid large deviations in calculating the head in the aquifer. Existing models considering only saturated flow can predict accurate results only when these two inequalities are satisfied.

The recharge efficiency initially equals zero, increases with time, and finally approaches a constant value (below or equal to unity) depending on the values of α(=aB) and $\stackrel{\mathrm{‾}}{b}(=b/B)$.

The ultimate recharge efficiency approaches unity when $\stackrel{\mathrm{‾}}{b}<\mathrm{0.05}$ or α>20 but less than 90 % when $\stackrel{\mathrm{‾}}{b}>\mathrm{0.1}$ and α<10. In other words, the surface source supplies more recharge water to the aquifer if the unsaturated zone has a large α and/or a small $\stackrel{\mathrm{‾}}{b}$.

The results of the sensitivity analysis indicate that the parameter a, l, or w causes positive influence but S_y, K_x, K_y, or K_z produces negative impact on the predicted head in the unsaturated zone. The head under the recharge area is sensitive to the changes in a and K_z, but that beside the area is not. Moreover, the head is rather insensitive to the change in S_s.