www.elsevier.com/locate/joes

Available online at www.sciencedirect.com

ScienceDirect

Journal of Ocean Engineering and Science 000 (2016) 1-10

Original Article

Approximation of surface-groundwater interaction mediated by vertical

streambank in sloping terrains

Rajeev K. Bansal

Department of Mathematics, National Defence Academy, Khadakwasla, Pune 411023, India

Received 15 April 2016; received in revised form 20 September 2016; accepted 5 October 2016

Available online xxx

Abstract

New analytical solutions are derived to estimate the interaction of surface and groundwater in a stream-aquifer system. The analytical model consists of an unconfined sloping aquifer of semi-infinite extant, interacting with a stream of varying water level in the presence of a thin vertical sedimentary layer of lesser hydraulic conductivity. Flow of subsurface seepage is characterized by a nonlinear Boussinesq equation subjected to mixed boundary conditions, including a nonlinear Cauchy boundary condition to approximate the flow through the sedimentary layer. Closed form analytical expressions for water head, discharge rate and volumetric exchange are derived by solving the linearized Boussinesq equation using Laplace transform technique. Asymptotic cases such as zero slope, absence of vertical clogging layer and abrupt change in stream-stage can be derived from the main results by taming one or more parameters. Analytical solutions of the linearized Boussinesq equation are compared with numerical solution of corresponding nonlinear equation to assess the validity of the linearization. Advantages of using a nonlinear Robin boundary condition, and combined effects of aquifer parameters on the bank storage characteristic of the aquifer are illustrated with a numerical example. ©2016 Shanghai Jiaotong University. Published by Elsevier B.V.

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Keywords: Stream-aquifer interaction; Streambank; Sloping aquifer; Boussinesq equation; Laplace transform.

1. Introduction

Mathematical formulation of subsurface seepage flow over horizontal beds is derived using conventional Dupuit-Forchheimer assumption in which streamlines are considered to be approximately horizontal and the gradient of hydraulic potential same as the absolute slope of the water table [18]. When the bed of an unconfined aquifer bed is sloping, the flow is constrained to the direction nearly parallel to the bed [10,13], i.e. streamlines are approximately parallel to the sloping bed, and equipotential are perpendicular to the bed. The consideration that the flow pattern is controlled by the bed slope yields the discharge rate, q, per unit width of the aquifer by the relation [12]:

q = -K h cos2 ß

--tan ß

E-mail address: bansal_rajeev31@hotmail.com.

where h is the height of the water table measured in the vertical direction from the sloping base; x is spatial coordinate; K is the hydraulic conductivity of the aquifer, and tan j is the bed slope. Considering a representative elementary volume and applying principle of mass balance on it, the flow of subsurface seepage in horizontal and vertical axes system is approximated by the following nonlinear Boussinesq equation:

d ( dh \ dh S dh

— h — - tan j— =------(2)

dx \ dx J dx Kcos2 j dt

where S is the specific yield of the aquifer. Eq. (2) is indeed a nonlinear advection-diffusion equation which plays key role in approximation of unconfined groundwater flow over sloping bed. Numerical and experimental studies [19] indicate that the results obtained from Eq. (2) are reasonably valid for bed slope up to 30%, i.e. tan j < 0.3. Due to non-linearity, Eq. (2) is analytically intractable. However, approximate analytical solutions of Boussinesq equation have been used in numerous studies to analyze the transient behavior of

http://dx.doi.org/10.1016/j.joes.2016.10.002

2468-0133/© 2016 Shanghai Jiaotong University. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

R.K. Bansal/ Journal of Ocean Engineering and Science 000 (2016) 1-10

Fig. 1. Schematic diagram of the stream-aquifer system consisting of an unconfined aquifer of semi-infinite extent overlying an impervious sloping base and interacting with a time-varying stream in presence of thin vertical sedimentary layer.

stream-aquifer models under varying hydrological conditions [1-4,11,15,5-7,23,8].

Most of the fundamental studies on surface-groundwater interactions focus on estimation of water exchange between streams and aquifers, without taking the effects of vertical clogging layer into account [14,16,22,24]. Often, if not always, the water interaction between streams and aquifers is mediated by a vertical layer of low permeable sedimentary deposits, henceforth referred to as streambank. Such layers are formed by deposition of silt particles or loading of waste along river/stream bank in rivers whose streamflow is sluggish [9,20,25]. These layers decelerate the seepage flow from river to aquifer and vice-versa. For horizontal strata, the approach to estimate the flow through the streambank is to consider the seepage velocity a function of difference between stream stage and aquifer's water table [17]; however, there is not much discussion in the literature on approximation of seepage flow through streambank when a fully penetrating stream interacts with an unconfined sloping aquifer.

Recently, Bansal et al. [8] derived an appropriate form of boundary condition to simulate the flow through streambank in an archetypical stream-aquifer model in which the aquifer is resting on sloping impervious bed. The problem of seepage flow between an unconfined aquifer of semi-infinite extent and a fully penetrating stream of time-varying water level in the presence of vertical streambank (see Fig. 1) is characterized by the following linearized form of Boussinesq Eq. (2).

d 2h tan p dh dx2

havg dx K havg cos2 p dt

where havg is the average saturated depth of the aquifer given by an iterative formula havg = (h , + ht )/2; h , is the initial depth and ht is the variable height at time t, at the end of which havg is calculated. Considering that the flow at stream-aquifer interface is made up of two terms, namely flow due to head variation across the streambank, and flow due to bed slope, they derived an appropriate nonlinear boundary condition at the stream-aquifer interface, given by

—Kcos2p h (x = 0+, t)( — )

\d x Jx=0

h2 (x = 0+, t ( — h2 (t )

where k and b are the hydraulic conductivity and thickness of the clogging layer. The stream-stage variation hs ( t) from an initial level h, to a final level h, controlled by a nonnegative parameter X, is given by the transient relation

ht (t ) = h( — (hf — h] (e

Moreover, an initial condition of depth independent water table reads

h (x, t = 0 ) = h,

and a far-end boundary condition is

h(x ^ to, t) = h,

Assuming that the head loss across streambank is small, Bansal et al. [8] used the approximation |h, (t) - h (x=0+, t)| << h(x=0+, t). With this consideration, they simulated the nonlinear boundary condition (4) by a simpler linear equation

—K cos2 p

h (x = 0+, t ( — ht (( )

The term (Kb cos2 p )/k = r denotes the streambed retardation coefficient or streambank leakance in sloping aquifer. This parameter controls the degree of hydraulic connection between the stream and aquifer, and depends on the physical properties of the sediments that make this layer. It is demonstrated in Bansal et al. [8] that the linearization of boundary condition (4) yields acceptable results only when the stream-stage variations and streambank leakance are in narrow range (ht-h, < 30% of h, and streambank leakance r < 10). When the ratios |h, - h, |/h, and streambank leakance are large, the analytical solutions developed by Bansal et al. [8] varied considerably from the numerical solution of the corresponding to nonlinear boundary condition (4). These facts highlight practical limitations of the results developed by Bansal et al. [8].

The aim of the present study is to develop new analytical solution of the linearized Boussinesq equation subjected to nonlinear boundary condition (4) along with initial condition (6) and a far-end boundary condition (7). Closed form analytical expressions are obtained using Laplace transform method to predict the spatio-temporal distribution of water head, rate of discharge and volumetric exchange between stream and aquifer for rise and decline in the stream-stage. Few previously known results are derived as special cases of the main results of this study. Numerical solution of the nonlinear Boussinesq Eq. (2) subjected to nonlinear condition (4) are also developed using an explicit Dufort-Frankel method to assess the validity of analytical solutions for large value of the ratios |h, - h, \/h, and streambank leakance ratio r,

2. Analytical solution of linearized Boussinesq equation

Assuming that the unconfined downward sloping aquifer is homogeneous and isotropic, and the streamlines are nearly parallel to the sloping bed, the flow of subsurface seepage is governed by nonlinear Boussinesq Eq. (2). This equation is a second order parabolic partial differential equation whose

R.K. Bansal/ Journal of Ocean Engineering and Science 000 (2016) 1-10

analytical solution is not tractable; however closed form solutions of the corresponding linearized Boussinesq equation are widely accepted for small variations in the stream-stage. It is worth mentioning that the Boussinesq equation can be linearized in terms of h or h2 , depending on the type of boundary conditions used in the study. Moreover, it has been brought out in several studies that linearization in terms of h2 yields better approximation of the solution than that of h. In this study, linearization in terms of h2 is preferred due to nonlinear nature of boundary condition (4). Rewrite Eq. (2) as

1 d h2

d /1 d h2 \ / 1 d h2 \ _ dx \2 ~dx J - tan 2h ~dx J = Kcos2 j V 2h dt

Now, replacing the term 1/h associated with second and third bracket by some judicious depth h avg which can be obtained using analogy of Marino [16] with the formula havg = (h i + ht)/2, where hi is the initial water table and ht is the water table at time t at the end of which, havg is calculated. Eq. (9) reads:

d (dh2 \ tan ß (dh

d X \ d X

Khavgcos2 ß \ dt

Water table in the stream varies (rise and decline) from an initial level hi to a final level hf by a known exponential function of time t given by Eq. (5). The depth independent initial condition is represented by Eq. (6), whereas the variations in the boundaries are simulated using conditions (4) and (7). Eq. (10) is linearized form of Boussinesq equation; primarily an advection-diffusion equation; which can be solved using Laplace transform technique. The main steps in the solution procedure include introduction of the following dimensionless variables h 2 h2

h2 - h2 '

x = —,

K havg COS2 ß

so that Eq. (10) is transformed to the form

d 2H dH _ dH

dxx2 - a~dx = ~3T

h i tan j a =-

The initial and boundary conditions become H (X, 0 ) = 0

H(X ^ œ, T) = 0

= H (X = 0+, t)

- 1 + p e-k12 + (1 - p)e-2À1 2

r K b cos2 ß

R = - =--; A.1 = 0

hi k h i K havgcos2 ß

A; p =

h f + hi

Now, Laplace transform is applied on Eq. (10), leading to the following ordinary differential equation

d2H dH

—T - 2a--sH = 0

dX 2 dX

where H denotes the Laplace transform of Hi Eq. (18) can be solved using standard method. The arbitrary constants contained in its general solution can be determined by invoking Laplace transform of boundary conditions (15) and (16). The resulting solution is

H (X, s ) =

(a-^ÖTSy - 1 Is + s + 2A1

(a-*J a2 +s )X

Using various properties of inverse Laplace transform, solution of Eq. (19) can be obtained in terms of variable t . After some simplification, the solution of Eq. (12) reads:

H (X, t)=2 {erfc ( - ^

e2aX ( X +——erfc

F ""'yii/r p e-klt [ e(a -a1 ) X

+ a^Tj

2 [1 - R (a - a\ ) \2„Jr

e (a+ai ) X ( X 2

+ --—:-r erfc ( ——= + ax „Jt

- a\„/r

1 - R(a + ai ) \2„Jr

( 1 - p) e-2A1 T

e (a-a' 1 ) X

1 - R (a - a' 1 ) \2„/r

- a 1 Vt

e (a+a' 1 ( X

+--erfc

1 - R (a + a' 1 ) p

+ ( 1-aR )

2 „Jr 1 - p

+ a' 1 „/r

F + A1R2 F + 2A1 R2 F

, X 1 - aR

x erfc | +--=— \l t

R + F „)

where ax = (a2-Xi)1/2 , ai' = (a2 -2A.i)1/2 , and F = 1 -2a R. It can be observed that for marginally sloping aquifers and rapidly rising streams, a2 < implying that the parameter a1 is pure imaginary. However, the complex numbers corresponding to complementary error function appear in conjugate pairs, and thus, the resultant value of the right-hand side is always real. For upward sloping aquifers, the corresponding results can be obtained by replacing a by -a. When the impervious bed is horizontal aquifer, we have a = 0, a1 = i (X1 )1/2 , a1' = i (2X1 )1/2 and F = 1. Using them in Eq. (20), the result can be expressed as

R.K. Bansalf Journal of Ocean Engineering and Science 000 (2016) 1-10

H (X, t) = erfc +

2 —fx P

x erfc

1 + X\R( X

1 — P

1 — 2X1R2

- 1 \e (R+R()

2 (/r + R (

—i (/Ï1X

2 [1 + i R (/11 '

i —x X / x

1 — i R (/Ti (1 — p) e—2XT

0—i (/2x1X

1 + i R (/2X1 \2 —/r

— i (/2X

e i (/2X1 X

1 — i R (/2X1 \2 —/r

+ i (/2X11

which are in complete agreement with results of Teloglou and Bansal [21] when the downward leakage in their study is considered to be zero. For instantaneous rise in the stream-stage (A.i ^ to), Eq. (21) reads

H (X, t) = erfc

X 2 —/r

-e (R + R( ( /

--1---

2 (fx R

which are identical to the expression for water head distribution derived by Marino [16].

3. Analytical expressions for discharge rate and volumetric exchange

Discharge rate per unit width of the aquifer is given by Eq. (1). The mass balance equation that leads to simulation of groundwater flow by Eq. (9) is based on the assumption that the linear expression for discharge rate is given by

K cos2 p

d h2 d x

Now, define a normalized discharge rate, Q, which can be expressed in terms of dimensionless water head, H, as

-q = —

--2a H — C

K(h2 - h2,cos2 P

where Q = Q(X, t ), q = q(x, t), H = H(X, t ), C = (h,3 P )/{havg (h, - h, 2 )}. The dimensionless discharge rate through unit width of streambank is given by

q0 (t) = Q(X = 0, t)

— — 2aH (X = 0+, x) — C]

d X /X =0

Evaluation of right-hand side requires use of Eq. (20), we obtain

Q0 (t) = —-H (X = 0+ , t) R

1 — p e—XT — ( 1 — p)e—2XT

+ R + C

The value of H(X = 0+, t ) can be obtained from Eq. (20). Noting that erfc(-z) = 2 - erfc(z), the expression for the flow rate can be simplified as

, «1 p F , t _

q0 (t ) = -a erfc[a^T] + e A1 t erfc(a1 Vt)

F + X1R2

a' 1 ( 1 — p)F T , —

+ —-—,r e—2X1 T erfc(a' 1Afx)

F + 2X1 R2 J v

( 1 —aR )F ~R

(1 — P)

F + X1 R2 F + 2X1 R2

* /1 — aR

x erfc I —-— Vt

+ 2a + C — ( 1 — P)

1 — R(a — a1 ) F

1 — R(a — a' 1 )

-2x1 t

In case the aquifer is underlain by a horizontal impervious base, the expression for flow rate reduces to

Q0 (t) = 1 2 e-"1 t erfc(i VX1T)

( 1 — P)

- 1 \ e R(

R [1 + X1 R2 1 + 2 X1 R2

x erfc (R —R {1—mbXr

—x1 t

( 1 — P) R

1 + R (/2X7

-2x1 t

Furthermore, if the variation in the stream stage is instantaneous, the flow rate at stream-aquifer interface can be obtained by setting X1 ^ to in Eq. (28). The resulting expression is

Q0 (x) = R e F t erfc VT^j

which is in full conformity with Marino [16].

Net volume, v(t), of water that enters in or leaves the aquifer through unit width of stream-aquifer interface can be obtained by integrating the flow rate q(x = 0, t) with respect to time t in the interval [0, t]. In terms of normalized variable, the dimensionless volumetric exchange can be defined as:

V (T) =

S h( (hf — h( )

= / Q0 ( t ' ) d t 1

where t' is the variable of integration. Using various properties of complementary error function and its product with

e r2 (

—x1 t

R.K. Bansal/ Journal of Ocean Engineering and Science 000 (2016) 1-10

[m5+;November 29, 2016;9:48] 5

exponential function, we finally obtain

hm „+1 . That i

V (?) = J-e n

a12 pF

a [ 2 A1 (F + k1R2 )

+ + ( 1 - a R )2 Q}erf (a—?)

a1 p F e-A 1 T ,—

ar erfc(as/r) + ^ (F + Alr2) erfc(a1 s/?)

a' 1 ( 1 - p) F e-2A 1 T , _

+--~—erfc(a 1 Jt)

2A1 (F + 2A1R2 ) 7 V

T ( 1 - a R

+ R( 1 - a R)Q eR„„ erfc„ —=— V?

+ (2a + C )t + Q'

- P- (1__F_

R \ 1 - R(a - a1 ) ( 1 - p) ( F

1 - e-

where Q =

1 - R(a - a' 1 ) 1 - p 1

-2a1 t

F + A1R2 F + 2A1R2 F

a' 1 ( 1 - p)F

A1 (F + A^2 ) 2A4 (F + 2A1Ti2 ) 4. Comparison of analytical and numerical solutions

Numerical solution of the nonlinear Boussinesq Eq. (2) subjected to conditions (4)-(7) is obtained by employing the Du Fort and Frankel scheme. The proposed scheme is an explicit finite difference numerical scheme which proceeds in three time levels. Rewrite Eq. (2) as

— = A1 dt

h d 2 h (d h dX„ V dX

- a2—

where A1 = (K cos2 j)/S and A2 = K sin2j/(2S). In Du Fort and Frankel scheme, the first order derivatives of h with respect to x and t are approximated using central differences, i.e.

d hn hn — hn dJhn = hm+1 hm-1 + ot( AX )2 „ dX 2 Ax „ ) '

d h" "" m

h n+1 h n-1

+ O( ( AX )2 )

whereas, the second order derivative d2h/dx2 is approximated by

d 2 hn hn — 2hn -I- hn

d hm _ hm+1 2hm + hm-1

( AX )2

+ O( ( Ax )2 )

Furthermore, the middle term hm n in the right-hand side of Eq. (37) is replaced by the arithmetic mean of hmn-1 and

d 2hm _hm+1 - (hm+1+hm-1)+hm-1

(AX )2

where m = 1 and n = 1 correspond to x = 0 and t = 0 respectively. The above finite difference approximations are second order accurate. Using them in Eq. (34), we get

h n+1 h n-1

h n h n+1 h n-1 h hn I hm+1 hm hm + h-.

( AX )2

h„ — hn

m+1 m-1

"m+1 ' lm-1

Expressing hn+1 explicitly in terms of hn and hn 1 , we obtain

2 At (Ax)2

h"m+1 = A1

hm+1 - hm-1+h

hn — hn

m+1 - m-1

hn — hn

m+1 - m-1

It is worth noting that Eq. (40) is two-step method requiring values of h at time steps n-1 and n to determine the value at current step n + 1. Consequently, the above formula cannot be used for determination of water head at time step n = 2. Thus, for n = 2, we use an alternative discretization of Eq. (34), which consist backward difference approximation for dh/d t i That is

hnm+\ - 2hm-1 + h

hn-1 m

n-1 m-1

(AX )2

"m+1 ' m-1

' hn-1 - hn-1 m+1 - m-1

Setting n = 2 in it, we get hi = hi + A1 At

h 1 ihm+1 2hm+hm -1

hm V ( Ax )2

h1 h 1

"m+1 m-1

- A2 At

"m+1 m-1

The initial and boundary conditions are discretized as follows:

hl = h „

hn h nM —

R.K. Bansalf Journal of Ocean Engineering and Science 000 (2016) 1-10

-t=1 -t=2 t=5 -t=10 -t=15 t=30

7 6 5 4

2 1 0 -1

Fig. 2. RPD between analytical and numerical solution under nonlinear boundary condition. The rise in the stream-stage is: (a) 100% of h, , and (b) 200% of h, .

where M is a large number such that the product MAx can be treated as infinity for computational purpose. To maintain second order accuracy in the discretization of nonlinear boundary condition (4), we use the following central difference approximation

—3h + 4hn2 — h3 2Ax

+ O (( Ax )2

so that, Eq. (4) reads

—3h1 + 4h2 — h3 \ (h" (2 — h2 (t )

where R' = (K b cos2 p)/k, Eq. (46) can be further simplified to obtain the values of water head at stream-aquifer interface by solving the quadratic equation

R+A )2—

4h2 — h3 Ax

h1 —

hf (t ) R'

Numerical experiments indicate that the method is stable when At,(Ax)2 < (0.05)2 . Analytical solutions are compared

with numerical solution of the nonlinear Boussinesq equation using L2 and Tchebycheff norms. The L2 norm determines the average distance between the numerical and analytical solutions using the formula

v—i2 ihfum hana I — -s/L

{hnum (x 2 2 — ana (x, t)}2 dx

where hnum (x, t) and hana (x, t) respectively denote the numerical and analytical solutions. L is the range of spatial coordinate in which the L2 norm is calculated. The maximum difference between the analytical and numerical solutions in the domain 0 < x < L is determined using Tchebycheff norm, defined as follows:

hana I I — max Ihnum hana i

0 < x < L

R.K. Bansal/ Journal of Ocean Engineering and Science 000 (2016) 1-10 7

2 10 a

6 4 2 0

-t=1 -t=2 t=5 -t=10 -t=15 t=30

75 x (m)

Fig. 3. RPD between analytical and numerical solution under linear boundary condition. Rise in the stream-stage is: (a) 100% of hi , and (b) 200% of hi . Table 1

Numerical values of RPD between analytical and numerical solutions under nonlinear boundary condition. Values of t are in days.

RPD for 100% rise

RPD for 200% rise

t = 2 t = 5 t = 10 t = 15 t = 30 t = 2 t = 5 t = 10 t = 15 t = 30

0 -0.02 -0.04 -0.06 -0.09 -0.12 -0.06 -0.06 -0.08 -0.13 -0.19

25 0.06 0.54 0.59 0.26 -0.23 0.26 1.87 1.28 0.43 -0.41

50 0 0.21 1.17 1.17 0.08 0.01 1.07 3.97 2.68 0.01

75 0 0.02 0.64 1.50 0.81 0 0.11 3.25 5.38 1.35

100 0 0 0.15 0.87 1.60 0 0.01 0.9 4.52 3.77

125 0 0 0.02 0.29 1.86 0 0 0.13 1.77 6.24

150 0 0 0 0.07 1.36 0 0 0.01 0.41 6.47

5. Discussions of results

Applications of the new solutions developed in this study are illustrate with the help of a numerical example

with aquifer parameters as: K = 2.5 m/day, S = 0.25, j = 5°, k=0.25 m/day and b = 1 m. In order to show the effective-

ness of the solutions developed in this study, two cases of large variations in the stream-stage are considered. In the first case, the rise in stream water level is 100% of the ini-

tial level (hi = 10 m, hf=20 m, X = 0.2 day-1 ), while in the second case, the rise is 200% of the initial level (h ; = 10 m, hf=30 m, X = 0.2 day-1 ). Mean saturated depth havg of the aquifer is determined using analogy of Marino [16] with an iterative formula havg = (ht + h t )/2. Spatial and temporal distributions of water head obtained from analytical solution are compared with numerical solution and found to be in good agreement during initial stages. As time progresses, difference between these solutions grows. For quantitative assessment of

8 R.K. Bansal/ Journal of Ocean Engineering and Science 000 (2016) 1-10

Table 2

Numerical values RPD between analytical and numerical solutions under linear boundary condition. Values of t are in days.

x (m) RPD for 100% rise RPD for 200% rise

t = 2 t = 5 t = 10 t = 15 t = 30 t = 2 t = 5 t = 10 t = 15 t = 30

0 0.25 0.83 1.04 0.94 0.59 0.88 2.04 2.05 1.70 1.03 0.25

25 0.12 1.81 3.79 3.78 2.48 0.54 6.25 9.03 7.65 4.41 0.12

50 0 0.44 3.37 5.17 4.49 0.01 2.00 11.42 13.33 8.71 0

75 0 0.04 1.32 3.92 6.01 0 0.18 5.94 13.73 13.39 0

100 0 0 0.28 1.78 6.25 0 0 1.41 8.01 16.86 0

125 0 0 0.04 0.53 5.03 0 0 0.20 2.72 16.9 0

150 0 0 0 0.12 3.15 0 0 0.02 0.62 13.10 0

3.5 3 2.5 2 1.5 1

Fig. 4. Comparison of analytical and numerical solutions under linear (dotted curves) and nonlinear (continuous curves) boundary condition (a) L2 norm, and (b) Tchebycheff norm.

the difference between analytical and numerical solutions, we use Relative Percentage Difference (RPD) which is defined as follows:

hnum — hana RPD = -num-— x 100

For both the aforementioned cases, RPD is determined in the range 0 < x < 150 at t = 1, 2, 5, 10, 15 and 30 days (Fig. 2 and Table 1). To demonstrate the efficiency of the solution developed in this study, RPD is also determined for the same set of aquifer parameters using solutions of Bansal et al. [8] (Fig. 3 and Table 2).

-t=1 d -t=2 d -t = 5 d -t = 10 d -t = 15 d t = 30 d

75 x (m)

t=1 d t=2 d t = 5 d -t = 10 d -t = 15 d t = 30 d

75 x (m)

Fig. 5. Comparison of water head distribution predicted by the solutions of this study (continuous curves) and Bansal et al. [8] (dotted curves) for (a) 100% rise, and (b) 200% rise.

It is clear from these figures and tables that the overall RPD between analytical and numerical solutions of this study is much lesser than that of Bansal et al. [8], indicating higher level of accuracy of the present results as compared to Bansal et al. [8] mainly due to linearization in term of h2 instead of h] For 100% and 200% rise in the stream stage, the maximum RPD yielded by the present solution are merely 1.85% and 6.9% respectively, whereas in case of Bansal et al. [8], these values are more than 6.5% and 17.4% respectively. One clear trend that can be observed from Figs. 2 and 3 is that the peak of RPD rises and drifts toward higher values of x as time t increases. Higher accuracy of current results is also reflected

R.K. Bansal/ Journal of Ocean Engineering and Science 000 (2016) 1-10 9

■ r=10

■ r=15 r=20

■ r=25

■ r=30 r=35

Fig. 6. RPD between analytical and numerical solutions at t = 15 days for various values of r.

Fig. 7. Comparison of water head heights predicted by analytical and numerical solutions: (a) L2 norm, and (b) Tchebycheff norm.

in the plots of L2 and Tchebycheff norms presented in Fig. 4. The integral involved in the right-hand side of Eq. (48) is evaluated using trapezoidal rule in the range 0 < x < 150 m.

Distribution of water head in the unconfined sloping aquifer predicted by the solutions of this study and that of Bansal et al. [8] is plotted in Fig. 5(a) (for 100% rise) and Fig. 5(b) (for 200% rise). It can be concluded from these figures that the analytical solutions based on the linear boundary condition predict significantly lower water table than that the corresponding results with non linear boundary condition.

Another major limitation of linear boundary conditions is that it yields acceptable results only when the streambed leakance ratio, r = (K b cos2 j)/kt is small. It is highlighted in Bansal et al. [8] that when the value of r is large, the difference between analytical and numerical solutions grows significantly. In Fig. 6, the RPD between analytical and numerical solutions of this study for r = 10, 15, 20, 25, 30 and 35 at t = 15 days is presented.

It can be observed from Fig. 6 that the differences between analytical and numerical solutions are very small. Moreover, the RPD decreases as r increases, establishing the accuracy of present solutions. Average distance between solutions (L2 norm presented in Fig. 7(a)) decreases almost linearly; however, the maximum difference (Tchebycheff norm, presented in Fig. 7(b)) increases with r.

6. Conclusions

Analytical solutions developed by Bansal et al. [8] for estimation of stream-aquifer interaction in presence of vertical sedimentary layer are based on a linearized boundary condition, and exhibit significant limitations when the stream-stage variations and streambed leakance are large. To overcome these limitations, new analytical solutions are developed and tested by incorporating the effects of nonlinear boundary condition for estimation of surface-groundwater interaction. The solution methodology uses Laplace transform to derive system response functions for water head distribution, flow rate and volumetric exchange. It is shown in the study that these results are consistent with numerical solution even when the stream-stage variations are extremely large (200% of the ini-

R.K. Bansal/ Journal of Ocean Engineering and Science 000 (2016) 1-10

tial depth), and the streambed leakance ratio is high (r < 40). Water table fluctuations predicted in this study are compared with that of Bansal et al. [8] and observed that except for the initial time duration, the results predicted by the Bansal et al. [8] are significantly lower than the actual values.

References

[1] E. Akylas , A.D. Koussis , J. Hydrol. 338 (2007) 85-95.

[2] R.K. Bansal, S.K. Das , J. Hydrol. Eng. 14 (8) (2009) 832-838.

[3] R.K. Bansal , S.K. Das , J. Hydrol. Eng. 15 (11) (2010) 909-917.

[4] R.K. Bansal, S.K. Das , Water Resour. Manage 25 (2011) 893-911.

[5] R.K. Bansal, Transp. Porous Med. 92 (2) (2012) 817-836,

[6] R.K. Bansal , J. Hydraul. Eng. 139 (11) (2013) 1165-1174,

[7] R.K. Bansal, Appl. Water Sci. (2015), doi:10.1007/s13201- 015- 0290- 2.

[8] R.K. Bansal, C. Lande, A. Warke, J. Hydrol. Eng. (2016), doi:10.1061/ (ASCE)HE.1943-5584.0001362.

[9] A.L. Berestov , H.J.S. Fernando , P. Fox , Water Resour. Res. 34 (11) (1998) 3025-3033.

[10] J. Boussinesq, Mem. presents par divers savants Acad. Sci. Institut Fr. 23 (1) (1877) 252-260 (in French),

[11] W. Brutsaert, Water Resour. Res. 30 (10) (1994) 2759-2763,

[12] T.G. Chapman , Water Resour. Res. 16 (6) (1980) 1114-1118,

[13] E.C. Childs , Water Resour. Res. 7 (1971) 1256-1263,

[14] R.E. Glover, G.G. Balmer, Trans. Am. Geophys. Union 35 (1954) 468-470 .

[15] J. Mas-Pla, J. Montaner, J. Sola, J. Hydrol. 216 (1999) 197-213,

[16] M. Marino, J. Hydrol. 19 (1973) 43-52,

[17] A.F. Moench , P.M. Barlow, J. Hydrol. 230 (2000) 192-210,

[18] P.Ya. Polubarinova-Kochina, Theory of Groundwater Movement, Princeton University Press, NJ, 1962,

[19] K.N. Shukla, H.S. Chauhan, V.K. Srivastava, J. Irrig. Drain. Eng. ASCE 116 (1) (1990) 107-113,

[20] M. Sophocleous , Hydrogeol. J. 10 (2002) 52-67,

[21] I.S. Teloglou, R.K. Bansal, J. Hydrol. 428-429 (2012) 68-79,

[22] C.Y. Theis , Trans. Am. Geophys. Union 22 (3) (1941) 734-738,

[23] N.E.C. Verhoest , P.A. Troch , Water Resour. Res. 36 (3) (2000) 793-800.

[24] S.R. Workman, S.E. Serrano, K. Liberty, J. Hydrol. 200 (1997) 149-164,

[25] V.A. Zlotnik, H. Huang, Ground Water 37 (4) (1999) 599-605,