Technical Note: Three-dimensional transient groundwater flow due to localized recharge with an arbitrary transient rate in unconfined aquifers

Most previous solutions for groundwater flow induced by localized recharge assumed either aquifer incompressibility or two-dimensional flow in the absence of the vertical flow. This paper develops a new three-dimensional flow model for hydraulic head variation due to localized recharge in a rectangular unconfined aquifer with four boundaries under the Robin condition. A governing equation describing spatiotemporal head distributions is employed. The first-order free-surface equation with a source term defining a constant recharge rate over a rectangular area is used to depict water table movement. The solution to the model for the head is developed with the methods of Laplace transform and double-integral transform. Based on Duhamel’s theorem, the present solution is applicable to flow problems accounting for arbitrary time-dependent recharge rates. The solution to depth-average head can then be obtained by integrating the head solution to elevation and dividing the result by the aquifer thickness. The use of a rectangular aquifer domain has two merits. One is that the integration for estimating the depth-average head can be analytically achieved. The other is that existing solutions based on aquifers of infinite extent can be considered as special cases of the present solution before the time when the aquifer boundary had an effect on head predictions. With the help of the present solution, the assumption of neglecting the vertical flow effect on the temporal head distribution at an observation point outside a recharge region can be assessed by a dimensionless parameter related to the aquifer horizontal and vertical hydraulic conductivities, initial aquifer thickness, and the shortest distance between the observation point and the edge of the recharge region. The validity of assuming aquifer incompressibility is dominated by the ratio of the aquifer specific yield to its storage coefficient. In addition, a sensitivity analysis is performed to investigate the head response to the change in each of the aquifer parameters.


Introduction
The water table rises due to localized recharge, such as rainfall, lakes, and agricultural irrigation, into the regional area of the aquifer.Excess recharge may cause soil liquefaction or wet basements of buildings.Groundwater flow behavior induced by recharge is therefore crucial in water resource management.The Boussinesq equation has been extensively used to describe horizontal flow without the vertical component in unconfined aquifers (e.g., Ireson and Butler, 2013;van der Spek et al., 2013;Yeh and Chang, 2013;Chor and Dias, 2015;Hsieh et al., 2015;Liang and Zhang, 2015;Liang et al., 2015).The equation can be linearized by assuming uniform saturated aquifer thickness for developing its analytical solution.Marino (1967) presented quantitative criteria for the validity of the linearized Boussinesq equation.The criteria are introduced in the next section.
The rate of localized recharge can be a constant for the long term but should be dependent of time for the short term (Rai et al., 2006).An exponentially decaying function of time is usually used for recharge intensity decreasing from a certain rate to an ultimate one.An arbitrary time-dependent recharge rate is commonly approximated as the combination of several linear segments of time to develop analytical solutions for water table rise subject to the recharge.Analytical models accounting for water table rise due to the recharge region of an infinite-length strip are reviewed.One-dimensional (1-D) flow perpendicular to the strip is considered while the flow along the strip is assumed ignorable.These models deal with aquifers of infinite or finite extent with various types of outer boundary conditions.Hantush (1963) considered an aquifer of infinite extent without a lateral boundary.Rao and Sarma (1980) considered an aquifer of finite extent with two constant-head (also called Dirichlet) boundaries.Later, they developed a solution (Rao and Sarma, 1984) for a finite-extent aquifer between no-flow and constant-head boundaries.Latinopoulos (1986) deliberated on a finite-extent aquifer between two boundaries, one of which is under the Robin condition and the other is under either the Dirichlet or no-flow condition.The recharge rate is treated as a periodical pulse consisting of constant rates for rainy seasons and zero for dry seasons.Bansal and Das (2010) studied an aquifer extending semi-infinitely from a Dirichlet boundary and overlying a sloping impervious base, and indicated that the change in groundwater mound induced by a strip-shaped recharge region increases with the base slope.
A variety of analytical models were presented to describe water table rise for 2-D flow induced by rectangle-shaped recharge into unconfined aquifers.The differences between these solutions are addressed below.Hantush (1967) considered an infinite-extent aquifer with localized recharge having a constant rate.Manglik et al. (1997) handled an arbitrary time-varying rate of recharge into a rectangular aquifer bounded by no-flow stratum.Manglik and Rai (1998) investigated flow behavior based on an irregularly time-varying rate of recharge into a rectangular aquifer with the lateral boundary under the Dirichlet condition.Bruggeman (1999) introduced an analytical solution for steady-state flow induced by localized recharge into a vertical strip aquifer between two Robin boundaries.Chang and Yeh (2007) considered one localized recharge and multiple extraction wells in an anisotropic aquifer overlying an impervious sloping bed.They indicated that the aquifer anisotropy and bottom slope notably influence water table distributions.Bansal and Teloglou (2013) explored the problem of a groundwater mound subject to multiple localized recharges and withdrawal wells in an unconfined aquifer overlying a semi-permeable base.They indicated that groundwater mound rises as the aquifer hydraulic conductivity decreases.
Some articles discussed water table rise near circle-shaped recharge regions and thus considered radial groundwater flow, which is symmetric to the center of the region.Rai et al. (1998) presented an analytical model describing water table growth subject to an exponentially decaying rate of recharge in a circle-shaped unconfined aquifer with an outer Dirichlet boundary.Illas et al. (2008) considered the same model but a leaky aquifer.They indicated that leakage across the aquifer bottom significantly influences spatiotemporal water table distributions despite a small amount of the leakage.On the other hand, some researchers considered radial flow having the vertical component near a circle-shaped recharge region of an infinite-extent unconfined aquifer.A first-order free-surface equation as the top boundary condition of the aquifer is applied to describe water table rise.Zlotnik and Ledder (1992) developed analytical models for describing the distributions of hydraulic head and flow velocity due to constant-rate recharge.They found that models neglecting aquifer compressibility overestimate the magnitudes of the head and flow velocity.Ostendorf et al. (2007) derived an analytical model for head variation due to an exponentially decaying rate of recharge.Predictions of their solution agreed well with the field data obtained in the Plymouth-Carver aquifer in southeastern Massachusetts given by Hansen and Lapham (1992).
Some studies developed a 3-D flow model based on the Laplace equation, which neglects the aquifer compressibility effect.Dagan (1967) derived an analytical solution for the velocity potential caused by regional recharge into an unconfined aquifer of infinite thickness.Zlotnik and Ledder (1993) also developed an analytical solution for the same model but considered finite thickness for the unconfined aquifer.Predictions of their solutions indicate that groundwater flow is horizontal in the area beyond 150 % of the length or width of a rectangular recharge region.
It would be informative to summarize the abovementioned models in Table 1.The solutions to the models are classified according to flow dimensions into 1-D, 2-D, 3-D, and radial flows and further categorized according to aquifer domain, aquifer boundary conditions, recharge region, and recharge rate.The table shows that those solutions assume neither vertical flow nor aquifer incompressibility.In addition, the Dirichlet and no-flow conditions considered by some of those solutions are not applicable to a boundary having a semipermeable stratum, but the Robin condition is.The former two conditions are indeed special cases of the third one.
The objective of this paper is to develop a new mathematical model for depicting spatiotemporal hydraulic head distributions subject to localized recharge with an arbitrary time-varying recharge rate in a rectangular-shaped unconfined aquifer.The four boundaries are considered under the Robin condition, which can reduce to the Dirichlet or no-flow condition.A governing equation describing 3-D transient flow subject to the effect of aquifer compressibility is used.A first-order free-surface equation with a source term representing recharge rate is chosen to describe the top boundary condition.The transient head solution of the model is derived by the methods of Laplace transform, double-integral transform, and Duhamel's theorem.The sensitivity analysis based on the present solution is performed to study the head response to the change in each of the hydraulic parameters.On the basis of the solution's predictions, the effect of the Robin boundaries on time-dependent head distributions at observation points is investigated.A quantitative criterion under which the Robin condition reduces to the Dirichlet or no-flow  condition is provided.In addition, quantitative criteria for the validity of two assumptions of aquifer incompressibility and no vertical flow are provided, and errors arising from the assumptions in the hydraulic head are also discussed.Temporal head distributions accounting for transient recharge rates are demonstrated as well.

Mathematical model
A mathematical model is developed for describing spatiotemporal hydraulic head distributions induced by localized recharge in a rectangular unconfined aquifer as illustrated in Fig. 1a.The four boundaries of the aquifer are considered under the Robin condition.The aquifer has the widths of l and w in x and y directions, respectively.The recharge uniformly distributes over a rectangular region having widths a and b in x and y directions, respectively.The lower left corner of the region is designated at (x 1 , y 1 ).The shortest distances measured from the edge of the region to boundaries 1, 2, 3, and 4 are denoted as d 1 , d 2 , d 3 , and d 4 , respectively.The shortest distance between the edge of the region and an observation point at (x, y) is defined as , where (x e , y e ) is a coordinate on the edge.The initial aquifer's thickness is B as shown in Fig. 1b.The governing equation describing 3-D transient head distributions in a homogeneous and anisotropic aquifer is expressed as where t is time, h(x, y, z, t) represents the hydraulic head, K x , K y , and K z are the hydraulic conductivities in x, y, and z directions, respectively, and S s is the specific storage.The initial static water table is chosen as the reference datum where the elevation head is set to zero.The initial condition is therefore written as The Robin conditions specified at the four sides of the aquifer are defined as where subscripts 1, 2, 3, and 4 represent the boundaries at x = 0, x = l, y = 0, and y = w, respectively, and K and b are the hydraulic conductivity and width of the medium at the aquifer boundary, respectively.Note that Eqs.(3)-( 6) reduce to the Dirichlet condition when b (i.e., b 1 , b 2 , b 3 , or b 4 ) is set to zero and the no-flow condition when K (i.e., K 1 , K 2 , K 3 , or K 4 ) is set to zero.The aquifer lies on an impermeable base denoted as The first-order free-surface equation describing the response of water table to recharge over the rectangular region can be written as (Zlotnik and Ledder, 1993) where S y is the specific yield, I is a recharge rate, and u( ) is the unit step function.Equation (8a) involves the assumption of I K z and the simplification from non-uniform saturated aquifer thickness below z = h to a uniform thickness Hydrol.Earth Syst.Sci., 20, 1225-1239, 2016 www.hydrol-earth-syst-sci.net/20/1225/2016/ below z = 0 (Dagan, 1967).Marino (1967) indicated that the simplification and assumption are valid when the water table rise is smaller than 50 % of the initial water table height (i.e., |h|/B < 0.5) and the recharge rate is smaller than 20 % of the hydraulic conductivity (i.e., I /K z < 0.2).
On the other hand, the effect of unsaturated flow above the water table on the model's predictions can be ignored when σ B ≥ 10 3 , where σ is a parameter to define the relative hydraulic conductivity as k 0 = exp (−σ z) in the Richards' equation (Tartakovsky and Neuman, 2007).Tartakovsky and Neuman (2007) achieved agreement on aquifer drawdown, which was evaluated by their analytical solution based on Eq. ( 1) for saturated flow and Richards' equation for unsaturated flow, and by the Neuman (1974) solution based on Eqs. ( 1) and ( 8a) with I = 0 when σ B = 10 3 (i.e., the case of κ D = 10 3 in Fig. 2 in Tartakovsky and Neuman, 2007).Dimensionless variables and parameters are defined as follows where the overbar denotes a dimensionless symbol.Notice that the variables in the horizontal and vertical directions are divided by d and B, respectively.According to Eq. ( 9), the mathematical model, Eqs. ( 1)-(8c), can then be expressed as

Analytical solution
The mathematical model, Eqs. ( 10)-(17c), can be solved by the methods of Laplace transform and double-integral transform.The former transform converts h(x, y, z, t) into h(x, y, z, p), ∂h/∂t into p h − h| t=0 , and ξ u x u y into ξ u x u y /p, where p is the Laplace parameter and h| t=0 equals zero in Eq. ( 11).After taking the transform, the model becomes a boundary value problem expressed as with boundary conditions and ∂ h/∂z + εp h/κ z = ξ u x u y /p at z = 0. We then apply the properties of the double-integral transform to the problem.
One can refer to the definition in Latinopoulos (1985, Table I, aquifer type 1).The transform turns h(x, y, z, p) into ĥ(α m , In addition, u x u y defined in Eqs. ( 17b) and ( 17c) is transformed into U m U n given by with where χ = x 1 + a and ψ = y 1 + b.Equation ( 18) then reduces to an ordinary differential equation as Two boundary conditions are expressed, respectively, as and Solving Eq. ( 25) with Eqs. ( 26) and ( 27) results in where Inverting Eq. ( 28) to the space and time domains gives rise to the following analytical solution: φ 0,m,n = −2λ 0,m,n cosh (1 + z)λ 0,m,n exp −γ 0,m,n t /η 0,m,n , (30c) φ j,m,n = −2λ j,m,n cos (1 + z)λ j,m,n exp −γ j,m,n t /η j,m,n , (30d) η j,m,n = γ j,m,n (1 + 2εκ z) λ j,m,n cos λ j,m,n + 1 − εγ j,m,n sin λ j,m,n , (30f) where j ∈ 1, 2, 3, . . .∞ and eigenvalues λ 0,m,n and λ j,m,n are determined, respectively, by the following equations: and Notice that Eqs. ( 19), (20), and (31) have infinite positive roots owing to the trigonometric function tan( ) while Eq. ( 32) has only one positive root.The method to find α m , β n , λ j,m,n , and λ 0,m,n is introduced in Sect.2.3.One can refer to Appendix A for the derivation of Eq. (30a).The first term on the right-hand side (RHS) of Eq. ( 30a) is a double series expanded by α m and β n .The series converges within a few terms because the power of α m (or β n ) in the denominator of φ m,n in Eq. ( 30b) is two more than that in the nominator.The second term on the RHS of Eq. ( 30a) is a double series expanded by α m and β n , and the third term is a triple series expanded by α m , β n , and λ j,m,n .They converge very fast due to exponential functions in Eqs.(30c) and (30d).Consider (m, n) ∈ (1, 2, . . ., N = 30) and j ∈ (1, 2, . . ., N j = 15) for the default values of dimensionless parameters and variables in Table A1 for calculation.The number of terms in one or the other double series is 30 × 30 = 900 and in the triple series is 30 × 30 × 15 = 13 500.The total number is therefore 900 × 2 + 13 500 = 15 300.We apply Mathematica FindRoot routine to obtain the values of α m , β n , and λ j,m,n and Sum routine to compute the double and triple series.It takes about 8 s to finish calculation for t = 10 5 with a personal computer with Intel Core i5-4590 3.30 GHz processor and 8 GB RAM.In addition, the series is considered to converge when the absolute value of the last term in the double series of φ m,n is smaller than 10 −20 (i.e., 10 −50 < 10 −20 in this case).That value in the other double or triple series may be even smaller than 10 −50 due to exponential decay.The use of finite aquifer domain has two merits.One is that the solution to depth-average head, defined as 0 −1 h(x, y, z, t) dz, can be analytically integrated.The integration variable z appears only in the functions of cosh[(1 + z)λ m,n ] in Eq. (30b), cosh[(1 + z)λ 0,m,n ] in Eq. ( 30c) and cos[(1 + z)λ j,m,n ] in Eq. (30d).The solution to depthaverage head therefore equals Eq. (30a), where these three functions are replaced by sinh λ m,n /λ m,n , sinh λ 0,m,n /λ 0,m,n , and sin λ j,m,n /λ j,m,n , respectively.The other merit is that the present solution is applicable to head predictions in aquifers of infinite extent before the dimensionless time to have lateral aquifer boundary effect on the predictions.Wang and Yeh (2008) reported a time criterion defined as t cr = 0.03(1 + ε)R 2 , where R = R/d denotes a shortest dimensionless distance from the lateral boundary to the edge of the recharge region.This criterion is, in effect, a boundaryeffect time when the hydraulic head is affected by the aquifer boundary.Existing solutions based on aquifers of infinite extent can therefore be considered as special cases of the present solution before the boundary-effect time.

Solution for time-varying recharge rate
The present solution, Eq. (30a), is applicable to arbitrary time-dependent recharge rates on the basis of Duhamel's theorem expressed as (e.g., Bear, 1979, p. 158) where h I t signifies a dimensionless head solution for a timedependent recharge rate ξ t (τ ) with t replaced by τ , h I 0 is Eq.(30a) in which ξ is replaced by ξ t (0), and h(t − τ ) is also Eq. (30a) with t replaced by t − τ .If Eq. ( 33) is not integrable, it can be discretized as (Singh, 2005) with where h N represents a numerical result of dimensionless head h at t = t × N, t is a dimensionless time step, ξ i and ξ i−1 are dimensionless recharge rates at t = t × i and t = t × (i − 1), respectively, and η(M), called ramp kernel, depends on Eq. (30a) in which t is replaced by M t − τ .The integration result of Eq. ( 34c) can be denoted as Eq.(30a) where φ m,n is replaced by φ m,n t and two exponential terms in Eqs.(30c) and (30d) are replaced, respectively, by exp(−M γ 0,m,n t)[−1 + exp(γ 0,m,n t)]/γ 0,m,n and exp(−M γ j,m,n t)[−1 + exp(γ j,m,n t)]/γ j,m,n .

Sensitivity analysis
The sensitivity analysis is administered to assess the change in the hydraulic head in response to the change in each of the hydraulic parameters.The normalized sensitivity coefficient of the hydraulic head to a specific parameter can be expressed as where P c is the cth parameter in the present solution, S c,t is the coefficient at a time to the cth parameter, and h is the present solution, Eq. (30a).The derivative in Eq. ( 35) can be approximated as where P c is an increment chosen as 10 −3 P c (Yeh et al., 2008).

Results and discussion
Previous articles have discussed groundwater mounds in response to localized recharge into aquifers with various hydraulic parameters (e.g., Dagan, 1967;Rao and Sarma, 1980;Latinopoulos, 1986;Manglik et al., 1997;Manglik and Rai, 1998;Rai et al., 1998;Chang and Yeh, 2007;Illas et al., 2008;Bansal and Das, 2010;Bansal and Teloglou, 2013).Flow velocity fields below groundwater mounds have also been analyzed (Zlotnik andLedder, 1992, 1993).This section therefore focuses on the transient behavior of hydraulic head at an observation point with the aid of the present solution.The default values of the parameters and variables for calculation are noted in Table A1.In Sect.3.1, transient head distributions subject to Dirichlet, no-flow and Robin boundary conditions are compared.In Sect.3.2, the effect of vertical flow on the head distribution is investigated.In Sect.3.3, errors arising from assuming aquifer incompressibility (i.e., S s = 0) to develop analytical solutions are discussed.In Sect.3.4, the response of the hydraulic head to transient recharge rates based on Eq. ( 33) is demonstrated.In Sect.3.5, the sensitivity analysis defined by Eq. ( 36) is performed.

Effect of lateral boundary
The Robin condition can become the Dirichlet or no-flow condition, depending on the magnitudes of κ 1 d 1 for Eq. ( 12), κ 2 d 2 for Eq. ( 13), κ 3 d 3 for Eq. ( 14), and κ 4 d 4 for Eq. ( 15).We consider a symmetrical aquifer system with l = w = 22, The magnitudes of κ 1 d 1 , κ 2 d 2 , κ 3 d 3 , and κ 4 d 4 are the same and defined as κ.The curves of h versus t plotted by the present solution, Eq. (30a), for κ = 10 −3 , 10 −2 , 10 −1 , 1, 10, 100, and 200 are shown in Fig. 2. The curves h versus t are plotted from the Manglik et al. (1997) solution with the noflow condition (i.e., κ = 0), the Manglik and Rai (1998) solution with the Dirichlet condition (i.e., κ → ∞), and the present solution with the Robin condition.Before t = 10 4 , these curves give the same magnitude of h at a fixed dimensionless time t, since the lateral aquifer boundary has been beyond the place where groundwater is affected by localized recharge.After t = 10 4 , the curves for the cases of κ = 10 −2 , 10 −1 , 1, 10, and 100 deviate from each other gradually as time increases.A larger magnitude of κ between κ = 10 −2 and κ = 100 causes a smaller h at a fixed t.On the other hand, the present solution for the cases of κ = 10 −3 and 10 −2 agrees well with the Manglik et al. (1997) solution based on κ = 0 and that for the cases of κ = 100 and 200 predicts the same result as the Manglik and Rai (1998) solution based on κ → ∞.We may reasonably conclude that the Robin condition reduces to the no-flow condition when κ ≤ 10 −2 and the Dirichlet condition when κ ≥ 100.

Effect of vertical flow
Dimensionless parameter κ z (i.e., K z d 2 /(K x B 2 )) dominates the effect of vertical flow on transient head distributions at an observation point.Consider The vertical flow prevails, and its effect should be taken into account when κ z < 1, indicating a thick aquifer, a small ratio of K z /K x , and/or an observation point near a recharge region.On the other hand, the present solution for the cases of κ z = 1 and 10 agrees well with Manglik and Rai (1998) solution, indicating that the vertical flow effect is ignorable when κ z ≥ 1.We can recognize from the agreement that existing solutions neglecting the vertical flow effect give good predictions when κ z ≥ 1.

Effect of specific storage
Some of the existing models use the Laplace equation as a governing equation by assuming S s = 0 (e.g., Singh, 1976;Schmitz and Edenhofer, 1988;Zlotnik and Ledder, 1993).
The assumption is valid when ε (i.e., S y /(S s B)) is larger than a certain value.This section quantifies the value.The Zlotnik and Ledder (1993)

Sensitivity analysis
Consider point A of (555, 500, −10 m) at a 3-D flow region (i.e., κ z < 1) and point B of (800, 500, −10 m) at a 2-D flow region (i.e., κ z ≥ 1) as discussed in Sect.3.2.Localized recharge distributes over the square area of 450 m ≤ x ≤ 550 m and 450 m ≤ y ≤ 550 m.The distance d herein is set to 5 m for point A and 250 m for point B. The aquifer system is of isotropy with K x = K y and symmetry with K 1 = K 2 = K 3 = K 4 for conciseness.To investigate the responses of the hydraulic heads at these two points to the change in each of a, b, S s , S y , K x (or K y ), K z , and K 1 (or K 2 , K 3 , and K 4 ), the sensitivity analysis is performed by Eq. ( 36).The curves of the normalized sensitivity coefficient S c,t versus t for these seven parameters are shown in Fig. 6a for point A and A first-order free-surface equation with a source term representing the recharge is employed for describing the water table movement.The analytical head solution of the model is obtained by applying the Laplace transform, the doubleintegral transform, and Duhamel's theorem.The use of rectangular aquifer domain leads to two merits.One is that the integration for the solution to the depth-average head can be analytically done.The other is that existing solutions based on aquifers of infinite extent are special cases of the present solution when the recharge time is less than the boundaryeffect time.The present solution is applicable under the con-  8a) neglecting the effect of unsaturated flow above water table (Marino, 1967;Tartakovsky and Neuman, 2007).The sensitivity analysis is performed to explore the response of the head to the change in each of hydraulic parameters.With the aid of the present solution, the following conclusions can be drawn: 1.With respect to affecting h at observation points, the Robin condition specified at x = 0 reduces to the Dirichlet condition when κ 1 d 1 ≥ 100 (i.e., K 1 d 1 /(K x b 1 ) ≥ 100) and no-flow condition when κ 1 d 1 ≤ 10 −2 .The quantitative criteria for κ 1 d 1 are applicable to κ 2 d 2 , κ 3 d 3 , and κ 4 d 4 for the Robin conditions specified at x = l, y = 0, and y = w, respectively.
2. The vertical flow causes a significant decrease in the hydraulic head at an observation point when κ z < 1 (i.e., K z d 2 /(K x B 2 ) < 1).When κ z ≥ 1, the effect of vertical flow on the head is ignorable, and conventional models considering 2-D flow without the vertical component can therefore predict accurate results.
3. The 3-D Laplace equation based on the assumption of S s = 0 can be regarded as a flow governing equation when ε ≥ 100 (i.e., S y /(S s B) ≥ 100) for the whole aquifer domain.Otherwise, head predictions based on the Laplace equation are overestimated.
4. The abovementioned conclusions are also applicable to problems of groundwater flow subject to recharge with arbitrary time-varying rates.

C
.-H.Chang et al.: Technical Note: Three-dimensional transient groundwater flow due to localized recharge

Figure 1 .
Figure 1.Schematic diagram of a rectangular-shaped unconfined aquifer with localized recharge (a) top view (b) cross section view.

Figure 2 .
Figure 2. Temporal distributions of the dimensionless head predicted by the Manglik et al. (1997) solution for a no-flow boundary, the Manglik and Rai (1998) solution for a Dirichlet boundary, and the present solution with κ z = 1 for a Robin boundary.
and eigenvalues α m and β n are the positive roots of the following equations: www.hydrol-earth-syst-sci.net/20et al.: Technical Note: Three-dimensional transient groundwater flow due to localized recharge

Figure 3 .
Figure 3. Temporal distributions of the dimensionless head predicted by the Manglik and Rai (1998) solution based on 2-D flow and the present solution for 3-D flow with various κ z .
100 for lateral aquifer boundaries under the Dirichlet condition as discussed in Sect.3.1.The temporal distributions of h predicted by the present solution, Eq. (30a), with κ z = 0.01, 0.1, 1, and 10 are demonstrated in Fig. 3.The temporal distribution of h predicted by the Manglik and Rai (1998) solution based on 2-D flow without the vertical component is taken in order to address the effect of vertical flow.The figure reveals that h increases with κ z when κ z ≤ 1.The difference in h predicted by both solutions indicates the vertical flow effect.The Manglik and Rai (1998) solution obviously overestimates the head.

Figure 4 .
Figure 4. Temporal distributions of the dimensionless head for (a) κ z = 10 and (b) κ z = 10 −2 predicted by the Zlotnik and Ledder (1993) solution based on the assumption of S s = 0 and the present solution relaxing the assumption.
Fig. 6b for point B. The figure shows that the hydraulic heads at both points are more sensitive to the changes in a, b, K x , and S y than those in the others.This may indicate that a flow model should include at least these four parameters.The figure also shows that the heads at points A and B are insensitive to the change in K 1 because of κ 1 d 1 = 4500 > 100 as discussed in Sect.3.1.In addition, S c,t to K z for point A is nonzero after t = 0.4 day due to κ z = 6.25 × 10 −3 < 1 as discussed in Sect.3.2.In contrast, S c,t to K z for point B is very close to zero over the entire period because of κ z = 15.625 > 1.Moreover, the heads at points A and B are insensitive to the change in S s due to ε = 500 > 100 as discussed in Sect.3.3.4ConclusionsA mathematical model is developed to depict spatiotemporal head distributions induced by localized recharge with an arbitrary time-varying rate in a rectangular unconfined aquifer bounded by Robin boundaries with different hydraulic parameters.A governing equation for 3-D flow is considered.

Table 1 .
Classification of existing analytical solutions involving localized recharge.