Modeling of water table profile variations owing to stream–aquifer interaction

In arid and semi-arid developing countries, such as India, Egypt, and Pakistan, canal irrigation is practised extensively. Providing water for irrigation, the implementation of a canal system also acts as a basis for rigorous seepage below the ground level (bgl) in highly previous tracts (Maheswaran et al. 2016). Seepage occurrence originates from the extensive lack of canal alignment, due to a lack of funding for its completeness (Elkamhawy et al. 2021). Due to continuing seepage from canals and percolation losses from irrigated fields, gradually over the years the water table (WT) has risen in some areas (Morway et al. 2013). On the other hand, excessive pumping and unplanned withdrawal of groundwater may lead to a lowering of the WT, which may cause problems such as a reduction in future water supply, depletion of aquifers, seawater intrusion in the coastal areas, land subsidence etc. (Yin et al. 2020). It is important to describe and understand the phenomenon of the rising or lower of the WT and its temporal and spatial distribution in a semi-infinite aquifer as a result of stream–aquifer interaction (SAI) and/or recharge from the land surface in order to develop policies for efficient and judicious use of water (Jeet et al. 2019). SAI is a proper understanding of surface water and groundwater interaction which is crucial for effective and sustainable water resource management (Prajapati et al. 2021). It may happen in three distinct ways: the stream is either (1) losing stream water infiltrates into the aquifer, (2) gaining – groundwater flows into the stream, or (3) disconnected – losing stream that is disconnected from the aquifer by an unsaturated zone. Mathematical modeling is an important technique to describe the physical situation in numerical terms and predict the rising/lower in WT in the semi-confined aquifer. The basic governing equation for describing WT behavior has been the continuity equation based on Dupuit–Forchheimer assumptions obtained by Boussinesq.

Analytical and numerical methods may be employed to solve the Boussinesq equation (BE) describing WT distribution in a horizontal semi-confined aquifer due to SAI. Edelman (1947) deliberates a problem in which the ditch water level rises or falls, causing water to drain out of, or into the adjoining aquifer until the WT and the ditch levels are again in equilibrium conditions. The effect of water level changes in the ditch was described by the linearized BE and expressed solution using the error function. Later on, Polubarinova-Kochina (1948); Polubarinova-Kochina (1949); Marino (1973); Lockington (1997); Upadhyaya & Chauhan (2001); Huang et al. (2010); Spanoudaki et al. (2010); Teloglou & Bansal (2012); Bansal et al. (2016) and Zhou et al. (2018) studied problem of SAI under assigned IBC considering horizontal/sloping aquifer with constant or no recharge from land surface by employing various analytical and numerical techniques. Saxena et al. (2021) derived and implemented a closed-form analytical equation for the simulation of hydraulic head (H) distribution in an SAI under the influence of stream-stage variations and percolation. The aquifer was assumed to be intrinsic by a sloping impermeable bed and variations in both the kinematic stage and recharge rate. They obtained that due to continuous recharge significant rise of the WT in the mid-part of the aquifer. Upadhyaya & Kankarej (2022) obtained AS and NS of one-dimensional (1D) BE to describe WT rise in a horizontal aquifer lying at a finite distance between canal systems having different heads. Considering Hybrid Finite Analytic Solution (HFAS) as the benchmark solution, Finite Element Solution (FES) and Fully Implicit Finite Difference Solution (FIFDS) were ranked at first and second place, respectively, and analytical solutions (ASs) of the BE, (supposed to be approximate due to linearization) were ranked lower than the NSs of nonlinear BE.

Spatial and temporal variation of WT in recharging and discharging aquifers has been studied by applying ASs and NSs of the 1D BE. ASs may not be as accurate as NSs due to the linearization of the BE. In order to study the effect of linearization done in ASs a comparison in spatial and temporal variation as obtained by AS and NS is essential. The objective of the study is to develop HFAS, FIFDS, FES, AS I, and AS II in order to obtain spatial and temporal WT distribution in semi-confined aquifer due to SAI and compare WT profiles obtained by employing AS and NS with WT profile obtained employing HFAS (where nonlinear BE is linearized locally and solved analytically, unsteady term is approximated by simple finite difference equation and overall nonlinear effect is preserved by the assembly of analytic solutions), considering it as the benchmark solution.

In Werner's transformation-based AS II, in nonlinear Equation (1), the term (∂h/∂x)² is absorbed by substituting z = h² and ‘h’ attached with the derivative term is replaced by D, the average depth of flow.

In this sub-section, NSs of nonlinear BE and AS of linearized BE have been presented.

Chen (1988) presented the concept of HFAS, where the nonlinear BE is linearized and solved analytically after approximating unsteady terms by a simple finite difference equation to approximately preserve the overall nonlinear effect by the assembly of analytic solutions.

Assuming the terms 1/D (∂v/∂x)² and 1/D (∂v/∂t) equal to constants C₁ and E₁, respectively. Again by applying inverse transformation Equation (3), the value of h at a given space and time is obtained.

h^A (x,t) is an approximation for h(x, t). Hence, its substitution in Equation (1) leaves a residual L(h^A) which is used to determine the coefficients z_i (t). In Galerkin's finite element method, the coefficients z_i (t) are determined by forcing the residual L(h^A) to be orthogonal to the basic functions N_i (x), i = 1, 2, 3, …, N.

A numerical solution to nonlinear BE incorporating constant recharge and constant evaporation with IBC, in a semi-infinite flow region was obtained by applying FIFDS. While analyzing the flow problem in a semi-infinite region with the help of any numerical technique, it is assumed that the infinite distance means an extensively large distance where the WT profile becomes asymptotic. So in place of a theoretically infinite value, a large value of distance may be considered in the analysis. Since the location of the point of infinity is not fixed and its position also changes with time, the non-dimensionalization of the governing partial differential equation was not feasible and the solution was obtained for its original dimensional form.

The system of algebraic equations formed at a given time step is a tri-diagonal matrix for which solution is obtained by standard algorithm available in the study by Jain et al. (1994) of numerical analysis and computed.

The above solution is in the form of z(x, t). The solution in the form of h (x, t) can be obtained.

The analytical and NSs explaining the variation of WT in the horizontal aquifer due to sudden rising or lowering of water level in the adjoining stream or canal were studied by considering a numerical example. The effect of constant infiltration and ET from the land surface on WT in both conditions (recharging and discharging) of STI in a semi-infinite flow region was also studied. A numerical example, which was considered for comparison of results obtained from different solutions, is given below.

The water flow in shallow aquifer having hydraulic conductivity K = 20 m/day and specific yield f = 0.27 was assumed. The aquifer was supposed to be underlain by a horizontal impervious barrier and to initially have a uniform WT elevation h₀ = 2 m. The water level was supposed to be instantaneously raised in the adjoining trench/canal to the elevation h₁ = 3 m to provide for an interacting recharging aquifer. Likewise, when WT in the aquifer was at an elevation of h₀ = 3 m, and the level in the canal/trench was at an elevation of h₁ = 2 m, it provided an interacting discharging aquifer. A constant surface infiltration (SI) of 0.05 m/day and ET of 0.008 m/day were assumed in order to study the effect of these factors on the WT in the aquifer. In NSs, values of time increment (Δt) and space increment (Δx) were assumed as 0.0025 days and 2 m, respectively. The resulting WT profiles in the aquifer for both after 1 day and after 5 days cases were determined using various solutions.

Results of HFAS, FES and FIFDS for zero SI and ET, in that instance recharging aquifer the WT elevations are higher at 5 days than at 1 day whereas in case of discharging aquifer, the WT elevations are higher at 1 day than at 5 days. For both conditions of the aquifer, the WT head near the interface is steeper at 1 day than at 5 days. Due to the combined effect of constant SI @ 0.05 m/day and ET @ 0.008 m/day from the land surface, i.e. due to net downward flux, the WT in both conditions of the aquifer for 1 day and 5 days are observed to be higher than the case when zero SI and ET have been considered. It may also be observed from Figure 2(a)–2(c) that only in the case of recharging aquifer does the groundwater mound form near the stream at 5 days, which decays with an increase in distance away from the stream and becomes parallel to the impermeable barrier.

Transient WT in horizontal recharging and discharging aquifer due to sudden rise or fall of water level in the adjacent stream or canal for 1 day and 5 days predicted by HFAS were compared with the WT estimated by Edelman (1947), Polubarinova-Kochina (1948) and Lockington (1997) solutions.

It may be observed that WT in both conditions of the aquifer estimated by Polubarinova-Kochina (1948) solution are very close to the WT estimated by HFAS. The absolute maximum difference in values of WT estimated by Lockington (1997) and HFAS is more than the difference between Polubarinova-Kochina (1948) and HFAS but less than the absolute difference between WT estimated by Edelman (1947) solution and HFAS.

Transient WT in horizontal recharging and discharging aquifers as a result of sudden rising or lowering of water level in the adjacent stream or canal with and without infiltration and ET from land surface were computed for 1 day and 5 days by HFAS and AS I. It may be observed that for recharging and discharging horizontal aquifers, without considering SI and ET from land, up to a certain distance, the WT heights estimated by AS I are marginally different (higher or lower) than those estimated by HFAS (Tables 1 and 2). Beyond this distance both the solutions predict similar constant values of WT head for 1 day and 5 days. If the effect of SI and ET with net downward flux is taken into account almost a similar trend of WT heads in the case of zero SI and zero ET is observed (Tables 3 and 4). Due to net downward flux, the WT heights in recharging and discharging horizontal aquifers estimated by both the solutions are relatively greater at all the values of distance and time than those obtained without considering constant SI and ET.

Transient WT in both conditions of aquifer having zero and 10% slope as a result of sudden rising or lowering of water level in the adjacent stream or canal with and without infiltration and ET from the land surface were computed for 1 day and 5 days. It may be observed that for recharging and discharging horizontal aquifers, without considering infiltration and ET from land, up to a certain distance, the WT heads estimated by AS II are marginally higher than those estimated by HFAS (Tables 1 and 2). Beyond this distance, both the solutions predict almost similar constant values of WT heads for both 1 day and 5 days. If the effect of constant SI and constant ET with net downward flux is taken into account, in a horizontal recharging aquifer, for 1 day and 5 days WT elevations estimated by AS II become higher than those obtained from HFAS after 20 and 70 m distance from the stream, respectively (Table 3). Similarly, in a horizontal discharging aquifer, for 1 day and 5 days, WT heads estimated by AS II are lower than those obtained from HFAS after 40 and 10 m distance from the stream (Table 4). Due to net downward flux, the WT heads in recharging and discharging horizontal aquifers predicted by both the solutions are relatively more at all the distances than those obtained without considering constant SI and ET.

Transient WT in horizontal recharging and discharging aquifers as a result of sudden rising or lowering of water level in the adjacent stream or canal with and without SI and ET from land surface were computed for 1 day and 5 days by HFAS and FIFDS. It may be observed that for horizontal recharging and discharging aquifers, without considering SI and ET from land, for 1 day and 5 days, the WT heads estimated by FIFDS are almost identical to those estimated by HFAS at all distances from the stream or canal (Tables 1 and 2).

If the effect of SI and ET with net downward flux is taken into account, up to a certain distance, in a horizontal recharging aquifer for 1 day and 5 days WT elevations estimated by FIFDS are marginally lower (Table 3). Similarly, in a horizontal discharging aquifer, WT elevations are marginally higher than those obtained from HFAS (Table 4). Beyond this distance WT heads in the recharging and discharging horizontal aquifers estimated by both the solutions become identical. Due to net downward flux, the WT heads in horizontal recharging and discharging aquifers estimated by both the solutions are relatively more at all distances and times than those obtained without considering constant infiltration and ET.

Transient WT in horizontal recharging and discharging aquifers as a result of sudden rising or lowering of water level in the adjacent stream or canal with and without SI and ET from land surface were computed for 1 day and 5 days from HFAS and FES. It may be observed that for recharging and discharging horizontal aquifers, without considering SI and ET, the WT heads estimated by FES are identical to those predicted by HFAS for 1 day and 5 days at all distances from the stream or canal (Tables 1 and 2).

If the effect of constant SI and constant ET with net downward flux is taken into account, in a horizontal recharging aquifer for 1 day and 5 days WT elevations predicted by FES are marginally lower, up to a certain distance (Table 3). Similarly, in a horizontal discharging aquifer WT elevations are marginally higher than those obtained from HFAS (Table 4). Beyond this distance WT heads in the recharging and discharging horizontal aquifers predicted by both the solutions become identical. Due to net downward flux, the WT heads in recharging and discharging horizontal aquifers predicted by both the solutions are relatively more at all distances and times than those obtained without considering constant SI and ET.

Prenter (1975) described L2 and Tchebycheff norms, which indicates differences between the two solutions were measured to theoretically rank the performance of various ASs and NSs with respect to HFAS. The computed values of L2 and Tchebycheff norms for horizontal recharging and discharging aquifers consider zero or constant SI and ET for 1–5 days.

It may be observed that for a horizontal recharging aquifer without considering SI and ET from land, L2 and Tchebycheff norms were minimum for FES and FIFDS followed by Polubarinova-Kochina (1948) solution. Among Lockington (1997) solutions, AS I and AS II these two norms could not consistently establish the order of performance of these solutions because the smaller value of the L2 norm was observed with AS I and the smaller value of Tchebycheff norm was observed for Lockington (1997) solution. Edelman's (1947) solution was observed to have the lowest rank of performance because of the highest values of L2 and Tchebycheff norms.

If constant SI and ET from land are taken into account for studying the WT elevations in a horizontal recharging aquifer, it may be also observed that the values of L2 and Tchebycheff norms were minimum for FES with successively increasing values for FIFDS, AS I and AS II, respectively. Thus, the performance of FES with respect to HFAS was ranked at first place whereas the performance of AS II was ranked in fourth place. It may be observed that for a horizontal discharging aquifer without considering SI and ET from land, L2 and Tchebycheff norms were minimum for FES and FIFDS followed by Polubarinova-Kochina (1948). Among Lockington (1997) solutions, AS I and AS II these two norms could not consistently establish the order of performance of these solutions because the smaller value of the L2 norm was observed with AS I and the smaller value of Tchebycheff norm was observed for Lockington (1997) solution. Edelman's (1947) solution was observed to have the lowest rank of performance because of the highest values of L2 and Tchebycheff norms. If constant SI and ET are taken into account for studying the WT elevations in a horizontal discharging aquifer, it may be observed that the values of L2 and Tchebycheff norms are minimum for FES and FIFDS with successively increasing values for AS I and AS II. Thus, the performance of FES with respect to HFAS was ranked at the first place whereas the performance of AS II was at the fourth place.

1D BE after incorporating constant SI and ET terms with appropriate IBC was found appropriate to describe SAI in a semi-infinite flow region. ASs were developed by AS I, AS II, FIFDS, and FES. These solutions were compared with HFAS considering it as a benchmark solution. For a selected example, the WT profile in the case of recharging and discharging aquifers on day 1 and day 5 was obtained by applying AS and NS. Comparison of spatial and temporal variation of WT estimated by various ASs and NSs as a result of SAI in both the cases of recharging and discharging aquifers subject to zero or constant SI and ET as well as L2 and Tchebycheff norm values showed that FES values were quite close to HFAS. So, FES ranked at first place and FIFDS at second place. For the case of zero SI and ET from land, Polubarinova-Kochina (1948) solution ranked at third place and Edelman (1947) ranked at last place. As evident from L2 and Tchebycheff norms, among ASs, the performance of AS I was better than AS II in both recharging and discharging aquifers subjected to constant SI and ET, so the performance of AS I and AS II were ranked at third and fourth place, respectively.

From the above comparative study, it may be concluded that HFAS (in which overall nonlinear effect is approximately preserved due to assembly of local ASs of nonlinear BE after approximating unsteady term by a simple finite difference equation) is mathematically simple, faster, robust, and more accurate than other techniques in studying the spatial and temporal WT variation in recharging and discharging aquifers in semi-infinite flow region.

The authors thankfully acknowledge the 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.