Numerical groundwater modelling for studying surface water-groundwater interaction and impact of reduced draft on groundwater resources in central Ganga basin

1.   Introduction

India is experiencing recurring water challenges of both extreme deficit (droughts) and surplus (floods) as a result of climate change and weak water resources management system ^{[1, 2, 3]}. Increasing sectoral demands for water resources due to increasing population, industrial manufacturing, and agriculture production has put tremendous pressure in managing the freshwater resources. Further, unequal distribution of water resources across geographies and seasons, climate change, and unplanned urban growth are putting additional pressure in resource management ^{[4, 5]}.

According to one of the surveys conducted by the Tata Institute of Social Sciences (TISS), most metropolitan centers are already water-stressed and primarily focusing on groundwater resources for meeting their water demands. This results in decreasing groundwater tables at an alarming pace in the major cities. In addition to this, about 39 million ha of irrigated land in India is accomplished through supply from groundwater ^[6].

According to the Central Ground Water Board (CGWB) report (2020), India's total annual replenishable groundwater resource is roughly 436 billion cubic meters (BCM), with 398 BCM net annual groundwater availability, of which India withdraws 245 BCM (62%) yearly. Moreover, around 39% of the wells show a decrease in groundwater level. Out of the country's 6965 assessment units, 1114 units have been labelled as "overexploited" due to groundwater depletion and long-term reduction in groundwater levels. Although statutory bodies such as the Central Ground Water Board & Central Water Commission are set-up to make policies for sustainable development and management of groundwater and surface water resources, but understanding the best practices from other nations and India's own community-based intervention models is essential in improving the governance structures and understanding critical indicators that can aid in data-driven decision-making.

The middle Gangetic plain (MGP) region is one of the most significant eco-regions which form the core portion of the Indo-Gangetic plains (IGP) ^[7]. Because of its dominating rice-wheat farming system, the MGP provides a living for millions of people ^[8]. In the MGP, several obstacles hinder the implementation of the critical measures that boost food production. Some of the key concerns existing in the middle Ganga basin are decreased net cultivated area, decreased irrigation water availability, and an increased strain on soil fertility ^{[9, 10]}. The use of conservative agriculture techniques and other farm-scale water-saving measures are the potential remedies known for combatting the unsustainable use of groundwater. However, these issues are normally considered in isolation ^[11]. Moreover, a prior scientific approach is imperative to understand the regional groundwater and surface water interaction vis-a-vis its' impact on the groundwater levels. Numerical groundwater modelling, now-a-days, is increasingly used to understand the groundwater availability and use, surface water-groundwater interaction and its impact on the water resources management and planning ^{[12, 13, 14, 15]}.

In the present study, numerical groundwater modelling has been performed for studying the impact of surface water-groundwater interaction and reduced draft on groundwater resources in the middle Ganga river plain region of India with the following major two objectives: i) trend analysis of groundwater level data to understand recharge and discharge processes in the study area; ii) to study surface water and groundwater interaction based on the available and monitored data. Findings from this study will enable the water resources manager and policy makers particularly working in the middle/central Ganga basin to efficiently plan and manage the water resources through conjunctive water use.

2.   Study area

The Bhojpur district of Bihar, India was selected as study area which lies between 25º 10' to 25º 40' N latitude and 84º 10' to 84º 50' E longitude, covering geographical area of 2395 Km². The monsoon season starts during the month of June and continues up to end of September. Two major rivers cover the district, viz., Ganga and Son in the north and eastern part of the study area respectively (Figure 1). The study area has a warm and humid climate. The alluvial formation (younger and older alluvium) occurs in the study area (Figure 1), in which north and north-eastern part are covered by younger alluvium whereas central and southern parts are formed by older alluvium plains ^[16]. The areas covered by different crops in Kharif, Rabi and Summer season are indicated in Table 1 and net irrigation area under Kharif and Rabi season are shown in Table 2.

3.   Materials and methods

For achieving the objectives framed under study, the following methodology was adopted for the study: i) a review of the literature was done to understand the impact of agricultural water-saving measures on hydrology; ii) various thematic map such as drainage, geological, depth to water level, LULC map for the study area were prepared; iii) GWL data were analysed to study spatio-temporal variation and trend analysis was done using Mann- Kendall test; and iv) based on the availability and generated data, a coarse groundwater modelling was attempted to study the SW-GW interaction. The detailed methodology for trend analysis and groundwater modelling is described below.

3.1.   Mann-Kendall test for groundwater level trend analysis

Long term trend in groundwater levels is best-known by using non-parametric Mann-Kendall (MK) test ^{[17, 18, 19, 20]}. The MK test is commonly used to detect monotonic trends. Assumptions of a normal distribution are not required in the MK test, which only indicates the significance and direction of the trend, not its quantitative measure ^[21]. Assumptions of the MK test include a stable, independent, and random time series with an equal probability distribution ^[22]. The time-series data used for the MK test should be uncorrelated; otherwise, it might result in an erroneous rejection of the null hypothesis due to serial correlation ^{[21, 22, 23, 24]}. The computational procedure for the MK test is described in ^[25].

The MK test statistic, S, is defined as

Where, ${T}_{1}, {T}_{2}, {T}_{3}, .....{T}_{n}$ represents $n$ data points; ${T}_{j}$ represents data at the time $j$ of data time series; and ${T}_{j}-{T}_{i}=\theta$.

The statistic S is almost normally distributed when n ≥ 0, with the variance (${\sigma ^2}$) as given below:

Where, t_i denotes the number of ties to the extent i. The standardized test statistic (Z_s), is given by:

Equation (4) is useful when the data length is ≥ 10 and the number of tied data is less ^[26]. A positive and negative Z_s values indicate an upward and downward trend, respectively. If |Z_s| is greater than Z_α/2, where α represents the chosen significance level (usually 5%, with Z_0.025 = 1.96), the null hypothesis is rejected, implying that the trend is significant.

MK test holds good in case of non-auto correlated time series data. A modified Mann-Kendall (MMK) test proposed by ^[27] is used if there exist auto-correlated data. The modified variance of S is given by Eqs (5) and (6).

Where, n is the actual number of observations and ${\rho }_{S}\left(i\right)$ is the autocorrelation function of the ranks of the observations. The autocorrelation between ranks of observations ${\rho }_{S}\left(i\right)$ is firstly evaluated after subtracting a suitable non-parametric trend estimator.

3.2.   Groundwater model development

The groundwater (GW) model is used to understand the behaviour of the GW flow system and predict the system's response due to forcing function. It also acts as an assessment tool for evaluating recharge, discharge, and quantifying the sustainable aquifer yield. Typically, a GW flow model comprises of an equation (algebraic or differential) governing the groundwater flow and numerical solution to compute the head distribution in space and time.

The mathematical model of a groundwater system is defined by a governing equation (partial differential equation), and initial and boundary conditions. The governing equation is essentially an expression of the continuity equation i.e., the difference between inflow and outflow rates equals the rate of change of storage in any selected domain of saturated flow ^[28]. The governing equation is derived by using the continuity equation and Darcy's law.

The mass flux (mass/time) in through the back face X-axis of the elementary boundary is given by

Where ${\rho }_{w}\left(\mathrm{x}\right)$ is the water density and ${q}_{\mathrm{x}}\left(\mathrm{x}\right)$ is the specific discharge at coordinate $\mathrm{x}$.

The mass flux out through the front face of the element is

As per mass balance, the mass influx minus mass outflux is the rate of change of mass stored in the element.

Rearranging the Eq (9) and using Darcy's law, the general GW flow equation can be written as

If aquifer is assumed to be homogeneous (K is independent of x, y, and z), then the general equation can be written as

If the aquifer is assumed to be homogeneous (K independent of x, y, and z) and isotropic (K_x = K_y = K_z), the general equation can be written as

The partial differential equation is converted into set of algebraic equation and then it is solved by computer program or codes. The groundwater head can be calculated by solving three dimensional groundwater flow equation using MODFLOW software ^[29]. The finite difference method (FDM) is used to solve the equations.

3.3.   Modelling process

The modelling process involves several steps, which are shown in Figure 2. The first step in the modelling process is to define the purpose. The modeller needs to answer the questions, viz., whether the model is to be used for prediction, system understanding, generic exercises, etc, before starting the modelling exercise.

3.4.   Conceptual and numerical model construction for the study area

The conceptual model is the base for model analysis. A conceptual model simplifies the actual field conditions so that groundwater system may be analyzed more willingly ^[30]. The conceptualization of model includes creating and framing data related to geology, hydrogeology, hydrology, meteorology, aquifer geometry, geologic information, etc. The numerical model is designed based on the developed conceptual model. The design of the numerical model comprises the size of model grids, defining time steps and stress periods, and defining initial and boundary conditions of the model. The grid size and time steps depend on the purpose of the modelling, as the memory and computing time of computers/computer code may have limitations. The present model in this study was conceptualized with a one-layer model (unconfined aquifer system); the unconfined aquifer thickness is taken as 80 m. The four stress periods (i. July to September; ii. October to December; iii. January to March; iv. April to June) were taken to simulate GW flow daily. The calibration was done for one year (July 2015 to June 2016), validation for one year (July 2016 to June 2017), and then prediction was made with the help of a developed model. The input data viz. specific yield, hydraulic conductivity, initial groundwater head, recharge and groundwater draft used in this model has been taken from CGWB report and other literature, and available river data has been gathered from the website of Water Resources Department, Govt. of Bihar.

The geographical area of the study domain is 2395 Km^2, and the grid size was taken 1 Km × 1 Km with 64 rows and 59 columns. The model domain is defined as

The developed model has certain limitations which have been mentioned here. The calibration and validation period was less due to the non-availability of required data. This is a scoping study, and a coarse GW modelling was attempted based on the available data with certain assumptions. We have assumed that the aquifer is homogeneous and isotropic, however, inclusion of heterogeneity may be desired for investigating the local problem. An improved future model should incorporate heterogeneity in specific yield and hydraulic conductivity. The recharge values have been taken from CGWB report and these recharge estimates were not available for the simulation period and therefore recharge values of year 2013 were used for the entire simulation.

4.   Results

4.1.   Hydroclimatology

Bhojpur district lies in a sub-humid region with a warm climate. Figure 3 shows the mean monthly precipitation for the duration (1991–2020). The average monthly rainfall occurs highest in the month of July, accounting for 266 mm. Almost ~87% of average annual precipitation occurs during the Kharif season (July to October), and most of the rainfall (about 85%) occurs in the South-West monsoon ^[31]. The cumulative rainfall from January to December, reaching its maximum limit of 950 mm. The hottest month is May when the average monthly temperature reaches 35℃ whereas January experiences the lowest mean monthly temperature as ~9℃. The average monthly temperature ranges between 25 to 35℃ in the Kharif season while it varies from ~15 to 25℃ in the Rabi season.

4.2.   Hydrogeology

The Bhojpur district is occupied by Quaternary alluvium, forming the potential aquifer. It is covered by alluvial formation; the northern part is enclosed by younger alluvium, whereas the central and southern parts are covered by older alluvium, as discussed in the study area section. The older alluvium of the study area consists of dark colored clay and silt rich with lime nodules locally known as kankars. The unconsolidated younger alluvium occurs along the flood plain of the Ganga and Son rivers, and it is characterized by sandy clay, and loam and contains less calcerous matter. The fence diagram (prepared based on borelog data) shows that the top layer of the geological stratum (within 30 mbgl) is an aquitard (Figure 4), which supports dug wells and shallow hand pumps. It works as an unconfined aquifer. From 30 m to approximately 80–100 m bgl, medium to coarse sand forms the aquifer, and after that, a thick layer (20–30 m) of clay is present. The aquifers which are at depth > 130 m bgl are under either semi-confined or confined conditions ^[31].

The groundwater level data analyses (data collected from CGWB) show groundwater level variation across the study area. The groundwater level was found to be higher in the North and North-East part of the study area. The unconfined aquifer has an average groundwater level ranging from 3 m bgl to 9 m bgl (pre-monsoon season, May 2018); the hydraulic gradient is about 0.0005, and groundwater flow directions is towards North and North-East (Figure 5) towards the river Ganga ^[32]. The groundwater level fluctuation can be noticed between pre and post monsoon season (Figure 6 has been shown as representative for a well hydrograph), indicating that natural recharge is good in the study area ^[33].

4.3.   Land use/land cover

Landsat-8 satellite imagery (year 2021, 10 m resolution, USGS) data was used for the preparation of land use land cover (LULC) map of the study area. The LULC classification (Figure 7) shows the major land use/land cover types as agricultural land (64.07%), barren land (5.65%), built-up area (24.75%), scrub/shrub (2.13%), vegetation (1.09%), and water (2.31%).

4.4.   Trend analysis of groundwater level

Groundwater level data for the present study comprises pre-monsoon, monsoon and post-monsoon monitored data from 13 wells within the study area, as shown in Figure 8. The data length for the monitoring wells ranges from 5 to 18 years as per the availability. Therefore, 3 time series (pre-monsoon, monsoon, and post-monsoon) for each well were prepared for trend analysis. A total of 39-time series data were analyzed using the Mann-Kendall test method, and the results are given in Table 3. The limitation of the trend analysis is in the short duration of data availability. Nevertheless, the present analysis provides an overview of the trends of rising and declining groundwater levels in the study area. It can be seen from Table 3 that only one station, namely, Kulharia-1 shows a statistically significant (@ 5% significance) declining trend of groundwater level during pre-monsoon as well as monsoon period. The trend remains declining during the post-monsoon season also, but it is not significant. The trends are not statistically significant for the rest of the 12 wells. The reason could be the short length of available data. Five more wells (Jagadishpur-1, Udwantnagar, Pirro-2, Milki, Pirro-1) show a declining trend (not significant) of groundwater levels, mainly in the pre-monsoon and monsoon periods with no trend during the post-monsoon period. Two wells (Bihiya, Garhani-1) in the study area show rising trends of groundwater level while two other wells (Garhani-2, Chandawa-1) show increasing trends in pre-monsoon period and declining trends in monsoon period with no trends in the post-monsoon period. There are three wells (Kulharia-2, Chandawa-2, Jagdishpur-2) that show neither rising nor declining trends of ground water levels.

4.5.   Groundwater modelling

4.5.1.   Model inputs

The inputs to the groundwater model include initial (initial heads) and boundary conditions, hydrogeological parameters such as transmissivity/conductivity and storability/specific yield, and hydrological stresses. The initial heads were interpolated from observed groundwater heads (May 2015) for the unconfined aquifer. Hydraulic conductivity (K) and storage coefficient (specific yield, Sy) are the two parameters of an unconfined aquifer that define the movement and storage of groundwater. The value of conductivity, and Sy were taken from the literature ^[31] and assigned to a single-layer unconfined aquifer. The hydraulic conductivity values for the younger and older alluvium were assigned as 40 m/day and 30 m/day, respectively, as shown in Figure 9. A specific yield of 0.10 was taken for the study domain ^[33].

4.5.2.   Boundary conditions

The physical boundaries in North and East of the study area are the Ganga and Son rivers and hence river boundary conditions were used. No flux boundary was assigned in the West (as plotted water table contours are parallel, and therefore flow is not possible from west to east), and flux boundary was assigned in the southwest (Figure 10) part of the study area. The flux for an unconfined aquifer was calculated based on the hydraulic gradient and hydraulic conductivity. The river heads and bed bottom elevations in the upstream and downstream of river Ganga (in the study area) are 52 and 48 meter above mean sea level (m amsl) and 48.75 and 44.75 m amsl, during 1^st stress period; 51 and 48 m amsl and 47.75 and 44.75 m amsl, during 2^nd stress period; 50.5 and 48 m amsl and 47.25 and 44.75 m amsl, during 3^rd stress period; 50 and 48 m amsl and 46.75 and 44.75 m amsl, during 4^th stress period respectively. Similarly, the river heads and bed bottom elevations at the start and end point of river Son (in the study area) are 74 and 72 m amsl and 46.75 and 44.75 m amsl, during 1^st stress period; 73.5 and 72 m amsl, 46.25 and 44.75 m amsl, during 2^nd stress period; 73 and 72 m amsl, and 45.75 and 44.75 m amsl, during 3^rd stress period; 72.5 and 72 m amsl and 45.25 and 44.75 m amsl, during 4^th stress period respectively. The river head data during monsoon season were taken from Water Resources Department, Govt. of Bihar website and it was assumed for other periods. The river widths were taken from Google earth.

The hydrological stresses considered in this study area are recharge and groundwater draft. The annual recharge value of Bhojpur district (block-wise), has been taken from CGWB (2013) report and has been distributed in the stress period (75, 15, 5 and 5% during 1^st, 2^nd, 3^rd and 4^th stress period). The block-wise recharge values were assigned to the model for different stress periods as shown in Table 4.

The groundwater draft was also taken from the CGWB report (year 2013) for the Bhojpur district. The annual GW draft was distributed block-wise for four stress periods, as shown in Figure 11 and Table 5. A total of 288 wells were estimated to be operated in the study area. The withdrawal pumping rate for stress periods 1^st, 2^nd, 3^rd & 4^th are -1800, -3022, -3445, and -4607 m³/day/pump, respectively. The well package of MODFLOW was used for incorporating the GW draft in the model.

4.5.3.   Calibration of the model

The model was calibrated by tweaking the input parameters (K, Sy, recharge, draft, etc.) to reproduce the field measured groundwater heads. The alteration of aquifer parameters or stresses was done by trial-and-error methods. The calibration for the unconfined aquifer was made for one year (1^st to 365 days), and simulated head values were compared with observed values of monitoring wells.

As reported in the literature, hydraulic conductivities were used as initial values for the simulation. By trial and error calibration, the conductivity values for younger and older alluvium were estimated to be 30 m/day and 20 m/day. The specific yield was adjusted as 0.08. The pumping rates were adjusted within 15% for calibrating the model. The scatter plots between observed and calculated heads for the aquifers are presented in Figure 12. The accuracy between simulated and observed groundwater levels is measured by model performance indicator such as residual mean, standard error of estimate, root mean square (RMS), etc. The error statistics with correlation coefficient (R) for the present model is given in Table 6. The performance in terms of RMS and standard error is very good and R value is above 0.9 which is excellent. These results demonstrate that model is very well capable of making predictions.

The simulated GW heads for the unconfined aquifer during different stress periods are shown in Figure 13. The GW heads are higher on the southern side and lower in the northern part, and overall GW head fluctuation is 2 m from the 1^st stress period to the 4^th stress period. The water budget for different stress periods of the groundwater system is shown in Table 5. It is observed that groundwater contributes more in monsoon season to rivers than in non-monsoon season. It is also evident that recharge is more in the monsoon season, and the pumping rate is maximum in the pre-monsoon season (stress period 4).

4.5.4.   Validation of the model

The model has been verified with another set of field data (366 days to 730 days) to check the model capability for predicting groundwater head. The calibrated model should be able to predict head and should match with observed values otherwise re-calibration may be required. For the present case, the validation period was from July 2016 to June 2017. The scatter plots between observed and calculated heads confirm a good match for different stress periods, as shown in Figure 14.

4.5.5.   Prediction of groundwater heads

The calibrated and validated model was used to forecast the future changes that can occur under different stresses. The impact of reduction of pumping (GW withdrawal) on GW heads was investigated. The 10 % reduction in groundwater withdrawal was made during all stress periods of the year. It was observed that the groundwater heads show pattern similar to calibration and validation periods. However, the groundwater head rises to 2 m in the southern part of the study area, and 0.2–0.5 m in the middle and northern part of the district. The predicted heads for different stresses are shown in Figure 15.

5.   Conclusions

The study investigates the surface water-groundwater interactions and impact of reduced draft on groundwater resources in a part of central Gangetic plain. The alluvial aquifer in the study area (Bhojpur district, Bihar) exhibits depth to groundwater in the range 3 to 9 m bgl in the pre-monsoon season during the year 2018. The trend in groundwater is found to be insignificant at annual scale as revealed by MK-test, however the groundwater levels vary considerable between seasons. A considerable rise during monsoon is observed due to recharge from rainfall and inflow from river, while the groundwater utilization during post-monsoon season leads to declined groundwater level. The dominating land uses in the study area are agricultural land (64.07%), followed by built-up area (24.75%), barren land (5.65%), waterbodies (2.31%), shrub (2.13%), and vegetation (1.09%). A transient groundwater model was developed for the study area. The results show that the groundwater head is higher in the southern part as compared to the northern part of the study area and GW is contributing to the river Ganga. The effect of better agricultural water management practices is simulated by reducing the groundwater extraction by 10%. It revealed that such reduction in groundwater withdrawal could help sustainable management of aquifers. The simulation results showed that a reduction of groundwater withdrawal by 10% would lead to a rise of up to 2 m in the southern part and rise of 0.2–0.5 m in the middle and northern parts of the district. It may be suggested that an integral agricultural water management practices, including increased water use efficiency, cropping pattern and rainwater harvesting, should be adopted for sustainable use of groundwater resources.

Acknowledgments

The authors thankfully acknowledge the data supports from CGWB, Ministry of Jal Shakti, Department of Agriculture, Govt. of Bihar and Minor Water Resources Department, Govt. of Bihar. The data and information obtained from different reports published by various organizations on this area are also duly acknowledged. We are also thankful to the Director, National Institute of Hydrology, Roorkee for providing all administrative and financial support for the successful completion of the study.

Conflict of interest

The authors declare there is no conflict of interest.