Boussinesq Approach

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

1

Environment Protection Program ADAM SZYMAŃSKI Kierzkowo 22A, 84-210 Choczewo

Poland [email protected] [email protected]

On A Generalization Of The Boussinesq Approach For Modeling The Flow In Unconfined Aquifers

Figure 1: Schematic illustration of the model geometry

We present an extension to the Boussinesq approach for modeling the flow in unconfined aquifers. The proposed technique is suitable for predicting the geometrical location of the free-surface elevation and both velocity components at the free-surface. Our formulation allows a closed-form solution of the nonlinear problem and an exhaustive discussion of its existence. In terms of the model proposed the standard Boussinesq approach based on the Dupuit assumption is treated as an approximation for the regional gradient approaching zero as the limit. Comparison with the closed-form solution developed by Schmitz and Edenhofer is also presented. TECHNICAL REPORT 2008/2009

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

2

1. Introduction We consider the steady state flow in a homogeneous, isotropic, unconfined aquifer represented by a 2-D vertical slab of length L, bounded by two-constant-head boundary conditions HS and HE (see Figure 1). Because the geometrical location of the free-surface is not known a priori, the boundary value problem describing the system above-mentioned is nonlinear (cf., Kufner and Fučik 1980). To the best of our knowledge, this problem has not been solved in a closed-form as yet. We note in passing that Schmitz and Edenhofer 2000 developed the closedform solution for the flow in a horizontally infinite, and vertically semi-infinite unconfined aquifer with non-symmetric recharge/drainage conditions. They used conformal mapping, transformation procedures, and an inversion of the resulting integral equation for calculating the location of the free-surface elevation. In this note we follow the line presented by Baiocchi and Capelo 1984, showing that the free-surface elevation of the problem in the Figure 1 is described by a strictly decreasing function. Some formal elements of the proof presented are used to construct the model that generalized the well-known Boussinesq approach in such a way that the vertical flow conditions are incorporated. Furthermore, we also extend the technique proposed by Schmitz and Edenhofer 2000 for calculating the velocity field at the free-surface. Thus the present note is organized as follows: in section 2 we rediscussed some basic issues related to the Boussinesq approach, pointing out its formal inconsistence, section 3 shows relations describing the behavior of the velocity field at the free-surface for the nonlinear problem shown in Figure 1, in section 4 we demonstrate the closed-form solution of the model proposed and discuss the issues related to its existence, section 5 presents the comparison of the model proposed with the approach suggested by Schmitz and Edenhofer 2000. Finally, in section 6 conclusions are drawn. 2. Standard Boussinesq approach The common way to overcome the formal complications related to the nonlinearity of the

above problem is the use of the Boussinesq approximation which is based on the so-called Dupuit assumption (see Bear 1979). Unfortunately, as was already pointed out by Baiocchi and Capelo 1984 the Dupuit assumption is formally inconsistent. Thus, all applications based on it should not necessarily be considered as guidelines. The Boussinesq flow equation and appropriate boundary conditions can be written in the following form, L hB(x) = 0,

0≤x ≤L

(1)

hB(0) = HS, hB(L) = HE,

HE ≠ 0

(2)

where L ≡ d2/dx2 + ε (1 – ε fS(x))-1[d/dx]2, ε = HS-1 and z = hB(x) = HS – fS(x) denotes the approximate geometrical location of the free-surface elevation obtained in terms of Boussinesq approach. Here and throughout the rest of this work the subscript B is used to denote an approximation based on the Dupuit assumption. Note, that in equation (1) the linear operator is perturbed by the nonlinear one, where ε can be considered as a perturbation parameter (see Kato 1966 ). For simplicity the vertical component of the infiltration/exfiltration velocity at the freesurface elevation has been assumed to be equal to zero. Equation (1) immediately follows from the dynamic and kinematic boundary conditions at the free-surface (cf., Szymanski 1999), and is often used for simulating the interface evolution between unsaturated and saturated zone (see Calugaru et al. 2003 a, b). The solution of (1-2) takes the form;

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

hB(x) = (A + B x)1/2 ,

0≤x ≤L

3

(3)

where A = HS2, and B = (HE2 – HS2)/L. Note, that d2hB(x)/dx2 < 0 for x ∈ [0, L] and the graph of (3) is convex upwards. Additionally, the function hB(x) has the compact support (cf., BolzanoWeierstrass theorem). Furthermore, the horizontal velocity component obtained in terms of Boussinesq approach is given by the relation vxB(x) = - k dhB(x)/dx, where k denotes the hydraulic conductivity. Hence from (3) follows that vxB(x) = - 0.5 k B (A + B x) -1/2 ,

0 ≤ x ≤ L.

(4)

We refer to the problem (1-2) as the Standard Boussinesq Approach, henceforth abbreviated to (SBA). Its solution has been presented by Marino 1973, and in an approximate form using Adomian's decomposition method by Serano 1995. Moreover, the Boussinesq approach is often used for verifying numerical models (e.g., Onyejekwe et al. 1999, and Taigbenu 2003). The simplified form of (1) is applied by Haneberg 1995 for simulating the groundwater flow across vertical faults, and Schorghofer et al. 2004 consider the above equation for approximating the length scale for the width of the subsurface drainage area. According to Hogarth et al. 1999 the solutions of the Boussinesq equation are applicable in catchment hydrology and base studies, as well as in agricultural drainage problems. As seen the SBA and its different extensions are still in use. From the heuristic point of view, the SBA is assumed to be valid in mildly sloping aquifers characterized by nearly horizontal flow conditions. This assertion is based on the flow being slowly-varying in the x direction, and the pressure being nearly hydrostatic. However, it is clear that the regional gradient m = (HE – HS)/L = tg(Δ) < 0 should be treated as an arbitrarily small quantity approaching zero as a limit to make the Boussinesq approach consistent. Further, vzB(x) is formally undefined for x ∈ [0, L], where vzB(x) denotes the vertical velocity component defined in terms of SBA. Usually, one heuristically assumes that vzB(x) ≈ 0 for x ∈ [0, L], where the sign ≈ denotes that the quantity is approximately equal to zero. To show the inconsistence of the Boussinesq approach, let us assume for a moment that the regional gradient is described by the small, finite quantity, i.e., │m │<< 1. From the examination of SBA one immediately obtains αB = vxB(0 )/vxB(L ) = HE/HS = 1 + γ,

0 < αB ≤ 1

(5)

where γ = m β. According to the Dupuit assumption the relation (5) is simply a volume-balance equation. Dimensionless coefficients αB and β = L/HS describe the kinematic and geometric deepness of an aquifer, respectively. We have introduced the term kinematic deepness for combining together the flow and geometrical characteristics. In fact, taking │m │<< 1 and β = constant we obtain that αB → 1. Further, using the Maclaurin expansion relations (3) and (4) reduce to the following representations, z = hB(x) = HS[1 + γL-1(2 + γ)x]1/2 = mx + HS – 0.5m2βL-1(x2 – Lx) + ... αB→1 (6) vxB(x) = - 0.5 km (2 + γ) [1 + γ L-1 (2 + γ)x]-1/2 = - km [1 - m β L-1 (x – 0.5L) + ... ] αB →1

(7)

Therefore from the SBA follows that for a kinematically deep aquifer the free-surface elevation is very well approximated by linear function. Incidentally, one obtains the same result from (1-2) putting simply ε = 0. Thus from the relation hB(x) ~ m x + HS follows that dhB(x)/dx ~ m. However, at the free-surface should be fulfilled the following; dhB(x)/dx = vzB(x)/vxB(x) (see

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

4

section 3 for details). Because vzB(x) is formally undefined and vxB(x) ≠ 0 for x ∈ [0, L] one obtains that dhB(x)/dx is undefined or heuristically dhB(x)/dx ≈ 0. Hence, the assumption that the regional gradient is treated as a small, finite quantity makes the Boussinesq approximation formally inconsistent. It is not always clear when the SBA will faithfully reproduce the physical phenomena that one attempts to model. Thus our motivation is to develop such an extension of SBA that allows the regional gradient to be considered as a small, finite quantity. We call it the Extended Boussinesq Approach, henceforth abbreviated to (EBA).

3. Velocity components at the free-surface of the nonlinear problem The most ingenious method for reaching the goal mentioned above has been provided by

Baiocchi and Capelo 1984 for proving that the free-surface elevation, considered in section 1, is formally described by a strictly decreasing function. It is motivated by the fact that the velocity components vx(x,h(x)) and vz(x,h(x)) can be defined as nonlinear mappings (cf., Appendix A for details), vx(x,h(x)) = - k dh(x)/dx [1 + (dh(x)/dx)2] -1

0≤x ≤L

(8)

vz(x,h(x)) = - k [dh(x)/dx]2 [1 + (dh(x)/dx)2] -1

0 ≤ x ≤ L,

(9)

where the function z = h(x) describes the free-surface elevation and the functions vx(x,h(x)), vz(x,h(x)) define the horizontal and vertical velocity components for the nonlinear problem, respectively. From (8) and (9) one obtains vz(x,h(x)) = [dh(x)/dx] vx(x,h(x)),

0≤x ≤L

(10)

Expression (10) shows an elementary geometric relation that the velocity vector is tangent to the free-surface elevation. Assuming further, that │dh(x)/dx│<< 1 for x ∈ [0, L] (e.g., in mildly sloping aquifers) we obtain the relations characterized by SBA, namely; vx(x,h(x)) ≈ - k dh(x)/dx, and vz(x,h(x)) ≈ 0. Moreover, from the formal viewpoint, relations (8) and (9) provide a complete description of velocity components at the free-surface elevation.

4. Extended Boussinesq Approach 4.1 Introducing remarks

The proposed extension of SBA is suitable for predicting the geometrical location of the

free-surface elevation, and both velocity components at z = h(x). The EBA is a nonlinear model that is based on the flow characteristics at the free-surface elevation. Although we do not make any simplifying assumptions, our formulation allows the closed-form solution of EBA and an exhaustive discussion of existence problems. Note, that the problems related to the extension of the Boussinesq approach are not new. For instance, Dagan 1967 has applied the perturbation technique to obtain the asymptotic expansion for improving the approximation based on the Dupuit assumption. His method is similar to Friedrichs' technique used in the theory of shallow water-waves. Some remarks on Dagan's approach are presented by Keady 1988. Parlange et al. 1984 used second-order theory to describe the propagation of steady periodic motion in a porous

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

5

medium, and Jeng et al. 2005 obtained the higher-order perturbation solution for calculating the groundwater head fluctuations in coastal aquifers. Recently, Knight 2005 for improving the Dupuit–Forchheimer groundwater free-surface approximation has assumed that the vertical velocity component is zero at the bottom of an aquifer, and is proportional to height above its bottom. He obtained the solution that is more accurate than the Dupuit–Forchheimer expression for the free surface. Let us now describe the basic difference between the SBA and the EBA. It is easy to observe that the solution of SBA formally exists for an arbitrary value of the regional gradient. In the case of EBA there exists an admissible set of values for the regional gradient that depends on the geometrical characteristics of an aquifer. It means that if a regional gradient considered for predicting the flow in an aquifer does not belong to the above-mentioned set, the solution of EBA does not exist. Analogous features arise when studying different branches of nonlinear hydrodynamics. For example, using a perturbation technique Cokelet 1977 has shown that the height of surface water wave cannot be arbitrary and is strictly related to its stability. Schmitz and Edenhofer 2000 found the restriction to the solution between parameters of their model (see section 5 for details). Thus it is not surprising, that for the EBA one also observes similar phenomenon. 4.2 Proposed model First of all let us summarize the assumption required: 1. the regional gradient m is considered to be finite and sufficiently small i.e, │m │<< 1, 2. the value z = 0, which is treated as a horizontal impervious bottom of an aquifer, is only used as a reference datum for the definition of potential energy, and 3. the components of the infiltration/exfiltration velocity at the free-surface elevation have been assumed to be equal to zero. The physical meaning of the second assumption can be explained in the following way. We are interested to develop the model using parameters related to flow characteristics at the freesurface elevation. In this way we only need an information that the bottom of an aquifer does exist, but we are not obliged to consider the exact value for the parameter HS (a characteristic depth of an aquifer). For the description of the bottom of an aquifer the idea of a kinematic deepness is adopted. Let us now introduced the EBA. The exact nonlinear differential equation governing the horizontal velocity component of the groundwater flow at the free-surface elevation is given by (8), vx(x,z) [dh(x)/dx]2 + k dh(x)/dx + vx(x,z) = 0

at z = h(x), 0 ≤ x ≤ L

(11)

where the functions vx(x,h(x)) and h(x) are unknown. Incidentally, similar idea was used for describing the interface between fresh and salt groundwater in heterogeneous aquifers (see, Scheid and Schotting 1997). Equation (11) is subject to the formal (cf., the second assumption) boundary conditions, h(0) = HS, h(L) = HE,

HE ≤ z ≤ HS

(12)

According to the first assumption the form of an unknown function vx(x,h(x)) is defined by vx(x,h(x)) = - 0.5 k b (a + b x) -1/2 , 0 ≤ x ≤ L

(13)

where a > 0 and b < 0 are unknown, real constants. The existence of the bottom of an aquifer

6

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

can be described by the coefficient α; α = Vc(0,HS)/Vc(L,HE),

0≤ α ≤1

(14)

where Vc(x,z) = [vx(x,z)2 + vz(x,z)2]1/2 at z = h(x). According to (5) we can call α the kinematic deepness of an aquifer. However, from the formal point of view α can also be considered as a parameter describing the kinematic characteristics at the free-surface (see section 5). For calculating the vertical velocity component one should use relations (9) or (10). Furthermore, from the heuristic assumption that the word “mildly” is equivalent to the term “finite and small” we suppose that (13) is a reasonable approximation for the horizontal velocity component at the free-surface. In conclusion, we note that for the use of EBA two parameters should be known; m and α. Next, putting dh(x)/dx = c the equation (11) reduces to the quadratic equation with real coefficients. Using the fundamental theorem of calculus its solution takes the form x

h(x) = 0.5∫ [k/vx(ξ,h(ξ))] {-1 + μ [1 – 4 (vx(ξ,h(ξ))/k)2]1/2} dξ + HS,

0≤x ≤L

(15)

0

where an integration constant has been determined by the boundary condition (12), 0 ≤ vx(x,h(x)) ≤ 0.5k, and μ ∈ {−1,1}. Furthermore, the sign of the square root (μ = 1) was chosen in such a way that for vx(x,h(x)) << 1 for x ∈ [0, L] from (15) follows that; vx(x,h(x)) ≈ - k dh(x)/dx. The formal problems related to the parameter μ are described in Appendix B. Note in passing that the inequality presented above does not depend on the relation (13). 4.3 Linearized solution of EBA If vx(x,h(x)) << 1 for x ∈ [0, L] one obtains the first-order approximation of the problem (11-14) by expressing the square root (see (15)) in the form of the power series. In fact, using the ordinary Maclaurin's series (1- c)1/2 = 1 – (1/2)c + O(c2) for │c │≤ 1 and assuming that c = 4[vx(x,h(x))/k]2 one can rewrite (15) in the following way x

h(x) = - k

-1

∫ vx(ξ,h(ξ))dξ

+ H S,

0≤x ≤L

(16)

0

Next, comparing (16) to (11) we conclude that the linearization procedure proposed implies; (dh(x)/dx)2 ≈ 0. Thus from (9) follows that vz(x,h(x)) ≈ 0 for x ∈ [0, L]. This results in simpler expression (see (14)) α = vx(0,HS))/vx(L,HE),

0< α ≤1

(17)

Note that from (16) one obtains the linearized version of (11), namely; k dh(x)/dx + vx(x,z) = 0 at z = h(x). Combining the latter differential equation with relations (12), (13) and (17) one obtains the linearized version of EBA. It is easy to find its closed-form solution. In fact, substituting (13) into (16), carrying out an integration and using (12) we obtain HE = (a + b L)1/2 – a1/2 + HS

(18a)

Further, from (17) follows that α = (a + b L)1/2/a1/2

(18b)

7

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

The solution of the nonlinear system (18a,b) is given by a = (HE – HS)2/(α -1)2 > 0,

b = [(α +1)/(α -1)][(HE – HS)2/L] < 0,

α≠1

(19)

Next, assuming in (19) that α = αB (see (5)) one obtains that a = A and b = B. In conclusion, the linearization procedure proposed implies that the EBA reduces to the SBA. 4.4 Solution of EBA (the iterative Newton-Raphson method) The procedure for solving the EBA is relatively simple. One can proceed in much the same fashion as in the preceding section, except that the calculation is more complicated. Hence, inserting (13) into (15) we carry out an integration. Next, using the boundary condition (12) we obtain the first equation needed for determining the unknown coefficients a and b. The second one is obtained by means of (14). Finally, we have F1(c1,f1) = 0

(20a)

F2(c1,f1) = 0

(20b)

where

F1(c1,f1) = (c1 + f1)3/2 – c13/2 – (c1+ f1 – f12)3/2 + (c1 -f12)3/2 - 1.5m f12, F2(c1,f1) = α2 {1 – [1 – f12/(c1 + f1)]1/2} -1 + [1 – f12/c1]1/2, c1 = a/L2 and f1 = b/L.

(20c)

The nonlinear system (20a,b) is more complicated as the one considered above. Using the power series procedure applied above we establish its approximate solution as follows (cf., Olver 1974); c1 ≈ m2/(1 – α)2 > 0,

f1 ≈ - m2 (1 + α)/(1 – α) < 0,

α≠1

(21)

Note that (21) is equivalent to (19), because we only used the first-order approximation. The procedure for solving (20a,b) described in this section is based on the iterative Newton-Raphson method which is locally quadratic converged, assuming that the initial guess is sufficiently close to the solution (e.g., Dahlquist and Björck 1974). We find after some straightforward algebraic manipulations that the roots of the system (20a,b) are given by f1 (i+1) = [PF1c1 (i) Le2 (i) - PF2c1 (i)]/[PF1c1 (i) PF2f1 (i) – PF2c1 (i) PF1f1 (i)]

(22a)

c1 (i+1) = [Le1 (i) – PF1f1 (i) f1 (i)]/PF1c1 (i)

(22b)

i = 1, 2, 3, ..., N

where PF1c1(i) = α2{-0.5[1 – f1(i)2/(c1(i) + f1(i))]-1/2 f1(i)2 (c1(i) + f1(i))2} + 0.5[1 - f1(i)2/c1(i)2]-1/2 (f1(i)2/c1(i)2), PF1f1(i) = α2 {0.5[1 – f1(i)2/(c1(i) + f1(i))]-1/2 [2f1(i) (c1(i) + f1(i)) – f1(i)2]/(c1(i) + f1(i))2} – (f1(i)/c1(i)) [1- f1(i)2/ c1(i)2]-1/2, PF2c1(i) = 1.5 [(c1(i) + f1(i))1/2 – c1(i)1/2 – (c1(i) + f1(i) – f1(i)2)1/2 + (c1(i) – f1(i)2)1/2], PF2f1(i) = 1.5 [(c1(i) + f1(i))1/2 – (c1(i) + f1(i) – f1(i)2)1/2(1- 2f1(i)) + (c1(i) – f1(i)2)1/2(-2f1(i))] – 3m f1(i) , Le1(i)

= - F1(c1(i),f1(i)) + c1(i) PF1c1(i) + f1(i) F1(c1(i),f1(i),

Le2(i) = - F2(c1(i),f1(i)) + c1(i) PF2c1(i) + f1(i) F2(c1(i),f1(i)), and

8

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

c1(1), f1(1) are described by (21). The (22a,b) define a family of approximations of increasing accuracy that are indexed with the integer i. The algorithm terminates when the residual measure is smaller than some predetermined criterion. Additionally, in (22a,b) the notation, for example, PF1c1(i) ≡ ∂F1(c1(i),f1(i))/∂c1(i) convenient for partial differential equations has been adopted. From extensive calculations the following observation may by noted; for │m │< 0.1 and 0.5 ≤ α < 1 the solution of (20a,b) is well approximated by (21). Thus it is reasonable to use (21) as an initial guess for starting the iteration procedure proposed. We postpone the description of this problem until the next section. However, in practice one can also use the computer algebra package Mathematica (see Wolfram 1988, the command FindRoot). Furthermore, having defined the roots of (20a,b), by means of (13) and (15) we can calculate the free-surface elevation for the EBA as follows; h(x) = (2/3)b-2 [(a + b x)3/2 – a3/2 – (a – b2 + b x)3/2 + (a – b2)3/2] + HS,

(23)

and the vertical velocity component is given by (10), where dh(x)/dx can be obtained from (15) using the fundamental theorem of calculus. 4.5 Definition of an admissible set To find the admissible set D for the regional gradient it is necessary to consider the domains of the real-valued functions F1(c1,f1) and F2(c1,f1), that is, the sets on which the functions are defined. In this sense the domain of the real-valued square root, as a component of the above-mentioned functions, cannot exceed the non-negative reals. Thus taking into account all functional components of F1(c1,f1) and F2(c1,f1) one obtains the following inequality c1 + f1 – f12 ≥ 0

(24)

Considering now the limiting case c1 + f1 – f12 = 0 the equation (20b) can be rewritten as follows; f1g = - [c1g(1 – H)]1/2

(25)

where H = (1 – α2)2. Substituting (25) into the limiting case of (24) one obtains c1g = (1 – H)/H2.

(26)

Relations (25) and (26) are the closed-form solutions of (20a,b) for the limiting case. Here and throughout the rest of this work the subscript g is used to distinguish the solution of (20a,b) from the limiting one. Next, from (26) and (25) follows that f1g = - (1 – H)/H

(27)

As (27) and (26) should fulfill (20a) we can rewrite the latter in the form 1.5 mg(α) f1g2 = (c1g + f1g)3/2 – c1g3/2 + (c1g - f1g2)3/2

(28)

The mg(α) denotes the limiting regional gradient. Finally, substituting (26) and (27) into (28) we obtain the relation for the limiting regional gradient in a simple form mg(α) = (2/3) [(1 – H)-1/2H-1][(1 – H )3/2 – 1 + H3/2]

(29)

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

9

Further, from α = 1 follows that H = 0. Hence using l' Hospital's rule one obtains from (29) the limits; lim mg(α) = -1, (i.e., Δ = 135o) α→1

(30a)

and from α = 0 follows that H = 1, yielding lim mg(α) = 0, (i.e., Δ = 0o) α→0

(30b)

Finally, the admissible set D for the regional gradient is given by D: = {(m,α): 0 ≤ α ≤ 1 │m │ ≤ │mg(α)│}

(31)

In conclusion, we now require condition (31) to be satisfied since otherwise the solution of EBA does not exist. It means that in terms of EBA the regional gradient cannot be “arbitrary” as in the case of SBA. Graphic representations of the admissible set D and the relation (29) are shown in Figure 2.

Figure 2: Admissible set for the regional gradient

The gray color area over the curve presented is called the admissible set for the regional gradient. Although great attention has been paid to D, it must be emphasized that D obtained should not be seen as intrinsic to the problem of the flow in unconfined aquifers; the form of D is purely the result of the assumption made in section 4.2 (see (13)). Further, to find the limiting values for the horizontal velocity component we should combine (20c) with (25) and (26). This implies that bg = - (1 – H)H-1L, and

ag = (1 – H)H-2L2

(32)

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

10

Next, substituting (32) into (13) we obtain for x ∈ [0, L] that vxg(x,hg(x))/k = 0.5 (1 – H)H-1L [(1 – H)H-2 L2 – (1 - H)H-1 L x]-1/2

(33)

Evaluating (33) we find vxg(L,HE) = 0.5 k, and

(34a)

for α → 1

vxg(0,HS) → 0.5 k, and

(34b)

for α → 0

vxg(0,HS) → 0

(34c)

Note that for the case x = L the result obtained does not depend on the coefficient α. In addition, for α → 1, the horizontal velocity component is equal at the free-surface elevation at both locations x = 0 and x = L. Consider now the behavior of the vertical velocity component. The relation (10) for the limiting case can be rewritten as follows; vzg(x,hg(x)) = [dhg(x)/dx] vxg(x,hg(x))

(35)

where dhg(x)/dx is obtained from (15) giving dhg(x)/dx = 0.5 [k/vxg(x,hg(x))] {-1 + [1- 4(vxg(x,hg(x))/k)2]1/2}.

(36)

Further, from (33) and (35) one obtains that at x = L vzg(L,HE) = - 0.5k

(37)

Again, the result obtained does not depend on the coefficient α and Vcg(L,HE) = (0.5)1/2k ≈ 0.71k. At x = 0 we find that vzg(0,HS) = 0.5 k (-1 + H1/2)

(38)

Finally, from (38) we have vzg(0,HS) = - 0.5 α2k, and

(39a)

for α → 1

vzg(0,HS) → - 0.5 k, and

(39b)

for α → 0

vzg(0,HS) → 0

(39c)

In addition, from (34b) and (39b) follows that for α = 1; Vcg(L,HE) = Vcg(0,HS). Let us now check the behavior of the function dhg(x)/dx in limiting cases at x = L and x = 0, respectively. From (36) we obtain that dhg(L)/dx = -1

(40a)

dhg(0)/dx = - α2/(1 - H)1/2

(40b)

Considering further (40b) and applying the l' Hospital's rule one obtains the following limits; if α → 0 then dhg(0)/dx → 0, and if α → 1 then dhg(0)/dx → -1. Hence from our considerations follows that, for α → 1, the free-surface elevation of EBA is very well approximated by the

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

11

function h(x) = m x + HS. The numerical examples are presented in Appendix C. Next, using (8) and (9) one calculates vx(x,h(x) → m x + HS) → - k m (1 + m2) -1, α→1

0≤x ≤L

(50a)

vz(x,h(x) → m x + HS) → - k m2 (1 + m2) -1, α→1

0≤x ≤L

(50b)

Therefore, for α = 1, the EBA reduces to the so-called Local Subsurface Solution (LSS) introduced by Szymanski 1998, 1999. Having defined the admissible set D and the closed-form solution (26-27) of the nonlinear system (20a,b), we can now examine the issues related to the convergence of the iteration procedure proposed in section 4.4. It can be proved that the Newton-Raphson method is usually convergent if the initial estimates to the roots are properly chosen. Unfortunately, if │m│→ │mg(α)│ initial estimates (21) cannot guarantee the convergence of the method used. In fact, calculating the distances │c1(1) – c1g│and │f1(1) – f1g│ for m = mg(α) we find that; │c1(1) – c1g│=│(1 - α2) -1 (1 - α ) -1 (- mg(α)2 (1 + α) - (1 – H) H-1(1 + α)-1)│

(51a)

│f1(1) – f1g│ = │(1 + α) (1 - α ) -1 mg(α)2 - (1 – H) H -1)│

(51b)

From (51a,b) one obtains: if α → 1 then │c1(1) – c1g│→ ∞, and │f1(1) – f1g│→ ∞. Thus to overcome this technical difficulty we suggest the following strategy for finding the solution of (20a,b) when │m │→ │mg(α)│. We should start the iteration procedure, for a given α, using such a small value of │m│that the iteration procedure suggested is convergent. Next, slightly increasing │m│, one should use the solution previously obtained as an initial estimate for c1(1) and f1(1). Note that the technique proposed above is based on homotopy methods. The idea is that as the parameter changes, the solutions trace out a path from a solution of the starting parameter to a solution of the original one. At each stage, the current system is normally solved by the Newton-Raphson method suggested above. 4.6 Solution of EBA (2-D bracketing search method) The method we introduce in this section can be called the 2-D bracketing search technique because one starts with the knowledge that the roots of (20a,b) lie in some region of the admissible set D for the regional gradient, so that it is only necessary to refine this knowledge. Thus it is clear that we do not look for the solution of the nonlinear system (20a,b) on the (f1,c1)-plane but on the (m,α)-plane. The search technique mentioned above is based on minimizing the Cartesian distance defined as; d = {(m - mp)2 + (α - αp)2}1/2

(52)

where mp and αp are the parameters of EBA (see section 4.2), and m, α are given by m = [(c1 + f1)3/2 – c13/2 – (c1+ f1 – f12)3/2 + (c1 -f12)3/2 ]/1.5f12,

(53a)

α = {{1 - [1 – f12/c1]1/2}/{1 – [1 – f12/(c1 + f1)]1/2}}1/2,

(53b)

where f1 ∈ [f1g, 0) (see (27)) and c1 = c1g (see (26)). For the given values of α ∈ (0,1), f1 ∈ [f1g, 0) and c1 = c1g we can calculate the distance (52) using relations (53a,b). Hence, the procedure combining (53a,b) with the partitions of intervals [f1g, 0) and (0,1) generates the system of curves

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

12

on the (m,α)-plane which forms a covering of the admissible set D. An example of such curves is presented in the Figure 3. Three curves of the mentioned system are shown for α = 0.5, 0.9 and 0.1, respectively. The set called mg.dat is defined by the relation (29).

Figure 3: An example of elements of the covering set The above-described geometrical approach was implemented as a numerical procedure for solving the nonlinear system (20a,b). Unfortunately, it is not the most efficient method. Thus we use the bracketing search technique for localizing the initial estimate for c1(1) and f1(1) with appropriate accuracy, so that the Newton-Raphson method can directly be applied.

5. Comparison of EBA with Schmitz and Edenhofer 2000 5.1 Example of Weber 2003

Following

Schmitz and Edenhofer 2000 and Weber 2003 we consider the flow in a horizontally infinite, and vertically semi-infinite 2-D slab representing unconfined aquifer. The flow is induced by the infiltration/exfiltration velocity at the free-surface elevation defined as; N (x) := {[p for x ∈ (- ∞, -b)], [r for x ∈ [-b, -a]], [0 for x ∈ (-a, a)], [-r for x ∈ [a, b]], [p for x ∈ (b, ∞)]} where a, b, p, r are real numbers. The free-surface elevation is representing by the following integral (see Schmitz and Edenhofer 2000), ∞ -1 hSE(x) = (- k π) PV ∫ N (ξ) ln │u (ξ) – u(x) │dξ (54) -∞ where PV denotes the principal value integral, and x u (x) = x – k

-1

∫ N (ξ) dξ

-∞

∞ and

∫ N (ξ) dξ -∞

= 0.

(54a)

13

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

Here and throughout the rest of this work the subscript SE is used to distinguish the solution obtained by Schmitz and Edenhofer 2000 from the EBA. Next, following Weber 2003, we assume that a = 2 [length], b = 3 [length], k = 2 [length/time], p = 0 [length/time], and r = 1 [length/time]. Note, that the inequality N(x) < k for x ∈[-b, -a] and x ∈ [a, b] should be fulfilled to guarantee the existence of the solution. For the comparison mentioned the interval [3, ∞) is used because in terms of EBA the flow is only induced by the regional gradient and the graph of the free-surface elevation is convex upwards. Using relations (54) and (54a), we obtained, hSE(x) = 0.5 π -1 {3-1 (6 – 2x) ln│6 – 2x │- (6 + 2x) ln│6 + 2x │- 3-1 (3 – 2x) ln│3 – 2x│ (5 + 2x) ln│5 + 2x │} for x ∈ [3, ∞)

(55)

where hSE [length] and for x ∈ [3, ∞) hSE(x) < 0. From (55) follows that dhSE(x)/dx = 0.5 π -1{(-2/3) ln│6 – 2x │- 2 ln│6 +2x │+ (2/3) ln│3–2x │+2 ln│5 + 2x │} for x ∈ [3, ∞)

(56)

Having defined the function dhSE(x)/dx in a closed-form, we can calculate both velocity components at the free-surface using relations (8) and (9). Further, it is clear that the regional gradient is generated by the accretion function N(x) defined outside of the interval [3, ∞). Thus we choose four intervals for the direct comparison, namely; x ∈ [200, 206], x ∈ [20, 26], x ∈ [100, 500] and x ∈ [600, 1000] to test the influence of the parameters m and α. on the results obtained. Additionally, for simplifying the comparison the following transformations have been introduced, HSE(x) = hS + hSE(x), and x1 = (lE - lS) - (x - lS)

(57)

where hS, lS, lE describe the so-called transformation parameters. 5.2 Comparison for the interval [200, 206] The parameters used in this case are as follows; hS = 10, lS = 200, lE = 206 and L = 6. Further, m = ( hSE(200) - hSE(206))/L = -1.92736E-05. According to (14) we calculated α = [vx(206, hSE(206))2 + vz(206, hSE(206))2]1/2/[vx(200, hSE(200))2 + vz(200, hSE(200))2]1/2 = 0.94265. Table 1. Comparison of the free-surface elevations obtained from (57) for the Schmitz and Edenhofer model, and from (23) for the EBA.

x1

0

1

2

3

4

5

6

HSE(x1)

9.99614

9.99612

9.99610

9.99608

9.99606

9.99604

9.99602

h(x1)

9.99614

9.99612

9.99610

9.99608

9.99606

9.99604

9.99602

0

0

0

0

0

0

0

0

0

0

0

0

0

0

9.99614

9.99612

9.99610

9.99608

9.99606

9.99604

9.99602

Eh [%] EhB[%] hB(x1)

14

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

where the relative errors are defined as follows; Eh = │(HSE(x1) – h(x1))/hSE(x1)│100%, EhB = │(HSE(x1) – hB(x1))/hSE(x1)│100%. Results obtained from the SBA are presented in italics on last line. We used (6) where β = 0, m = -1.92736E-05 and HS = 9.99614.

Table 2. Comparison of the normalized horizontal velocity component at the free-surface obtained from (56) and (8) for the Schmitz and Edenhofer model, and from (13) for the EBA.

x1

0

1

2

3

4

5

6

S-1 [vxSE(x1,.)/ 1.87128 k]

1.88957

1.90812

1.92694

1.94605

1.96544

1.98513

S-1[ vx(x1,.)/k] 1.87210

1.88973

1.90786

1.92653

1.94575

1.96556

1.98599

Evx [%]

0.05

0.01

0.01

0.02

0.02

0.01

0.04

EvxB [%]

3.00

2.16

1.00

0.02

1

1.94

2.91

S-1[ vxB(x1)/k]

1.92736

1.92736

1.92736

1.92736

1.92736

1.92736

1.92736

where S denotes the scaling factor, and S = 1.0E-05. Relative errors are defined as follows; Evx = │(vxSE(x1,.) – vx(x1,.))/vxSE(x1,.) │100%, EvxB = │(vxSE(x1,.) – vxB(x1,.))/vxSE(x1,.) │100%. Results obtained from the SBA are presented in italics on last line. We used (7) where β = 0, m = -1.92736E-05.

15

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

Table 3. Comparison of the normalized vertical velocity component at the free-surface obtained from (56) and (9) for the Schmitz and Edenhofer model, and from (9) for the EBA. x1

0

1

2

3

4

5

6

S-1 [vzSE(x1,.)/k]

-3.50170

-3.57046

-3.64091

-3.71312

-3.78712

-3.86297

-3.94073

Evz [%]

0.09

0.02

0.03

0.04

0.03

0.01

0.09

S-1 [vz(x1,.)/k]

-3.50476

-3.57107

-3.63993

-3.71151

-3.78595

-3.86344

-3.94417

where S = 1.0E-10. Relative error is given by Evz = │(vzSE(x1,.) - vz(x1,.))/vzSE(x1,.) │100%.

16

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

5.3 Comparison for the interval [20, 26] The parameters used in this case are as follows; hS = 10, lS = 20, lE = 26 and L = 6. Further, m = ( hSE(20) - hSE(26))/L = -1.51856E-03 . According to (14) we calculated α = [vx(26, hSE(26))2 + vz(26, hSE(26))2]1/2/[vx(20, hSE(20))2 + vz(20, hSE(20))2]1/2 = 0.59094. Table 4. Comparison of the free-surface elevations obtained from (57) for the Schmitz and Edenhofer model, and from (23) for the EBA.

x1

0

1

2

3

4

5

6

HSE(x1)

9.96957

9.96835

9.96704

9.96561

9.96405

9.96234

9.96046

h(x1)

9.96957

9.96833

9.96700

9.96559

9.96405

9.96236

9.96046

Eh [%]

0.0

0.09

0.12

0.06

0.0

0.05

0.0

EhB[%]

0.0

0.95

1.55

1.74

1.53

0.96

0.0

hB(x1)

9.96957

9.96805

9.96653

9.96501

9.96350

9.96198

9.96046

Results obtained from the SBA are presented in italics on last line. We used (6) where β = 0, m = -1.51856E-03 and HS = 9.96957. Relative errors Eh and EhB are defined in section 5.2 (cf., Table

17

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

1).

Table 5. Comparison of the normalized horizontal velocity component at the free-surface obtained from (56) and (8) for the Schmitz and Edenhofer model, and from (13) for the EBA. x1

0

1

2

3

4

5

6

S-1 [vxSE(x1,.)/k]

1.1675 2

1.26289

1.37048

1.49249

1.63163

1.79128

1.97571

S-1 [vx(x1,.)/k]

1.2079 7

1.27934

1.36507

1.47072

1.60543

1.78556

2.04414

Evx [%]

3.46

1.30

0.39

1.46

1.61

0.32

3.46

EvxB [%]

30.07

20.24

10.80

1.75

6.93

15.22

23.14

S-1 [vxB(x1)/k]

1.51856

1.51856

1.51856

1.51856

1.51856

1.51856

1.51856

where S = 1.0E-03. Results obtained from the SBA are presented in italics on last line. We used (7) where β = 0, m = -1.51856E-03. Relative errors Evx and EvxB are defined in section 5.2 (cf., Table 2).

18

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

Table 6. Comparison of the normalized vertical velocity component at the free-surface obtained from (56) and (9) for the Schmitz and Edenhofer model, and from (9) for the EBA.

x1

0

S-1 [vzSE(x1,.)/k] -1.36312

1

2

3

4

5

6

-1.59487

-1.87821

-2.22753

-2.66222

-3.20870

-3.90343

Evz [%]

7.05

2.62

0.78

2.90

3.19

0.64

7.05

S-1 [vz(x1,.)/k]

-1.45919

-1.63671

-1.86342

-2.16302

-2.57742

-3.18825

-4.17851

where S = 1.0E-06. Relative error Evz is defined in section 5.2 (cf., Table 3).

19

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

5.4 Comparison for the interval [100, 500] The parameters used in this case are as follows; hS = 10, lS = 100, lE = 500 and L = 400. Further, m = ( hSE(100) - hSE(500))/L = -1.58753E-05. According to (14) we calculated α = [vx(500, hSE(500))2 + vz(500, hSE(500))2]1/2/[vx(100, hSE(100))2 + vz(100, hSE(100))2]1/2 = 4.01242E-02. Table 7. Comparison of the free-surface elevations obtained from (57) for the Schmitz and Edenhofer model, and from (23) for the EBA. x1

0

50

100

150

200

250

300

350

400

HSE(x1)

9.9984

9.9982

9.9980

9.9977

9.9973

9.9968

9.9960

9.9947

9.9921

h(x1)

9.9984

9.9980

9.9975

9.9970

9.9965

9.9959

9.9951

9.9941

9.9921

Eh [%]

0.0

11.11

25.00

30.43

25.93

28.13

22.50

11.32

0.0

EhB[%]

0.0

33.33

60.00

73.91

74.07

71.88

57.50

33.96

0.0

hB(x1)

9.9984

9.9976

9.9968

9.9960

9.9953

9,9945

9.9937

9,9929

9.9921

Results obtained from the SBA are presented in italics on last line. We used (6) where β = 0, m = -1.58753E-05 and HS = 9.9984. Relative errors Eh and EhB are defined in section 5.2 (cf., Table 1).

20

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

Table 8. Comparison of the normalized horizontal velocity component at the free-surface obtained from (56) and (8) for the Schmitz and Edenhofer model, and from (13) for the EBA.

x1

0

50

100

150

200

250

300

350

400

0.31802 0.39258 0.49680 0.64878 0.88288 1.27100 1.98513 3.52682 S-1 [vxSE(x1,.)/k ]

7.92592

0.82562 0.88252 0.95308 1.04382 1.16666 1.34642 1.64726 2.32214 S-1 [vx(x1,.)/k]

20.576

Evx [%]

159.51

124.80

91.84

60.89

32.14

5.93

17.02

34.16

159.61

EvxB [%]

399.19

304.38

219.55

144.69

79.81

24.90

20.03

54.99

79.97

S-1 [vxB(x1)/ 1.587 k]

1.587

1.587

1.587

1.587

1.587

1.587

1.587

1.587

where S = 1.0E-05. Results obtained from the SBA are presented in italics on last line. We used (7) where β = 0, m = -1.58753E-05. Relative errors Evx and EvxB are defined in section 5.2 (cf., Table 2).

21

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

Table 9. Comparison of the normalized vertical velocity component at the free-surface obtained from (56) and (9) for the Schmitz and Edenhofer model, and from (9) for the EBA.

x1

0

50

100

150

200

250

300

350

-0.0101 -0.0154 -0.0246 -0.0421 -0.0779 -0.1615 -0.3941 -1.2438 S-1 [vzSE(x1,.)/k ]

Evz [%]

575.25

639.09

269.11

158.91

74.71

12.26

31.16

56.65

-0.0682 -0.0779 -0.0908 -0.1090 -0.1361 -0.1813 -0.2713 -0.5392 S-1 [vz(x1,.)/k]

where S = 1.0E-09. Relative error Evz is defined in section 5.2 (cf., Table 3).

400 -6.2820

573.97 -42.339

22

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

5.5 Comparison for the interval [600, 1000] The parameters used in this case are as follows; hS = 10, lS = 600, lE = 1000 and L = 400. Further, m = ( hSE(600) - hSE(1000))/L = -1.32548E-06. According to (14) we calculated α = [vx(1000, hSE(1000))2 + vz(1000, hSE(1000))2]1/2/[vx(600, hSE(600))2 + vz(600, hSE(600))2]1/2 = 0.36011. Table 10. Comparison of the free-surface elevations obtained from (57) for the Schmitz and Edenhofer model, and from (23) for the EBA. x1

0

50

100

150

200

250

300

350

400

HSE(x1)

9.99920 9.99916 9.99912 9.99906 9.99901 9.99894 9.99886 9.99878

9.99867

h(x1)

9.99920 9.99916 9.99911 9.99906 9.99900 9.99894 9.99886 9.99878

9.99867

Eh [%]

0.0

0.0

1.14

0.0

1.01

0.0

0.0

0.0

0.0

EhB[%]

0.0

3.57

5.68

6.38

8.08

6.60

5.26

3.28

0.0

hB(x1)

9.99920

9.99913

9.99907

9.99900

9.99893

9.99887

9.99880

9.99874

9.99867

Results obtained from the SBA are presented in italics on last line. We used (6) where β = 0, m = -1.32548E-06 and HS = 9.99920. Relative errors Eh and EhB are defined in section 5.2 (cf., Table 1).

23

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

Table 11. Comparison of the normalized horizontal velocity component at the free-surface obtained from (56) and (8) for the Schmitz and Edenhofer model, and from (13) for the EBA.

x1

0

50

100

150

200

250

300

350

0.79541 0.88132 0.98194 1.10082 1.24269 1.41385 1.62297 1.88217 S-1 [vxSE(x1,.)/ k]

400 2.20881

0.90142 0.95485 1.01908 1.09829 1.19940 1.33481 1.52968 1.84591 2.50317 S-1 [vx(x1,.)/k ]

Evx [%]

13.33

8.34

3.78

0.23

3.48

5.59

5.75

1.93

13.33

EvxB [%]

66.64

50.40

35.03

20.41

6.66

6.25

18.33

29.58

39.99

1.325 S-1 [vxB(x1)/k]

1.325

1.325

1.325

1.325

1.325

1.325

1.325

1.325

where S = 1.0E-06. Results obtained from the SBA are presented in italics on last line. We used (7) where β = 0, m = -1.32548E-06. Relative errors Evx and EvxB are defined in section 5.2 (cf., Table 2).

24

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

Table 12. Comparison of the normalized vertical velocity component at the free-surface obtained from (56) and (9) for the Schmitz and Edenhofer model, and from (9) for the EBA.

x1

0

50

100

150

200

250

300

-0.6327 -0.7767 -0.9642 -1.2118 -1.5443 -1.9990 -2.6340 S-1 [vzSE(x1,. )/k]

Evz[%]

28.43

17.38

7.71

0.46

6.84

10.87

11.17

350

400

-3.5426

-4.8788

3.82

28.43

-0.8126 -0.9117 -1.0385 -1.2062 -1.4386 -1.7817 -2.3399 -3.4074 S-1 [vz(x1,.)/ k]

where S = 1.0E-12. Relative error Evz is defined in section 5.2 (cf., Table 3).

-6.2659

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

25

6. Summary and Conclusions A new

method, called the Extended Boussinesq Approach (EBA), was developed for predicting the flow in regional, unconfined aquifers. The version presented here is simply an extension of the well-known Standard Boussinesq Approach (SBA) based on the Dupuit assumption. An important benefit of the proposed approach is that we can calculate not only the geometrical location of the free-surface elevation but also both velocity components at the freesurface and discussed the limiting case, describing the formal characteristics of the solution obtained. In the method proposed, there is no need to use any linearization procedure with respect to the boundary conditions at the free-surface, and the closed-form solution is obtained. The EBA is characterized by two parameters, namely; regional gradient m and the so-called kinematic deepness of an aquifer denoted by α . In Tables 1–12 the examples are presented to examine the accuracy of EBA against the closed-form solution developed by Schmitz and Edenhofer 2000. In this context the example proposed by Weber 2003 is used. Comparison with the SBA is also included. In all cases discussed the EBA supplies better results as the SBA. For small regional gradients, i.e., less than 0.001, and 0.6 < α < 1 the EBA is a reasonable approximation to the model suggested by Schmitz and Edenhofer 2000, and seems to be a simplest tool available to most hydrologists for predicting the flow characteristics in unconfined aquifers. For α → 0 the EBA strongly overestimates both velocity components in comparison to the model proposed by Schmitz and Edenhofer 2000. Note, that the functions h(x), vx(x,h(x)) and vz(x,h(x)) describing the physical characteristics of EBA have the compact supports. Therefore, from the philosophical point of view, the model proposed by Schmitz and Edenhofer 2000 can be considered as the universe with respect to the EBA. The EBA can also be used for analyzing the behavior of different numerical approaches

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

26

often applied in hydrological studies (see Knupp 1996). The formal implementation of EBA is very simple. In general, from a practical point of view, one can utilize any computer algebra package to overcome the problems related to the solution of the set of nonlinear equations. In this context, we suggest the application of the 2-D bracketing search technique described above. Thus, the EBA is only slightly harder to use in comparison to the SBA. However, the main contribution of this work is simply an assumption that the regional gradient should be considered as a finite quantity. It pushes the EBA into the class of nonlinear problems, where the linearization procedure denotes that the regional gradient is taken as an arbitrarily small quantity approaching zero as a limit. In this sense the SBA is considered as the linearized solution of EBA. Thus, formally, the solution of SBA exists for an arbitrary value of the regional gradient. The EBA, however, is characterized by an admissible set for the regional gradient that depends on the kinematic deepness of an aquifer. It means that the solution of EBA does not always exist, implying that the velocity components of the regional flow are bounded; 0 ≤ vx(x,h(x)) ≤ 0.5k and 0 ≤ |vz(x,h(x))| ≤ 0.5k where x ∈ [0, L]. In addition, for the original problem shown in the Figure 1, from (15) follows that the both velocity components are also bounded, and the free-surface elevation is described by a strictly decreasing function.

7. References Baiocchi, C., Capelo, A., 1984. Variational and Quasivariational Inequalities Application to Free Boundary Problems, John Wiley, New York. Bear, J., 1979. Hydraulics of Groundwater, McGraw-Hill, New York. Calugaru, C I., Calugaru, D-G., Crolet, J-M., Panfilov, M., 2003a. Stability of the interface between two immiscible fluids in porous media, Electronic Journal of Differential Equations, 25, pp. 1–10. Calugaru, C I., Calugaru, D-G., Crolet, J-M., Panfilov, M., 2003b. Evolution of fluid-fluid interface in porous media as the model of gas-oil fields, Electronic Journal of Differential Equations, 73, pp. 1–13. Cokelet,E. D., 1977. Steep gravity waves in water of arbitrary uniform depth, Phil. Trans. Roy. Soc. London, A286, 183-230. Dagan, G., 1967. Second-order theory of shallow free-surface flow in porous media, Quart. J. Mech. Applied Math., 20, 617-526. Dahlquist, A., Björck,G., 1974. Numerical Methods, Prentice-Hall. Haneberg, C.W., 1995. Steady state groundwater flow across idealized faults, Water Resour. Res., 31(7), 1815-1820.

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

27

Hogarth, W.L., Parlange, J-Y., Parlange, M-B., Lockington, D.A., 1999. Approximate analytical solution of the Boussinesq equation with numerical validation, Water Resour. Res., 35(10), 3193-3197. Jeng, D.-S., Seymour, B.R., Barry, D.A., Parlange, J.-Y., Lockington, D.A., Li, L., 2005. Steepness expansion for free surface flows in coastal aquifers, J. Hydrol., 309, 85–92. Kato, T., 1966. Perturbation theory of linear operators, Springer-Verlag, Berlin -Heidelberg-New York. Keady, G., 1988. The Dupuit approximation for the rectangular dam problem, Internal report, Mathematics Department, University of Western Australia, Nedlands 6009, Western Australia. Knight, J.H., 2005. Improving the Dupuit-Forcheimer groundwater free surface approximation, Adv. Water Resour., 28, 1048-1056. Knupp, P., 1996. A moving mesh algorithm for 3-D regional groundwater flow with water table and seepage face, Adv. Water Resources,10 (1), 83-95. Kufner, A., Fučik, S., 1980. Nonlinear Differential Equations, Elsvier Publ, AmsterdamOxford-New York. Marino, A.M., 1973. Water-table fluctuation in semi-pervious stream-unconfined aquifer system, J. Hydrol ., 19, 43-52. Olver, F.W.J., 1974. Introduction to asymptotics and special function, Academic Press Onyejekwe, O.O., Karama, A.B, Kuwornoo, D.K., 1999. A modified boundary integral solution of recharging and dewatering of an unconfined homogeneous aquifer,Water SA, 25(1), 1-13. Parlange, J-Y., Stagnitti, F., Starr J.L., Braddock R.D. 1984. Free surface flow in porous media and periodic solution of the shallow flow approximation. J. Hydrol., 70, 251–263. Scheid, J.F., Schotting, R.J., 1997. The Interface between Fresh and Salt Groundwater in Heterogeneous Aquifers: A Numerical Approach, CWI report, MAS-R9735.

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

28

Schmitz, G.H., Edenhofer,J., 2000. Exact closed-form solution of the two-dimensional Laplace equation for steady groundwater flow with nonlinearized free surface boundary condition, Water Resour. Res., 36(7), 1975-1980. Schorghofer, N., Jensen, B., Kudroll, A., Rothman, D.H., 2004. Spontaneous channelization in permeable ground: theory, experiment, and observation, J. Fluid Mech., 503, 357-374. Serrano, S.E., 1995. Analytical solutions of the nonlinear groundwater flow equation in unconfined aquifers and the effect of heterogeneity, Water Resour. Res., 31(11), 2733-2742. Szymanski, A., 1998. Analytical solution of the regional groundwater flow problem in unconfined aquifers, Groundwater Quality 1998: Remediation and Protection, Posters, Tübinger Geowissenschaftliche Arbeiten (TGA), Reihe C: Hydro-,Ingenieur-und Umweltgeologie, Eds. M. Herbert und K. Kovar, TGA, C36, 112-114. Szymanski, A., 1999. Comment on „Analytical solutions of the nonlinear groundwater flow equation in unconfined aquifers and the effect of heterogeneity ” by Sergio E. Serrano, Water Resour. Res., 35(1), 341-343. Taigbenu, A., E., 2003. A time-dependent Green’s function-based model for streamunconfined aquifer flows, Water SA, 29(3), 257-265. Weber, J., 2003. Analytische und numerische Untersuchungen zum nichtlinearen Problem der Grundwasserströmung mit freier Oberfläche, PhD thesis, Zentrum Mathematik der Technischen Universität München. Wolfram, S., 1988. Mathematica: A System for Doing Mathematics by Computer, Addison-Wesley.

8. Appendix A “The basic equations of EBA” Following

Baiocchi and Capelo

1984

we introduce here the notion of EBA by means of the

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

29

kinematic boundary condition. ∂Φ/∂z = ( ∂Φ/∂z)2 + (∂Φ/∂z)2

z = h(x)

(A1)

where Φ(x,z) represents the piezometric head function. It is known that the seepage velocity vector at the free-surface elevation can be defined by Darcy law as follows, → vs = - K grad (Φ),

kxz = kzx = 0, kxx = kzz = k > 0

(A2)

Moreover, the kinematic boundary condition (A1) can be rewritten in the form, ∂Φ/∂n = 0

z = h(x)

(A3)

where n describes the normal direction to the h-curve. Further, at z =h(x) we have, → n := [-dh/dx (1 + (dh/dx)2)-1/2, (1 + (dh/dx)2)-1/2]

(A4)

→ τ := [(1 + (dh/dx)2)-1/2, dh/dx (1 + (dh/dx)2)-1/2]

(A5)

→ → where n , τ denote the unit normal and tangent vectors to the h-curve, respectively. Thus ∂Φ/∂n = [ - ∂Φ/∂x (dh/dx) + ∂Φ/∂z ] (1 + (dh/dx)2)-1/2

z = h(x)

(A6)

∂Φ/∂τ = [ ∂Φ/∂x + ∂Φ/∂z (dh/dx)] (1 + (dh/dx)2)-1/2

z = h(x)

(A7)

∂Φ/∂x = ∂Φ/∂τ (1 + (dh/dx)2)-1/2

z = h(x)

(A8)

∂Φ/∂z = ∂Φ/∂τ (dh/dx) (1 + (dh/dx)2)-1/2

z = h(x)

(A9)

Using (A6,7) and (A3) we obtain,

By means of (A8,9) the relation (A2) can be rewritten as follows, → → vs = - k ∂Φ/∂τ τ

z = h(x)

(A10)

or using the well-known relation dh/dx = ∂Φ/∂x (1- ∂Φ/∂z )-1

(A11)

in the following form, vx (x,z) = - k (dh/dx) (1 + (dh/dx)2)-1

z = h(x)

vz(x,z) = -k (dh/dx)2 (1 + (dh/dx)2)-1

z = h(x)

(A12) (A13)

→ Further, since the direction of the tangent vector τ is the same as the flow direction, then dh/dx

30

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

< 0, and the function h(x) belongs to the class of strictly decreasing functions.

9. Appendix B “The problems related to the parameter μ” Using the fundamental theorem of calculus from (15) one obtains, dh(x)/dx = 0.5 [k/vx(x,h(x))] {-1 + μ [1 – 4 (vx(x,h(x))/k)2]1/2}

0≤x ≤L

(B1)

where vx(x,h(x)) > 0 for x ∈ [0, L], dh(x)/dx < 0 (see Appendix A). Further, assuming that μ = 1 from (B1) formally follows; d2h(x)/dx2 = - 0.5 dvx(x,h(x))/dx {f1[vx(x,h(x))/k]/vx2(x,h(x))}

0≤x ≤L

(B2)

where f1[vx(x,h(x))/k] = 4[vx(x,h(x))/k][1 – 4(vx(x,h(x))/k)2]-1/2+{-1 + [1 – 4 (vx(x,h(x))/k)2]1/2} and 0 < vx(x,h(x))/k < 0.5 (see (15)). The graph of the function f 1[vx(x,h(x))/k] is presented in the Figure B1.

Figure B1: Function f1 [vx(x,h(x))/k] Assuming that dvx(x,h(x))/dx > 0 (i.e., a strictly increasing function) from (B2) formally follows that dh2(x)/dx2 < 0 for x ∈ [0, L]. Thus the graph of the free-surface elevation (15) for μ = 1 is convex upwards. Let us consider the case μ = -1. The relation (B2) can be rewritten in the following way,

31

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

d2h(x)/dx2 = 0.5 dvx(x,h(x))/dx {f2[vx(x,h(x))/k]/vx2(x,h(x))}

0≤x ≤L

(B3)

where f2[vx(x,h(x))/k] = 4[vx(x,h(x))/k][1 – 4 (vx(x,h(x))/k)2]-1/2+ 1 + [1 – 4 (vx(x,h(x))/k)2]1/2 and 0 < vx(x,h(x))/k < 0.5 (see (15)). The graph of the function f2[vx(x,h(x))/k] is presented in the Figure B2. Assuming that dvx(x,h(x))/dx > 0 from (B3) formally follows that dh2(x)/dx2 > 0 for x ∈ [0, L]. Thus the graph of the free-surface elevation (15) for μ = -1 is convex downwards.

Figure B2: Function f2[vx(x,h(x))/k]

10. Appendix C “The numerical examples of EBA for α → 1” We

present the numerical examples showing the behavior of EBA for α → 1. The parameters used in this case are as follows; α = 0.95 and m = -0.05, -0.2, -0.5, -0.8, respectively. The results obtained are compared with the SBA (i.e., the linearized version of EBA presented in section 4.3), and the EBA, where the horizontal velocity component has been defined as the linear function. It means that instead of the function given by (13) we simply used the linear one. We denoted this model as the EBALIN.

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

32

For the presented example LSS model (i.e., α = 1, see (50a)) gives vx(x,h(x))/k = 0.04988, where x ∈ [0, L]. Note, that even for m = -0.05 the SBA model overestimates the horizontal velocity component.

From the LSS model (see (50b)) we obtain vz(x,h(x))/k = -0.00249, where x ∈ [0, L].

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

33

The LSS model gives vx(x,h(x))/k = 0.1923, where x ∈ [0, L].

From the LSS model we obtain vz(x,h(x))/k = - 0.03846, where x ∈ [0, L].

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

34

The LSS model gives vx(x,h(x))/k = 0.4, where x ∈ [0, L].

From the LSS model we obtain vz(x,h(x))/k = - 0.2, where x ∈ [0, L].

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

35

The LSS model gives vx(x,h(x))/k = 0.4878, where x ∈ [0, L].

From the LSS model we obtain vz(x,h(x))/k = - 0.3902, where x ∈ [0, L]. Note in passing that for α → 1 the model EBALIN is a very well approximation to the EBA presented in section 4. In addition, the models EBA, SBA and EBALIN give nearly the same results for reproducing the geometrical location of the freesurface elevation. Unfortunately, the SBA does not reasonable approximate the horizontal velocity component at the freesurface. In practice, the use of

On a generalization of the Boussinesq approach for modeling the flow in unconfined aquifers

36

SBA for calculating the horizontal velocity component can only be suggested for very small values of the regional gradient. Formally, in this case the regional gradient should be treated as an arbitrarily small quantity approaching zero as a limit.