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 1Schematic diagram of unsaturated–saturated flow in an unconfined aquifer system with localized recharge: (a) top view and (b) cross-sectional view.

Download

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

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)

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

and

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

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

The initial conditions for those two zones are

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

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

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

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

and

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:

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

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

with the property that

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

and

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

and

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

and

Similarly, the transformed boundary conditions are expressed as

and

The transformed continuity conditions are written as

and

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

and

with

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

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)

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 2Schematic diagram of finite-difference grids: (a) top view and (b) cross-sectional view.

Download

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

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

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

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

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

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

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

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

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

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 1Default values of variables and hydraulic parameters.

Download Print Version | Download XLSX

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

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.1 Effect of unsaturated flow on head distributions in aquifers

Figure 3Temporal 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.

Download

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 4Temporal distribution curves of the recharge efficiency (RE) for different pairs of (α, $\stackrel{\mathrm{‾}}{b}$) representing the effect of unsaturated flow.

Download

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