Evaluating groundwater depletion under natural and induced stresses: a numerical modeling approach toward aquifer sustainability

Groundwater is the largest source of freshwater on Earth. Its widespread availability, low capital cost, and dependability make it an ideal source of water supply for domestic, agricultural, and industrial needs (Omar et al. 2020; Gebere et al. 2021; CGWB 2022). The misconception of pumping rates not exceeding natural recharge rates for safe groundwater development has led to the overexploitation of groundwater resources (Zhou 2009; Demiroglu 2019; Gaur et al. 2023).

Quantification of the water budgets offers a comprehensive way to evaluate the availability and sustainability of water supply and allows for the investigation of natural or human-induced impacts on groundwater storage (Healy et al. 2007). Distinct approaches have been applied to quantify groundwater flow budgets, such as the water balance methods (Yang 2022; Bayat et al. 2023), remote sensing, and geographic information system (GIS)-based methods (Rashid & Ahmed 2018; Gaur et al. 2021; Nazari et al. 2023). Numerical groundwater models have been found to be one of the most effective methods in groundwater resources quantification and flow budget assessment studies, as they can effectively account for aquifer heterogeneity and spatiotemporal variability in source and sink data (Zhou & Li 2011; Omar et al. 2021a; Ma et al. 2023). A wide range of literature is available on application of numerical models in groundwater budget studies (de Graaf et al. 2017; Viaroli et al. 2018; Alattar et al. 2020; Omar et al. 2021b; Bayat et al. 2023). Yet, only few studies have discussed the surface and groundwater budget dynamics as a single system and the implications of natural and induced stresses on other budget components.

Several authors have utilized integrated models to quantify groundwater budget components and assess the impact of natural and induced stresses. Sisay et al. (2023) quantified the flow budget of Modjo River catchment in Central Ethiopia using the Soil and Water Assessment Tool (SWAT) and a steady-state MODFLOW-NWT model and observed that decrease in recharge affected the baseflow to rivers more severely than increase in pumping. Mohsenifard et al. (2023) applied the SWAT and MODFLOW 2000 model to study the influence of land use and cropping change on groundwater of Shazand plain, Iran. Nasiri et al. (2022) determined water balance components and groundwater–surface water interactions of Samalqan plain, Iran, using SWAT and MODFLOW model. Gebere et al. (2021) employed a steady-state MODFLOW model to quantify flow budget of Modjo River catchment, Central Ethiopia, and observed that decrease in recharge rates had more influence on groundwater contribution to rivers and baseflow than increase in abstraction rates. Eltarabily et al. (2018) assessed the impact of different recharge and discharge scenarios on Quaternary aquifer east of Nile Delta, Egypt, using a MODLFOW model and observed that the current rate of abstraction was sustainable, and an increase of 50% pumping would result in reversal of stable conditions. Demiroglu (2019) suggested that groundwater budget studies needed to be reassessed with new insights that consider all the budget components. However, most of the studies were focused on impact of increase/decrease in recharge and pumping stresses on groundwater levels and river–aquifer interactions and failed to account for other budget components, especially storage. Further, the spatiotemporal variations in storage dynamics have been rarely discussed.

In this study, SWAT and MODFLOW-NWT model was employed to quantify the water budget dynamics and assess the impact of recharge and abstraction stresses on groundwater budget components by taking up the case study of Lower Betwa River Basin (LBRB), India. LBRB is part of the Betwa River Basin (BRB), an important drainage basin in the Bundelkhand region of central India. The region has been affected by frequent drought occurrences in the recent past (Gupta et al. 2022) and the severe meteorological and hydrological drought of 2005–2007 season resulted in drying up of 70% of the tanks, ponds, and dug wells in the region (Gupta et al. 2014). Socioeconomically poor, the lack of water for irrigation has caused food insecurity and mass migration in the region (Gupta et al. 2014; Chaurasia & Chandra 2021). Previous studies in the region have evaluated the climate change impact (Gupta et al. 2023), the potential recharge zones (Pandey et al. 2021), groundwater potential zones (Jeet et al. 2019) using machine learning (Kumar et al. 2021), rainfall trend and drought jeopardy (Gupta et al. 2022), groundwater-based crop water productivity (Bhattacharjee et al. 2021), and surface water balance (Suryavanshi et al. 2017; Desai et al. 2021). However, no literature was available on the estimation of groundwater balance and groundwater depletion in the region, to the best of the authors' knowledge. Thus, the specific research objectives of this study are as follows: (1) to quantify regional groundwater flow budget dynamics of the LBRB from 2003 to 2020. (2) Investigate trends in groundwater levels over the 17-year study period. (3) Demarcate groundwater storage stressed zones to investigate management scenarios for sustainable groundwater development.

This study is focused on the LBRB, which is typically underlain by alluvial and weathered formations covering approximately 15,700 km² area. The LBRB is situated between 110 and 400 m above mean sea level (amsl), and distinctly divided from the upper BRB (358–713 m amsl) (Kumar et al. 2024). Three major districts, namely, Jhansi and Lalitpur in Uttar Pradesh and Tikamgarh in Madhya Pradesh, cover most of the LBRB (Figure 1). Agriculture and barren land occupy over 70% of the basin, while only 4.5% was under built-up area. The area is classified as having semi-arid climate with the potential evapotranspiration rates of over 1,600 mm (CGWB 2017a, b, c). Groundwater recharge is primarily through precipitation infiltration, irrigation return flow, and seepage from canal and surface water structures (Joshi et al. 2021). The general groundwater flow direction was from southwest toward the northeast.

Because the LBRB is geologically a hard rock area, the groundwater resources are limited. The major part of the study area is rain-fed, socioeconomically poor, and lacks access to infrastructure and technology (CGWB 2017a). Agriculture and livestock rearing are the main occupations in the region and are dependent on monsoon rainfall. To combat the uncertain and erratic monsoon, multiple dams and canals were constructed. However, during the drought periods these reservoirs were not able to meet the water requirements (CGWB 2017a, b, c). Per the 2011 census, the total population of the three districts was 46.65 lakhs. With 70% of the land under agricultural use, farming is the main contributor to the economy in the LBRB, and the cropping intensity stood at 165% (CGWB 2017a, b, c). According to the latest Dynamic Groundwater Assessment Report (CGWB 2022), the annual groundwater demand in the three major districts within the study area was at 1,055 million cubic meters (MCM), i.e. the average water demand in the basin was 76.56 mm/year, majority of which was allocated to irrigation (90%) and rest to domestic (9%) and industrial use (1%). The dynamic groundwater recharge was estimated at 1,376 MCM and the average stage of groundwater extraction in the three districts worked out to be 77% (CGWB 2022). In accordance with the stage of groundwater extraction, the LBRB aquifers were under semi-critical condition (CGWB 2022).

Data accessibility is a key aspect in the model development phase (Omar et al. 2023). The hydro-meteorological, hydrogeological, groundwater level, and groundwater demand data were obtained from different sources and are presented in Table 1. The lithological data for borewells were obtained from Central Ground Water Board (CGWB) office, Lucknow, and the National Project on Aquifer Management (NAQUIM) reports. The hydrogeological information of aquifers was collected from NAQUIM reports (CGWB 2017a, b, c). The groundwater level data were obtained from the India WRIS portal (https://indiawris.gov.in/wris/#/groundWater). The district-wise groundwater draft data for reference were collected from the Dynamic groundwater resources report (2004–2022) (CGWB 2022).

The watershed scale, semi-distributed, physically based, continuous hydrologic, and water balance model the SWAT (Arnold et al. 1998) was employed to simulate the surface water balance of the BRB. The SWAT model simulates the water balance based on two phases, the land phase that controls the amount of water, sediment, and nutrients being transported to the main channel and the channel phase that moves all of them through the channel network to the outlet (Neitsch et al. 2011). Discretization of the watershed is done based on the water network into sub-watersheds, which are further segmented into hydrologic response units (HRUs) based on identical land use, slope, and soil characteristics.

ArcSWAT 10.5 was employed to simulate the BRB hydrology. The input data details such as topography, land use, soil characteristics, meteorological data, and their source are provided in Table 1. The digitized stream network was burned into DEM with final outflow point serving as the drainage point. Monthly gridded precipitation and temperature data for 18 years (2003–2020) were obtained from the Indian Meteorological Department (IMD). The model was calibrated and validated using streamflow data of two Central Water Commission (CWC) gauge stations (Figure 1). The BRB was discretized into 182 sub-basins of which 78 were in the LBRB. The Sequential Uncertainty Fitting algorithm application (SUFI-2) pre-embedded in the SWAT-Calibration and Uncertainty Program (CUP) was used for model sensitivity, calibration, and validation. The SUFI-2 algorithm utilizes the highly efficient Latin Hypercube sampling scheme to obtain optimal results and is capable of handling large number of parameters with computational efficiency (Desai et al. 2021). The initial parameters and their range were based on previous studies in the BRB (Suryavanshi et al. 2017; Desai et al. 2021).

The conceptual flow model for the LBRB aquifer system was defined based on the available hydrological and hydrogeological information of the area (Figure 2(b)). The topographic and morphological information was retrieved from SRTM DEM data. Aquifer geometry was delineated based on borewell exploration data obtained from CGWB office, Lucknow. Hydrogeology was defined based on CGWB (2012) report on Principle Aquifer Systems India and refined using NAQUIM reports (CGWB 2017a, b, c). The linear quartz reefs in the field have been conceptualized as horizontal flow barriers (HFBs) in the model with very low K values as they prevent groundwater flow (CGWB 2017a). Initial values of hydraulic conductivity (K) and specific storage (S_s) were obtained from pumping data (CGWB, Lucknow), NAQUIM, and specific yield (S_y) values from Groundwater resource Estimation Committee (GEC) 2015 (CGWB 2017d). The hydrogeological boundaries of the basin were defined as no-flow boundary toward northwest, east, west, and south directions, based on the topographic highs on the surface that coincides with the groundwater divide below (Jain et al. 2021). Specified Flow boundary (SFB) was applied in the north, northeast, southeast, and southwest directions, estimated based on groundwater contour data (Tewari et al. 2024). GIS environment was used to organize, analyze, and synthesize the data collected from field and other sources.

Groundwater recharge is a very important boundary condition in the model input, as it has the greatest influence in the groundwater budget, but is also a difficult parameter to estimate (Healy et al. 2007). Its spatiotemporal variations are influenced by many factors such as land use, vegetation, hydrogeological conditions, and climate (Healy et al. 2007; Karki et al. 2021). SWAT output has been successfully implemented in the past, to quantify groundwater recharge into aquifers (Chinnasamy et al. 2018; Karki et al. 2021). The HRUs-based groundwater recharge estimates were averaged on each sub-watershed and incorporated into the MODFLOW model as specified flux using the recharge (RCH) package (Anderson et al. 2015). Similarly, averaged evapotranspiration (ET) rates for each sub-watershed were incorporated into the MODFLOW model using the evapotranspiration (EVT) package with a constant extinction depth of 3 m, similar to Yao et al. (2015).

Based on the lithological information, the LBRB aquifer was conceptualized as a two-layered system. The first layer was composed of clay and kankar toward north, weathered granite in the central region and non-porous sandstone toward the south (Figure 2(a)). The first layer was delineated as an unconfined aquifer with smoothed SRTM DEM as the top of the first layer and the thickness varied from 7 to 75 m. The second layer consisted of fractured fine to medium chips of granite, diorite, or dolerite (CGWB 2017b) and conceptualized as a semi-confined/convertible layer based on the concept of equivalent porous medium (EPM) (Anderson et al. 2015). The second layer varied from 0 to 40 m, taking the total thickness of the aquifers to a maximum of 115 m.

The model was discretized into 212 rows and 175 columns resulting in a total of 15,757 active cells in one layer. A transient model was built from 2003 to 2019 with quarterly stress periods. Initial values of aquifer properties such as K, S_s, and S_y were obtained from NAQUIM reports (CGWB 2017a, b, c) and adjusted during model calibration. The transmissivity values ranged from 12.87 to 145 m²/day and storativity values ranged from 2.8 × 10⁻⁴ to 1.57 × 10⁻³ (CGWB 2017a, b, c). Different K zones were deciphered based on the hydrogeological and thickness distributions in the study area (Figure 2(a)). Specific yield (S_y) values were adopted from GEC-2015 norms; the detailed description of the calibrated aquifer properties is given in Table 3.

Sensitivity Analysis was conducted to determine response of various parameters. Review of the previous literatures in the BRB (Suryavanshi et al. 2017; Desai et al. 2021) was used to select the 10 most sensitive parameters for calibration. The Sequential Uncertainty Fitting ver.2 (SUFI-2) in SWAT-CUP was used to determine sensitivity of the parameters, and calibrate and validate the model (Abbaspour 2015). The summary of the sensitivity analysis is provided in Table 2 and the detailed analysis for the same is available in Kumar et al. (2024).

The performance of the hydrological model was evaluated using the widely accepted Nash–Sutcliffe efficiency coefficient (NSE) and the coefficient of determination (R²). NSE objective function was initially optimized using mean monthly discharge data of Basoda and Mohana gauging stations from 2003 to 2013 and the results were validated for the 2014–2020 period (Table 3). NSE indicates the strength of the relationship between observed and simulated streamflow values while R² indicates the extent to which the model can account for the overall variation in the observed discharge data (Kausher et al. 2024). The NSE and R² values varied from 0.65 to 0.77 and 0.71 to 0.79 in the calibration phase, indicating a good to very good performance (Moriasi et al. 2007), while the NSE and R² values for the validation phase varied from 0.74 to 0.86 and 0.75 to 0.90, respectively, which indicated a very good model performance (Moriasi et al. 2007). The P factor indicates the percentage of data bracketed by 95% prediction boundary (95PPU) and R factor measures the thickness of the 95PPU model envelope (Kausher et al. 2024). A P factor of 0.61 indicates that 61% of the observational data is enclosed within the 95PPU (Kausher et al. 2024).

Initial manual calibration involved the adjustment of K, S_s, and streambed conductance for the STR and DRN package. The recharge and discharge (through pumping wells) values were corroborated with the NAQUIM reports (CGWB 2017a, b, c). After the initial process, K and S_s parameters were chosen for automated calibration. PEST is a widely supported tool for parameterized model calibration (Doherty & Hunt 2010). Tikhonov regularization was used to provide PEST with sufficient information to obtain an inverse solution. The broad range of K and S_s values obtained from pumping test reports of CGWB, Lucknow, and the NAQUIM reports was used to constrain the parameter estimation process while avoiding unrealistic parameter distributions. Further, Truncated Singular Value Decomposition (SVD) with parallel PEST was used to reduce the time necessary to complete the optimization process. The fitted aquifer properties after calibration are given in Table 4.

A positive (negative) value of Z_MK indicates that the data tend to increase (decrease) with time, at the Type I error rate, where 0 < α < 0.5. Note that α (significance level) is the tolerable probability that the MK test will falsely reject the null hypothesis. Then H₀ is rejected and H_a is accepted if Z_MK ≥ Z_1−α, where Z_1−α is the 100(1 ₋α)th percentile of the standard normal distribution.

Spatial trends showed that the western region of Tikamgarh and Chattrapur districts received the maximum amount of average annual rainfall and correspondingly generated larger values of surface runoff and groundwater recharge (Figure 4). The ET rates were maximum for regions with large water bodies. Overall, the northern region of Jhansi district received the least amount of rainfall, and the ET rates were relatively higher and accounted for poor surface runoff and groundwater recharge.

For the groundwater balance estimation, the simulation model was run for 69 stress periods using a 4-month time step each year (January 2003 to January 2020). The universal groundwater budget equation (Flow In–Flow Out = Change in Storage) terms are referred as groundwater components and written in terms of flux (volume per unit time) in this paper. The simulated groundwater flow budget components are tabulated as annual average volumetric flux (Mm³/year) (Table 5). For transient simulation the storage component is given as loss in inflow and gain in outflow in accordance with storage loss as positive and gain as negative (Jones & Torak 2006). Simulation results displayed that for the entire modeled area and study period, areal recharge was the largest inflow (64%) into the system, followed by river leakages (19%) and negligible lateral inflow from SFB (0.4%). Simulated outflow was largest for abstraction wells (41%), followed by the network of minor streams outflow (21%), discharge into rivers (14%), ET (8%), and SFB (0.6%). The difference between the inflow and outflow was accounted by a net cumulative loss in groundwater storage at the rate of 36.5 Mm³/year. This was at about 1.4% of the total budget. It is important to note that even when the average stage of groundwater development was at only 64% during the study period, the LBRB aquifers were gradually losing water from storage. The observation highlights the importance of critical analysis of groundwater flow budgets and that natural recharge alone cannot determine the sustainable yield in aquifers (Zhou 2009).

Analysis of the results indicated that during the 2003–2008 period, the majority of the LBRB experienced a decline in the groundwater levels with 73% of the basin under 0–1.2 m/year decline rate (Figure 7(a)). Additionally, about 6% of the LBRB had an average annual decline rate of 1.2–1.6 m and was concentrated on the north and central regions of the LBRB. The linear quartz reef acted as groundwater barrier (CGWB 2017a) and its influence was visibly observed with the right side of the HFB suffering greater decline in groundwater levels (Figure 7(a)). Such geological heterogeneities can be potential subsurface water storage structures (Senthilkumar & Elango 2011; Houben et al. 2018). The HFB along with model boundary was responsible for groundwater table rise in the west region of the LBRB. With the decreasing trend in aquifer recharge rates, over 63% of the basin showed a decreasing trend in groundwater levels, while 11% had an increasing trend, and the rest showed no significant trend (Figure 7(e)).

During the phase of increasing recharge patterns (2008–2014), over 77% of the LBRB responded with increasing trend in groundwater levels (Figure 7(f)). Groundwater head response rate in the south was higher than for the rest of LBRB. However, the northern blocks of Jhansi district with alluvial aquifers still showed a declining trend of 0 to −0.6 m/year (Figure 7(b)). The last phase of mixed recharge patterns (2014–2020) showed intensive decline rate (greater than 1.2 m/year) in the north and central regions of the LBRB (Figure 7(c)). While the north and central region groundwater levels declined, the rest of the LBRB showed no significant trend (Figure 7(g)).

Overall, greater than 1.2 m/year decline rates were observed only in grid cells occupied by extraction well points in the north and central regions (Figure 7(d)). The north and central regions that experienced larger decline rates (0.6–1.2 m/year) were the most affected regions in the LBRB. Over 53% of the basin was under 0–0.6 m/year decline rate, while 30% of the basin showed no significant trend in groundwater level change and these regions occupied the river stretch, the west, and south sections (Figure 7(h)). The over 1.2 m/year head rise observed in the western region of the LBRB was due to the HFB and model boundary conditions, which restricted the flow of the recharged water.

Based on the average storage change over the entire study period, a groundwater stressed zone map was prepared for the whole of LBRB (Figure 8(e)). The basin was classified into five categories ranging from very low (−419 m³/d) to very high (543 m³/d) storage stressed zones. It was observed that the majority of the basin area (78%) was under positive or low storage stress (–43 to 8 m³/d), which indicated that the aquifer systems were stable under the existent stress conditions (Eltarabily et al. 2018). In the north region, about 16% of the LBRB was gradually losing storage water and was under moderate storage stress (9–95 m³/d), most of it in the Jhansi district. The very high storage stressed zone was concentrated in the north region of Bamaur block in the Jhansi district and coincided with the alluvial zone of the LBRB aquifer (Figure 2(a)). It was about 1.6% of LBRB area and was losing storage water at the rate of 268–543 m³/d. Groundwater depletion was also maximum in the same area (Figure 7(d)) and it was the most vulnerable area in the whole LBRB, requiring immediate management plans for conservation and sustainable development.

Quantification of the flow budget of any assessment unit is crucial for effective water resource management planning. It provides a basis for assessing how natural or anthropogenic induced changes in one part of the cycle affect the other (Healy et al. 2007). A thorough investigation and quantification of flow budget dynamics for the 17-year period resulted in comprehending the spatial storage loss and groundwater depletion in the LBRB system. The northern blocks of Jhansi and the central region of Tikamgarh districts were the most affected region in the whole LBRB (Figures 7(d) and 8(e)). Two blocks in the Jhansi district under very high and moderate groundwater stress (Bamaur and Gursarai) (Figure 8(e)) were selected for the aquifer management planning. Zone budgeting was employed to quantify flow budget within the blocks. Zone budgeting allows for analysis of flow budget for a specific area of interest within the modeled area. The net annual average flow budget (inflow–outflow) for the selected blocks is presented in Table 6.

The annual average zone budgeting of Bamaur and Gursarai blocks for the 17-year study period revealed that the abstraction rates had been exceeding the natural recharge rates in both the blocks (Table 6). The Bamaur block with an area of 799 km² received about 30% more groundwater recharge than the Gursarai block, which had an area of 705 km², but the abstraction rates were also significantly higher (27%) in Bamaur block than the Gursarai block. Both agriculturally intensive blocks with cropping and irrigation intensity of 150 and 115%, respectively, were dependent on groundwater during the dry summer and winter seasons (CGWB 2017a). The groundwater contribution to irrigation in 2013 was 30 and 49% in Bamaur and Gursarai blocks, respectively, and the rest was through canals and ponds (CGWB 2017a). The intense requirement of groundwater to sustain crop productivity under poor aquifer replenishment due to low recharge (Figure 4) led to the decline in groundwater storage and water table. The resultant annual average storage loss was observed to be 14.1 and 6.3 Mm³/year in Bamaur and Gursarai blocks, respectively. The seasonal analysis of the two blocks showed similar trends in storage variation as the whole LBRB aquifer system (Figure 6). While other budget components showed similar variations in both the blocks, the river and boundary exchanges varied, with the Bamaur block that was bordered by rivers on three sides (Figure 8(e)) losing water to the rivers, and the aquifers drawing water from the river in the Gursarai block.

When abstraction rates exceed total recharge, the groundwater storage will be continuously used to balance the extraction rates and the groundwater levels will continue to deplete, which is also known as groundwater mining (Zhou 2009). The Bamaur and Gursarai blocks showed similar mining conditions and, subsequently, either the pumping rates must be reduced, or additional recharge must be identified to prevent further groundwater depletion and promote storage recovery. Two management scenarios were analyzed to assess the impact of natural recharge and abstraction stress on other budget components. First, the abstraction rates of Bamaur and Gursarai blocks were gradually reduced by 10% intervals till 50% of actual extraction. In the second scenario, the natural recharge in the two blocks was increased by intervals of 10% till 50% and the resultant change in budget components was recorded.

In the case of incremental increase in recharge rates, for every 10% increase in natural recharge the storage loss decreased by an average of 6,391 m³/d or 17% in the Bamaur block (Figure 9(c)). With the increase in inflow, the outflow through other budget components also increased, and it was the maximum for river exchanges. A similar observation was reported by Sisay et al. (2023) with regard to the Modjo River catchment, central Ethiopia, where reduction in recharge rates drastically reduced the baseflow into rivers. Due to the lower water table in Bamaur block no significant change in ET rates was observed and the low boundary inflow rates switched to outflowing conditions. In the Gursarai block, the 10% incremental increase in natural recharge rates resulted in 501 m³/d or mere 3% reduction in storage loss (Figure 9(d)). Most of the increased inflow of water was balanced off by decrease in inflow from boundary and river exchanges. Slight increase in ET rates was also observed, but it was not as significant as the change due to abstraction reduction.

The heterogeneity in response of budget components in the two blocks can be due to other groundwater head influencing factors such as aquifer properties, slope, and aquifer thickness (Tewari et al. 2024). Over half of the Bamaur block was made of alluvial aquifers with higher K and S_y rates. Higher Sy rates influenced higher changes in storage component (Nazari et al. 2023). Gursarai block was underlain mostly by weathered formations with low K and S_y properties (Figure 2(a)).

Similar observations were recorded by Sisay et al. (2023) and Gebere et al. (2021); when they increased abstraction rates in Mojdo aquifers, Ethiopia, the baseflow to rivers decreased. However, the reduction in baseflow contribution was more significant when the recharge rates were lowered indicating the sensitivity of Mojdo aquifers to recharge rates. In the case of LBRB aquifers, abstraction rates had a larger influence than recharge rates on flow budget components, especially storage, yet practically, it would be impossible to reduce pumping demands by over 50% without securing an alternative source of supply water. Techniques such as micro irrigation through drip and sprinkler irrigation can enhance water use efficiency and lower the water demand and they are currently being adopted in the region (CGWB 2017a). About 10–20% reduction in abstraction could be achieved using the above techniques, which could lower the storage loss by 29 and 16% in Bamaur and Gursarai, respectively. Further improvement in storage loss reduction could be achieved by targeted interventions such as managed aquifer recharge, which is more effective in managing groundwater depletion (Dillon et al. 2019). Practices such as rainfall-runoff harvesting have been successful in improving groundwater levels and increasing cropping area in the Bundelkhand region (Garg et al. 2020; Singh et al. 2021).

All groundwater flow models are limited by the conceptualization of their hydrologic system, assumptions taken, and uncertainty of how the system is represented (Jones et al. 2017). The lack of any previous modeling studies in the area and sufficient data on aquifer properties resulted in conduction of multiple model runs with various conceptualizations to arrive at an acceptable solution. Lack of data on reservoir and lake water level variations led to their exclusion from the LBRB groundwater flow model. Information on irrigation return flow could have improved estimated results of groundwater recharge. Zonal distribution of aquifer hydraulic properties gave abrupt storage variation results and would require detailed field investigation to improve spatial distribution of aquifer properties. Lack of data on smaller rivers of the LBRB limited their conceptualization to drain package in the model, restricting estimation of their exchanges.

This study examined the regional water balance in the LBRB using integrated surface water and groundwater flow models. The study area showed spatial variations in quantified budget components, influenced by the precipitation and groundwater recharge rates. Groundwater storage variations are a seasonal phenomenon in the LBRB, with the aquifers gaining water in storage during the post-monsoon and winter seasons and losing it during the summer and monsoon seasons. The overall change in storage during the 17-year study period was found to be negative with the LBRB aquifers losing 36.5 Mm³/year. The storage loss was more severe toward the northern Jhansi district of the LBRB where groundwater recharge was lowest. Zone budget analysis revealed that the Bamaur and Gursarai blocks in Jhansi district were losing storage water at the rate of about 0.11–2.1 Mm³/year. Groundwater trend analysis revealed the influence of linear quartz reefs in distribution of groundwater head in the region, with most of the LBRB showing a declining trend. The management scenario analysis showed that changes in abstraction rates had a greater influence on storage depletion than changes in natural recharge rates. However, to achieve aquifer sustainability, economic and social sustainability also requires consideration, and the supply (abstraction) cannot be reduced at free will. Management practices such as micro irrigation and conjunctive use can reduce the demand to some extent but needs to be further supported by managed interventions such as managed aquifer recharge to maintain aquifer sustainability in the LBRB region.

The application of loosely coupled integrated surface and groundwater models was effective to decipher the long-term water budget dynamics and analyze suitable management scenarios toward groundwater sustainability in the LBRB. Future studies can concentrate on storage stressed regions by development of more complex local-scale models based on the basin-scale model results. The scope of artificial recharge can be assessed using the local-scale models and an effective management plan can be prepared for sustainable development of groundwater resources in the Bundelkhand region.

The authors would like to express their gratitude to the National Resources Data Management System (NRDMS) (Department of Science & Technology, New Delhi-110016 Ref: NRDMS/01/259/018(G)) for their support. The authors would also like to acknowledge CWC (YBO), New Delhi, and CGWB, Northern Region (Lucknow), for provision of the requisite data.

A.T. contributed to conceptualization, modeling, and original draft preparation; P.K.S. contributed to supervision, review, and editing; S.G. contributed to conceptualization, review, and editing; R.K. contributed to modeling, data curation, and preparation of results; S.M. contributed to conceptualization and preparation of results.

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

All the authors have no conflict of interest.