Analytical and numerical solutions to describe water table fluctuations due to canal seepage and time-varying recharge

To improve agricultural production, planners and policymakers consider canal irrigation as one of the effective methods of providing irrigation to crops growing in the command area. But most of the canals and their distribution networks are unlined and act as a source of seepage. Before planning for any intervention to control seepage loss, it is essential to have knowledge of spatial and temporal variations of the water table profile in an unconfined aquifer, which can be achieved by representing physical situation in a mathematical term and solving a problem (consisting of governing equation and initial and boundary conditions) employing analytical and numerical solutions.

Many researchers have attempted to study the effect of recharge due to canal seepage and field irrigation resulting in water table rise in the affected areas. Kraijenhoff van de Leur (1958), Maasland (1959), Hantush (1967) and Marino (1974) have carried out initial pioneering studies on this aspect.

Gill (1984) presented analytical solution of 1D Boussinseq equation linearized using Werner's transformation to describe transient water table profiles in an unconfined horizontal aquifer as a result of seepage occurring from single or more than one canal. Later, Mustafa (1987) incorporated constant replenishment from land surface and developed analytical solution of a linearized Boussinesq equation using Laplace transformation to describe water table variation in a finite aquifer bounded by two recharging canals. Rai & Singh (1992) considered a variable rate of recharge and employing Laplace transformation developed analytical solution of the linearized Boussinesq equation. Ram et al. (1994) also obtained analytical solution of linearized solution of Boussinesq equation by devising unique transformation and showed that water table profiles between two canals were fairly closed with those obtained by Mustafa (1987).

All the above studies either considered no recharge or constant recharge from land surface. Bear (1979) advocated consideration of recharge rate, which is similar to infiltration rate (i.e. decreasing with time in exponential form due to sediment clogging of soil pores beneath the recharge basin). Later, Abdulrazzak & Morel-Seytoux (1983), Zomorodi (1991), Rai & Singh (1996), Manglik et al. (1997), Rai & Manglik (1999), Upadhyaya (1999), Upadhyaya & Chauhan (2001a) considered time-varying recharge rate.

To describe water table fluctuations in a sloping aquifer, Singh et al. (1991) obtained analytical solution of the linearized Boussinesq equation incorporating transient recharge function employing an eigenvalue–eigenfunction expansion method. Rai & Singh (1995) developed analytical solution of 1D Boussinesq equation to obtain water table fluctuation in a finite aquifer due to transient recharge from a strip basin. Ramana et al. (1995) studied water table fluctuations as a result of transient recharge in a sloping two-dimensional (2D) aquifer system. Upadhyaya & Chauhan (2002) obtained analytical solution of the linearized Boussinesq equation and a fully implicit finite difference numerical solution of the nonlinear Boussinesq equation to describe water table rise in unconfined sloping aquifer as a result of seepage from two canals located above the sloping impermeable barrier and constant recharge from land surface. They reported that numerical solution underestimates water table elevations compared to analytical solution, assuming it as the reference solution. Rai et al. (2006) used Fourier Cosine Transform and presented an analytical solution of a 2D-linearized Boussinesq equation to predict water table variations in a horizontal aquifer induced by time-varying recharge/withdrawal from any number of recharge basins, pumping wells and leakage sites. Singh & Jaiswal (2006) presented a numerical solution of 2D free flow of water subjected to time-varying recharge to aquifer underlain by a slanting impervious base and studied the impact of time-varying recharge and depth-dependent ET on water table profile. Song et al. (2007) used perturbation solution of the nonlinear Boussinesq equation for 1D tidal groundwater flow in a coastal unconfined aquifer. Teloglou et al. (2008) presented the analytical solution using Hankel Transform of a linearized Boussinesq equation to describe the water table fluctuation in an unconfined aquifer overlying a semi-impervious layer in response to transient recharge. Bansal & Das (2010) studied the analytical solution of the linearized Boussinesq equation using Laplace and its inverse to characterize transient groundwater flow in a downward sloping unconfined aquifer of semi-infinite extent. Bansal & Das (2011) studied the response of an unconfined sloping aquifer to constant recharge and seepage from the stream of varying water level. Sontakke & Rokade (2014) predicted the water table fluctuations in an unconfined aquifer due to time-varying recharge from the rectangular basin for one canal using the finite difference method. Moshirpanahi et al. (2016) employed the differential quadrature method (DQM) in the discretization of governing equation and compared water table rise between two canals overlying sloping impermeable barrier with those obtained from the discretization of governing equation employing explicit, implicit and Crank–Nicolson numerical scheme-based finite difference methods. The DQM of discretization was reported to be efficient, and results were exactly same as finite difference method. Saeedpanah & Azar (2017) derived analytical solutions for unsteady flow in a leaky aquifer between two parallel streams employing the Laplace transform method and observed that obtained results agreed very well with the results of MODFLOW. Kankarej (2021) developed analytical solution of Boussinesq equation linearized employing Werner's transformation in order to describe water table rise in a sloping unconfined aquifer as a result of seepage from two canals and constant recharge from land surface. Special case for the horizontal aquifer yielded quite close to water table rise compared to the values obtained from solution of Mustafa (1987). Doulgeris & Zissis (2021) studied the numerical solutions for 1D Boussinesq equation and 2D Richards equation by embedding infinite elements in the finite element analysis to discretize only the key subsurface flow region close to the stream.

In India and other countries, canals generally run on the ridge lines. Many a times situation occurs, where two canals separated apart run parallel for a long distance and seepage from these canals as well as recharge from land surface replenish the aquifer. The spatial and temporal rise of water table between two canals can be very well described by obtaining solution of 1D Boussinesq equation.

Analytical solutions of the linearized Boussinesq equation are always approximate and either over- or underpredict the water table profile due to the linearization of governing equation. Numerical solutions of the nonlinear Boussinesq equation are always better than approximate analytical solutions. So, various solutions like hybrid finite analytic solution (HFAS) presented by Chen (1988) and fully implicit finite difference solution (FIFDS) of 1D nonlinear Boussinesq equation and analytical solution of Boussinesq equation linearized by Baumann's transformation (analytical solution I, AS1) as well as analytical solution of Boussinesq equation linearized by Werner's transformation (analytical solution II, AS2) were obtained to describe transient water table variation in a horizontal unconfined aquifer lying between two canals located at finite distance having different elevations and subjected to various patterns of recharge, i.e. zero recharge, constant recharge, as well as time-varying recharge. For comparison of water table profiles, decision about bench mark solution is important.

Upadhyaya & Chauhan (2001b) reported that HFAS, where the nonlinear Boussinesq equation is locally linearized and solved analytically after approximating unsteady term by a simple finite difference formula to approximately preserve overall nonlinear effect by the assembly of locally analytic solutions, predicted fall of midpoint water tables between two drains in a horizontal/sloping unconfined aquifer quite close to the existing experimental results. So, HFAS was considered as benchmark solution in this study and water table profiles obtained from all other solutions were compared with it.

Figure 1 shows the definition sketch of the problem identified for the study. An unconfined aquifer lying over the horizontal impermeable barrier is receiving seepage from two canals located at different elevations h₁ and h₂ above the impermeable barrier as well as time-varying recharge from land surface.

The recharge rate is considered to be decreasing exponentially from an initial value of R₁ + R₀ to a lower value R₀ and becoming constant thereafter. Due to this, there is a water table rise in an unconfined aquifer. It has been assumed that (i) aquifer is homogeneous, isotropic and incompressible with time invariant hydraulic properties, (ii) the rate of recharge is small compared to hydraulic conductivity and vertically percolated water flows almost horizontally after meeting water table and (iii) flow is characterized by 1D Boussinesq equation derived using Dupuit's assumptions and Darcy's Law.

In this section, numerical solutions of the nonlinear Boussinesq equation and analytical solutions of the linearized Boussinesq equation have been presented.

The values of all other basis functions are zero over the elements (X_i−1, X_i) and (X_i, X_i+1). In Equation (22), the multiplier Z_i (T) associated with N_i(X) at node i is the value of H at i. Because there are only two nonzero basis functions over an element (X_i, X_i+1), the summation is performed only over two consecutive indices i and i + 1 in order to approximate the solution H^A (X,T) over the element.

The solution of this system of algebraic equations provided the values of V(T) at different nodes. This V(T) value at a particular node was added to the value of Z(T) at that node to get the value of Z(T + ΔT) at that particular node for the next time step.

This system of algebraic equations formed at a given time step is a tridiagonal matrix for which solution can be obtained and can be computed. To get the values at n + 1 time step, i.e. the values of are added into ⁠, respectively.

Here a = KD/f

To compare water table profiles obtained from analytical and numerical solutions with water table profile obtained from HFAS, L2 norm giving average difference and Tchebycheff norm giving maximum difference as mentioned by Upadhyaya & Chauhan (2002) were employed and accordingly results were interpreted.

The analytical and numerical solutions describing water table fluctuations between two canals in a horizontal aquifer subjected to various patterns of recharge, i.e. zero recharge, constant recharge, and time-varying recharge, were studied and analyzed. A numerical example considered to compare the results obtained from various solutions is given below.

The flow of water in an unconfined aquifer with hydraulic conductivity K = 450 m/day and specific yield f = 0.30, bounded by two canals spaced 1,000 m apart, located at the elevations of h₁ = 12 m and h₂ = 10 m, respectively, above the impermeable barrier was assumed. The aquifer was assumed to be underlain by a horizontal impermeable barrier initially having water table at the impermeable barrier with elevation, h₀ = 0 m. The aquifer was subjected to time-varying recharge, R(t) = R₀ + R₁e^−rt, where R₀ = 0.003 m/day, R₁ = 0.012 m/day and r = 0.5 per day. Consistency, convergence and stability are important issues in numerical solutions. Many values of the dimensionless time increment, ΔT and dimensionless space increment, Δx were considered before arriving at final values of ΔT and Δx as 0.00001 and 0.01, respectively. Chen (1988) has presented the effect of range of λ = (ΔT/(ΔX)²) varying from 0.01 to 10 on different finite difference (FD) scheme coefficients and finite analytic (FA) coefficients. From these curves also, it seems that the selection of λ = 0.1 is appropriate. The transient water tables were computed for t = 5 days and t = 10 days at every 100 m distance. The results obtained from various solutions have been presented below.

Dimensionless water table elevations in a horizontal aquifer between two canals were computed at t = 5 days and t = 10 days by the hybrid finite analytic method for different patterns of recharge, i.e. zero, constant and time-varying recharge. The variations of water table heights with distance expressed in dimensionless form for t = 5 days and t = 10 days are given in Figure 2.

It may be observed from Figure 2 that in a horizontal aquifer, the water table profiles, even with higher seepage from upper canal, tend to be symmetrical around the midpoint. Due to continuing seepage from canals and continuing recharge from soil surface, water table in a horizontal aquifer rises with increase in time. For any pattern of recharge, dimensionless water tables at all the space coordinates are higher for t = 10 days than for t = 5 days. The effect of various patterns of recharge for the horizontal aquifer as observed from Figure 2 indicates that the water table elevations at t = 5 days and at t = 10 days are the highest for constant recharge of 0.015 m/day followed in a decreasing order for time-varying recharge (with constant recharge component, R₀ = 0.003 m/day, and exponentially decreasing recharge component, R₁ = 0.012 m/day) and for no recharge, respectively.

Dimensionless water table elevations in a horizontal aquifer between two canals were computed at t = 5 days and t = 10 days by the finite element method for different patterns of recharge, i.e. zero, constant and time-varying recharge. The variations of water table heights with distance expressed in a dimensionless form for t = 5 days and t = 10 days are given in Figure 3.

It may be observed from Figure 3 that the water table profiles, even with higher seepage from the upper canal, tend to be symmetrical around the midpoint and rise with increase with time. For any pattern of recharge, dimensionless water tables at all the space coordinates are higher for t = 10 days than for t = 5 days. The effect of various patterns of recharges indicates that the water table elevations at t = 5 days and at t = 10 days are the highest for constant recharge of 0.015 m/day followed in a decreasing order for time-varying recharge and for no recharge, respectively. Dimensionless water table heights at t = 5 days and t = 10 days as obtained from FES are throughout marginally higher (with very little difference) at all distances between two canals than those obtained from HFAS.

Dimensionless water table elevations in a horizontal aquifer intercepted by two canals were computed at t = 5 days and t = 10 days using the fully implicit finite difference scheme for various patterns of recharge, i.e. zero, constant and time-varying recharge. The variations of water table elevations with distance expressed in a dimensionless form for t = 5 days and t = 10 days are presented in Figure 4.

It may be observed from Figure 4 that in a horizontal aquifer, the water table profiles even with higher seepage from the upper canal tend to be symmetrical around the midpoint and are almost similar to Figures 2 and 3. Due to continuing seepage from canals and continuing recharge from soil surface, water table in a horizontal aquifer rises with increase in time. The effect of various patterns of recharges is also similar as observed in the case of water table profiles obtained in Figures 2 and 3. Water table profiles computed by FIFDS are marginally higher (with a very little difference) at all the distances between two canals than those obtained from HFAS.

Dimensionless water table elevations in a horizontal aquifer between two canals were computed at t = 5 days and t = 10 days using ASI based on Baumann's method of linearization for various patterns of recharge, i.e. zero, constant and time-varying recharge. The variations of water table heights with distance expressed in a dimensionless form for t = 5 days and t = 10 days are presented in Figure 5.

Water table profiles as obtained by ASI as shown in Figure 5 also follow the same trend as shown in Figures 2–4, but there is difference in values of water table heights. On t = 5 days, water table heights in the midregion are higher and lower in other regions than those obtained from HFAS. But at t = 10 days, water table heights obtained from HFAS are always higher at all distances than those obtained from ASI, and this difference is attributed to the linearization of Boussinesq equation.

Dimensionless water table elevations in a horizontal aquifer between two canals were computed at t = 5 days and t = 10 days using ASII based on Werner's method of linearization for various patterns of recharge, i.e. zero, constant and time-varying recharge. The variations of water table heights with distance expressed in a dimensionless form for t = 5 days and t = 10 days are given in Figure 6.

It may be observed from Figure 6 that in a horizontal aquifer, the water table profiles as obtained by ASII are following the similar trend as shown in Figures 2–5.

Dimensionless water table elevations between two canals in the horizontal aquifer and receiving zero or constant recharge computed for 5 and 10 days by HFAS and ASI based on Baumann's method of linearization were compared graphically and are presented in Figure 7.

It may be seen from Figure 7 that at t = 5 days, water table elevations in a horizontal aquifer receiving zero or constant recharge as predicted by ASI are higher in the midregion and lower in other regions than those predicted by HFAS. At t = 10 days, water table elevations predicted by HFAS are always higher at all distances than those obtained from ASI.

Dimensionless water table elevations between two canals in a horizontal aquifer receiving zero or constant recharge computed for t = 5 days and t = 10 days by HFAS and ASII based on Werner's method of linearization were compared graphically as shown in Figure 8.

It may be seen from Figure 8 that at t = 5 days and at t = 10 days, water table elevations in a horizontal aquifer receiving zero or constant recharge as predicted by ASII are consistently higher at all distances (with more difference in the midregion) than those predicted by HFAS.

Dimensionless water table elevations between two canals in a horizontal aquifer receiving zero and constant recharges computed for t = 5 days and t = 10 days by HFAS and FIFDS were compared graphically as shown in Figure 9.

It may be seen from Figure 9 that at t = 5 days and at t = 10 days, water table elevations in a horizontal aquifer receiving zero or constant recharge as predicted by FIFDS are marginally higher (with a very little difference) at all the distances between two canals than those predicted by HFAS.

Dimensionless water table elevations between two canals in a horizontal aquifer receiving zero and constant recharge computed for t = 5 days and t = 10 days by HFAS and FES were compared graphically as presented in Figure 10.

It may be seen from Figure 10 that at t = 5 days and t = 10 days, dimensionless water table heights in a horizontal aquifer receiving zero or constant recharge as predicted by FES are throughout marginally higher (with a very little difference) at all the distances between two canals than those predicted by HFAS.

Prenter (1975) described about L2 and Tchebycheff norms and Upadhyaya & Chauhan (1998) employed these L2 and Tchebycheff norms to compute average and maximum differences, respectively, between the HFAS (benchmark) and other fully implicit finite difference numerical solution as well as linearized ASI and ASII. Values of L2 and Tchebycheff norms indicating average and maximum differences between two solutions were computed and are given in Table 1.

It may be observed from Table 1 that for both t = 5 days and t = 10 days, L2 and Tchebycheff norms were minimum for FES followed by FIFDS, ASI based on Baumann's method of linearization and maximum for ASII based on Werner's method of linearization. It indicates that the theoretical performance of various solutions with respect to HFAS in a decreasing order may be ranked as FES, FIFDS, ASI and ASII, respectively.

Water table rise in an unconfined horizontal aquifer lying between two parallel canals due to time-varying recharge from land surface and seepage from canals was studied by obtaining analytical and numerical solutions of 1D Boussinesq equation incorporating time-varying recharge term. Analytical solutions of Boussinesq equation linearized by Baumann's transformation (AS I) as well as Werner's transformation (AS II). FES and FIFDS of nonlinear Boussinesq equation were developed and compared with HFAS of nonlinear Boussinesq equation (considering it as benchmark solution). Comparison of water table profiles as well as L2 and Tchebycheff norm values showed that water table profiles at t = 1 day and t = 5 days as obtained by FES were quite close to the water table profiles obtained from HFAS, FIFDS followed by FES. Therefore, FES was ranked at the first place and FIFDS at the second place. As far as analytical solutions are concerned, at t = 5 days, both AS I and AS II predicted higher water table values, but at t = 10 days, AS I predicted lower and AS II predicted higher water table values at all distances due to the linearization effect. According to L2 and Tchebycheff norm values, AS I was ranked at the third place and AS II at the fourth place.

The authors thankfully acknowledge help and support provided directly or indirectly by individuals and respective institutions in carrying out this study.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.