Studies on groundwater recharge through surface drains

A mathematical model is developed for the estimation of groundwater recharge through surface drains for both free flow as well as detained flow conditions. The Dupuit-Forchheimmer equation is solved using the Crank-Nicolson central finite difference scheme, to obtain the mound height matrix. The Gaussian elimination method was used to solve the matrix, to obtain the mound height at different radial distances across the drain. Various hydrogeological parameters like hydraulic conductivity, specific yield, etc. are determined by field investigations. Surface runoff available due to a particular rainfall event is correlated with recharge rate available in the drain. The model gives volume of water recharged for various rainfall events under different antecedent moisture conditions for both free flow and detained flow conditions. The value of recharge rate computed by using the model for a particular depth of flow in the drain is matched with the observed values. The model is more sensitive to change in the value of specific yield than hydraulic conductivity.


List of symbols a
= area of cross-section of flow at any z (m) A x = 1 for zone (1) and 0 for zone (2) b = average bed width of drain (m) drl = total drain length (km) = 12 km D = L/2 e = specific yield of the aquifer h = height of the water table above the base of the aquifer after the incidence of recharge.h is termed as h 1 in Zone 1 and h 2 in Zone 2 h o = depth of water table above impermeable layer i = 0,1,2…….n for i∆x ≤ L IR z = infiltration rate along the drain bed (m/s) j = 0,1, 2……n for j∆t ≤ 24 hr K = hydraulic conductivity of the aquifer L = radius of influence along radial distance from drain, at which recharge volume is assumed to meet initial water table n d = Manning's roughness coefficient which depends upon surface roughness of drain = 0.03 p(z,t) = time varying recharge rate(m 3 /s), defined as recharge volume rate per unit area receiving it q z,t = net discharge in drain at any longitudinal distance z and at any time instant t (m 3 /s) = q p -SL z = q p -IR z *(b + 2√2y(z,t))*z q p = peak flow rate for longitudinal distance z (m 3 /s) q o(m) = discharge (m 3 /s) as overflow from m th pool into (m+1) th pool q(z,t) = area under the mound curve at any longitudinal distance z and time t Q(t) = volume recharged along the drain at any time t (m 3 ) r = Δt/(Δx) 2 R = hydraulic radius at any z (m) s = height of groundwater mound above the initial water table after the incidence of recharge.

Introduction
The over-exploitation of groundwater resources is posing a serious problem of declining water tables, which have been declining at an average rate of 0.23 m/a during the past 15 years almost all over the Punjab State, India (Gupta et al., 1995).This situation of a declining water table can be handled by one of the methods of artificial groundwater recharge.It has been assessed that 0.5 m ha-m surplus runoff is available for recharge in Punjab State.The 7 000 km length of surface drains available make the drains an attractive tool for recharging this surplus runoff, which usually goes waste during monsoons.Therefore a mathematical model simulating the hydraulics of groundwater recharge through surface drains and estimating groundwater recharge was developed.
In this model the unconfined flow from drain surface to groundwater is described by the non-linear Dupuit-Forchheimer equation with the transient recharge rate.The transient recharge rate is correlated with the depth of flow and discharge rate in the drain caused by runoff coming into the drain.The inflow of runoff into the drain is taken as spatially varied flow and is described by the Saint Venant equation.
In this paper the finite difference technique has been used to study the hydraulics of groundwater recharge through a surface drain considering recharge as transient for a finite aquifer.

Mathematical formulation
A diagrammatic representation of the flow system is shown in Fig 1 .A finite aquifer with horizontal impermeable base is considered with initial water table at depth h o above the impervious base of the aquifer.The groundwater movement in the aquifer is represented by the following non-linear Dupuit-Forchheimmer equation: where: K is the hydraulic conductivity of the aquifer h is the height of mound above the impervious base x is the radial distance across the centre of the drain e is the specific yield of the aquifer p(z, t) is time varying recharge rate(m 3 /s) A x = 1 for Zone (1) and 0 for Zone (2) t is time elapsed since start of incidence of recharge subject to the following boundary and initial conditions: where: h 1 and h 2 are the mound height in Zone (1) and Zone (2) respectively.
The centre of the drain is treated as a divide, a boundary across which no flow takes place.Therefore, the centre of the drain (x = 0) is considered as an impermeable boundary.Also because of the symmetry only one half of the groundwater ridge is considered.The governing equation reduces to: while the boundary and initial conditions remain the same.
Furthermore, the governing equation is linearised by using the following transformation, where s (mound height above initial water table) is assumed to be a small perturbation of h: The boundary and initial conditions change as: where: s 1 and s 2 are the mound height above initial water table in Zone (1) and Zone (2) respectively.
Substituting Eqs. ( 7) and (8) in Eq. ( 5) and after rearranging we get the terms: where: r = Δt/(Δx) 2 and α = h o /e The last expression together with the initial and boundary con- ditions (6a-e) provide a set of simultaneous equations giving values of mound height above initial water table (s) at different radial distance (x) and time (t) that can be solved by the Gaussian elimination method (Maron, 1987).

Recharge rate
The recharge rate is correlated with runoff volume rate coming into trapezoidal section of the drain as: where: z is distance along the length of the drain in meters, b is the bed width in meters y(z,t) is depth of flow in drain at any time t in meters q z,t = net discharge in drain at any distance z and at any time t considering seepage losses in m 3 /s.q z,t = q p -IR z *(b + 2√2y(z,t))*z (11) where: q p = peak flow rate for distance z (m 3 /s) and IR z = infiltration rate along the drain bed in m/s.

Flow rate (free flow)
A trapezoidal section of drain with side slopes 1:1 is considered.
The flow to drain is classified as spatially varied unsteady flow.Basic continuity equation, known as Saint Venant equation, which describes this type of flow is used to solve for y (Subramanya, 1996).

Computation of initial depth of flow
Initial depth of flow (y o ) at any longitudinal distance z is computed by using Manning's discharge equation as: The peak runoff rate coming into the drain is computed by curve number method (Murty, 1998).Substituting values of R (hydraulic radius) and a (area of cross-section of drain) in Eq. ( 15) and rearranging we get: where: n d is Manning's roughness coefficient which depends upon surface roughness of drain s o is slope of drain (m/m).
From this equation y o is calculated by bisection method (Grewal, 1998).

Flow rate (detained flow)
At the start of the incidence of recharge, i.e. when flow rate is maximum total discharge and m th pool consists of runoff contribution from watershed area between m th pool and (m-1) th pool and overflow from (m-1) th pool, if any.The overflow is measured by using a 90º V-notch weir.
It is assumed that with check detentions the water is stored as pool and there is no channel flow movement.Depth of water only changes due to infiltration losses.Thus the depth of flow is related to time by the relation: The initial depth of flow equation is obtained as: [(b+y m )y m ] 5 /(b+2√2y m ) 2 = [q o(m-1) + q p(m) ] 3 n d 3 / s o 3/2 (18) where: q o(m-1) is the discharge (m 3 /s) as overflow from (m-1) th pool into m th pool y (m) is the depth of flow(m) in m th pool q p(m) is the peak runoff rate (m 3 /s) coming into the drain at m th check.

Total volume recharged
At any cross-section along the length of drain, the area under the mound curve represents the amount of water recharged.Using Darcy's law the area under the s -curve at any distance z from start of drain and at any time t is written as: where: q(z,t) is area under the mound curve at any distance z along the drain and time t L is radius of influence (m) from centre of drain after which recharge mound meets the initial water table.
Using Simpson's rule of numerical integration Eq ( 19) is written as: where: s (0,t) , s (1,t) …..s (2n,t) are the mound height values at centre of the drain, 1 st node up to 2n th node respectively at z distance along the drain length and time t.
To get the total volume, q(z,t) is integrated over the whole length of drain, drl, as: Various values of parameters are determined by field investigations.A computer program is encoded in C++ language to solve for transient recharge rate available, mound profile and total volume of water recharged through the drain at any time t from the start of the incidence of recharge.A flow chart of the program is shown in Fig. 2.

Results and discussion
For the field investigations, the Raipur link drain in Ludhiana district of Punjab State, India is selected.The length of the drain for computational purposes is taken as 12.0 km.From the drain data the weighted bed width, weighted slope, check spacing and height of checks comes out to be 8.33 m, 0.000314 m/m, 3 km and 1.02 m respectively.Manning's roughness coefficient for drain is taken as 0.03.The values of various parameters are measured by field investigations and are given Table 1.
In the finite difference schemes used, the time step is taken as 4 min for free flow conditions and 1 h for detained flow conditions.Longitudinal distance step (Δz) for free flow conditions and detained flow conditions is taken as 4 km and 3 km respectively.Radial distance step (Δx) is taken as 16 m for both the flow conditions.
Rainfall value is given as input to the model to know for peak discharge rate and volume of runoff coming into the drain.This discharge rate is used to calculate variation of the depth of flow along the drain length with time, thus correlating it to the transient recharge rate term in the Dupuit-Forchheimer (D-F) equation.The finite difference scheme of the D-F equation is solved by using the Gaussian elimination method (Maron, 1987) to get mound profile at any drain length and at any time instant t.Darcy's law is employed over mound profile, using Simpson's numerical integration rule to give the volume of water recharged.Based on the model the results obtained are discussed below.

Recharge rate
Recharge rate available is computed at different depths of flow in the drain.The computed values are compared with the values as observed during the study conducted under ICAR National Professorship Project, Punjab Agricultural University, Ludhiana, India (Khepar, 2001).The comparison of the values is shown in Table 2.
Results of the paired t-test on the above values showed that there is no significant difference between computed and observed values.The integral is multiplied by two as we are considering only one half of the drain due to symmetry.For total volume q(z,t) has to be doubled.The integral is computed using Simpson's rule as: Q(t) = 2 ∆z/3[q(0,t) + 4{q(1,t) + q(3,t) +….+q(2n-1,t)} + 2{ q(2,t) + q(4,t) + ….+ q(2n-2,t)} + q(2n,t)] The average rise in water table is computed as: the flow velocity is quite low as it encounters resistance due to weeds and depressions in the drain.Thus the conditions become more or less similar to detained flow conditions.

Threshold rainfall
Test runs are made on the model to investigate the threshold rainfall, i.e. the minimum amount of rainfall which can cause significant recharge under different antecedent moisture conditions, viz.AMC-I, AMC-II and AMC-III (Murty, 1998).Based on the numerical simulations it is found that a minimum of 32 mm rainfall for AMC-I, 17 mm for AMC-II and 8 mm for AMC-III is necessary to cause any significant recharge.

Sensitivity analysis
A sensitivity analysis of the computed results was performed and checked for its sensitivity by changing the value of the hydraulic conductivity (K) and specific yield (e) by ±30%.The model is found to be more sensitive to change in specific yield than change in hydraulic conductivity as shown in Table 4.

Conclusions
A mathematical model for hydraulics of groundwater recharge through surface drains was developed.It was applied to a drain for the estimation of groundwater recharge under varying antecedent moisture conditions and different flow conditions viz.free flow and detained flow.It is concluded that the recharge rate increases with increase in depth of flow in the drain and the volume of water recharged under detained flow conditions varied between 1 and 9 times that under free flow conditions.

Volume recharged
Volume of water recharged for both free flow conditions and detained flow conditions is computed for selected rainfall events, i.e. 32 mm, 66 mm and 100 mm rainfall.It is observed that the volume of water recharged over 24 h after peak discharge rate of runoff for detained flow conditions varies between 1.2 and 9.4 times that for free flow conditions as shown in Table 3.The ratio of water recharged in detained flow conditions and in free flow conditions is lower for low rainfalls while it is greater for high runoff-causing rains.This is because for high runoff causing rains, the flow is routed through the drain in case of free flow conditions whereas it is detained almost completely in the case of detained flow conditions.In case of low runoff-causing rains across the centre of drain (m) y(z,t) = depth of flow in drain at any time t (m) y o = initial depth of flow y m,j = depth of flow at m th pool and at time node j (m) z = longitudinal distance along drain (m) ∆x = radial distance step ∆z = longitudinal distance step ∆t = time step i.e. time difference between two nodes along t-axis α = h o /e

Figure 1
Figure 1 Diagrammatical representation of groundwater mound

Figure 2
Figure 2Flow chart of the groundwater recharge model for surface drainage system