Hydrological Modeling of Aquifer’s Recharge and Discharge Potential by Coupling WetSpass and MODFLOW for the Chaj Doab, Pakistan

: The estimation of the groundwater (GW) potential in irrigated areas is crucial for the sustainable management of water resources in order to ensure its sustainable use. This study was conducted in a selected area of the Chaj doab, Punjab, Pakistan, to quantify the impacts of the pumping and the recharge on the aquifer therein. To that end, a groundwater ﬂow model (MODFLOW) and a groundwater recharge model (WetSpass) were coupled to assess the conditions of the aquifer. The model was calibrated manually on twelve-year data (2003–2014) against the observed groundwater levels, and it was validated with ﬁve-year data (2015–2019). Three main scenarios (divided into ten subscenarios) were simulated for the future prediction of the groundwater: Scenario-I (to assess the impact of the pumping if the prevailing conditions of the years from 2003 to 2019 were to continue until 2035); Scenario-II (to assess the impact of the pumping on the aquifer by increasing the pumping capacity by 25, 50, 75, and 100% for the coming 10 years); and Scenario-III (to assess the impact on the aquifer of the decrease in the average groundwater recharge from the river by 50% by following the same pumping trend). The Scenario-I results show that there would be an 18.1 m decrease in the groundwater table at the end of the year 2035. The Scenario-II results predict decreases in the water table by 2.0, 5.5, 9.8, and 14.3 m in the year 2029 as a result of increases in the pumping capacity of 25, 50, 75, and 100%, respectively. The results of Scenario-III show that, with the decrease in the recharge from the rainfall, there would be a 0.7 m decrease in the water table, and that, from open-water bodies, there would be a 2.4 m decrease in the water table. These results are very helpful for determining the recharge and discharge potential of the aquifer.


Introduction
Water availability is essential for meeting agricultural, domestic, and industrial needs.
In the past few decades, the burden on the groundwater resources of Pakistan has increased because of their unchecked, unguided, and unregulated exploitation [1]. At a global level, about 750-800 billion cubic meters (BCM) of groundwater is being consumed for agriculture [2,3]. Shortfalls in the surface water is the reason for the increasing use of groundwater [4]. The supply of surface water has decreased significantly in recent years, and the demand for groundwater has also increased proportionately [5]. Groundwater resource management in the Chaj doab, Pakistan, and especially in the Khadir canal subdivision, is necessary for improving agriculture and protecting the ecosystem, as well as for biodiversity [6]. condition, the layers, the starting and end times, the digital elevation of the model/study area layout, the climatic parameters (i.e., rain, ET, sunshine hours, wind speed), the recharge, the discharge, the groundwater level, and the grid size. However, on the basis of the provided initial conditions, the conceptual model is converted into the numerical model by applying the values of all of the parameters in the numerical model, "MODFLOW". A pictorial elaboration of the conceptual model is shown in Figure 2.  There is a gentle slope of 0.3 m/km, and an elevation range of 155-255 m throughout the majority of the region. During the months of June through October, the west monsoon is chiefly responsible for the region's significant rainfall. The region receives 65 to 70 percent of its total rainfall in the monsoon season, which ranges from roughly 1000 mm towards the northeast, to 230 mm towards the southwest of Punjab [32].

Conceptual Modeling Approach
The formulation of a conceptual model is a very crucial step for the development of the numerical model for the purpose of a modeling study. The intention of the conceptual model development is helpful for the simplification of the actual model. A successful conceptual model should be very close to the actual field conditions of the system. It is difficult to study actual systems because of the large number of complexities, and so the conceptual model is very important for the simulation of the actual field conditions. A good conceptual model carries less chances of errors in the numerical models. The designed conceptual model is based on the initial conditions of the model, the boundary condition, the layers, the starting and end times, the digital elevation of the model/study area layout, the climatic parameters (i.e., rain, ET, sunshine hours, wind speed), the recharge, the discharge, the groundwater level, and the grid size. However, on the basis of the provided initial conditions, the conceptual model is converted into the numerical model by applying the values of all of the parameters in the numerical model, "MODFLOW". A pictorial elaboration of the conceptual model is shown in Figure 2.

Climatic Data
The input data of the climatic parameters for the WetSpass-M were taken from the Pakistan Meteorological Department from 2000 to 2019. The minimum and maximum long-term annual rainfalls were 424 and 679.5 mm, respectively, with a mean value of 455.7 mm. A notable fluctuation in the long-term temperature was observed in the summer (April-September), with a temperature range between 19 and 41 °C , and, in the winter season (October-March), with a range between 4 and 28°C. The long-term monthly wind speed data show that the maximum wind speed was recorded on July 9 at 674 km/h, and that the minimum wind speed was 1.02 km/h.

Lithology
The lithological data for the study area were obtained from the Water and Power Development Authority (WAPDA). A comprehensive study was conducted by the Water and Soil Investigation Division [33] to determine the subsurface lithology and the geological characteristics of the Chaj doab, Pakistan. Only sixteen boreholes (Figure 1d) in the study area were found. The sampling depths of the data, on average, were taken at 20 cm. The sampling results show that the lithology was comprised of fine (silty clay, clay, and sandy clay) and moderately fine soil (clay loam, sandy clay loam, and silty clay loam), as is shown in (Figure 1c).

Topography
A map of the topography of the study area was attained from the Shuttle Radar Topography Mission (STRM-1), with a 30 m resolution, and by using the link: https://earthexplorer.usgs.gov/ (accessed on 24 October 2020). The elevation ranged between 155 and 255 m at the lowest and highest points, with an average elevation of 190 m from the mean sea level (Figure 1b).

Groundwater Data
The groundwater level data were collected from 19 piezometers from the Punjab Irrigation Department, the "Sargodha Irrigation Zone", for the period from 2003 to 2019. The groundwater level data were collected two times yearly: one time during the premonsoon (June), and the second time during the postmonsoon (October), and these data were available. The mean groundwater depth fluctuation in the area for the period from 2003 to 2019 was recorded as 5.32 m, as is shown in Figure 3.

Climatic Data
The input data of the climatic parameters for the WetSpass-M were taken from the Pakistan Meteorological Department from 2000 to 2019. The minimum and maximum long-term annual rainfalls were 424 and 679.5 mm, respectively, with a mean value of 455.7 mm. A notable fluctuation in the long-term temperature was observed in the summer (April-September), with a temperature range between 19 and 41 • C, and, in the winter season (October-March), with a range between 4 and 28 • C. The long-term monthly wind speed data show that the maximum wind speed was recorded on July 9 at 674 km/h, and that the minimum wind speed was 1.02 km/h.

Lithology
The lithological data for the study area were obtained from the Water and Power Development Authority (WAPDA). A comprehensive study was conducted by the Water and Soil Investigation Division [33] to determine the subsurface lithology and the geological characteristics of the Chaj doab, Pakistan. Only sixteen boreholes (Figure 1d) in the study area were found. The sampling depths of the data, on average, were taken at 20 cm. The sampling results show that the lithology was comprised of fine (silty clay, clay, and sandy clay) and moderately fine soil (clay loam, sandy clay loam, and silty clay loam), as is shown in (Figure 1c).

Topography
A map of the topography of the study area was attained from the Shuttle Radar Topography Mission (STRM-1), with a 30 m resolution, and by using the link: https: //earthexplorer.usgs.gov/ (accessed on 24 October 2020). The elevation ranged between 155 and 255 m at the lowest and highest points, with an average elevation of 190 m from the mean sea level (Figure 1b).

Groundwater Data
The groundwater level data were collected from 19 piezometers from the Punjab Irrigation Department, the "Sargodha Irrigation Zone", for the period from 2003 to 2019. The groundwater level data were collected two times yearly: one time during the premonsoon (June), and the second time during the postmonsoon (October), and these data were available. The mean groundwater depth fluctuation in the area for the period from 2003 to 2019 was recorded as 5.32 m, as is shown in Figure 3.

Groundwater Recharge and Evapotranspiration
The spatial scattering of the evapotranspiration and the groundwater recharge was taken from [34]. The WetSpass-M [35] is widely applied to evaluate the distribution of the groundwater balance component model on a spatial basis [35][36][37][38][39][40]. The WetSpass-M model takes the spatial distributions of the soil, the slope, the groundwater depth, the land use, and the climatic data. The seasonal water balance distribution was considered to be: P stands for precipitation; S is used to indicate the surface runoff; ET is for evapotranspiration; and R denotes the groundwater recharge. The runoff, recharge, interception, and evapotranspiration are all included in the WetSpass-M ArcGIS integrated model as input parameters. The weather/climatic data (potential evapotranspiration, temperature, wind speed, and rainfall), the terrain, the soil texture, the land use/land cover, the leaf-area index data, and the groundwater depth are all basic inputs for this model. It also takes into account a number of other variables [14,28]. ArcGIS was used to prepare all of the input data, which were then imported into the model with the same cell size in ascii format. To prepare the input data from 2000 to 2019, ArcGIS 10.3.1 was used, and the raster cell size was maintained at 100 m × 100 m. Each parameter has a total of 72,713 grids. WetSpass is based upon the principle of the water balance. As a result, the average annual recharge from the rainfall was between 99.10 and 177.2 mm, with a mean value of 119.4 mm per year, as is displayed in (Figure 4a). The long-term annual evapotranspiration values for the years from 2000 to 2019 ranged between 134.4 and 287.2 mm/year, with a mean value of 231 mm/year, as is depicted in (Figure 4b).

Groundwater Recharge and Evapotranspiration
The spatial scattering of the evapotranspiration and the groundwater recharge was taken from [34]. The WetSpass-M [35] is widely applied to evaluate the distribution of the groundwater balance component model on a spatial basis [35][36][37][38][39][40]. The WetSpass-M model takes the spatial distributions of the soil, the slope, the groundwater depth, the land use, and the climatic data. The seasonal water balance distribution was considered to be: P stands for precipitation; S is used to indicate the surface runoff; ET is for evapotranspiration; and R denotes the groundwater recharge. The runoff, recharge, interception, and evapotranspiration are all included in the WetSpass-M ArcGIS integrated model as input parameters. The weather/climatic data (potential evapotranspiration, temperature, wind speed, and rainfall), the terrain, the soil texture, the land use/land cover, the leaf-area index data, and the groundwater depth are all basic inputs for this model. It also takes into account a number of other variables [14,28]. ArcGIS was used to prepare all of the input data, which were then imported into the model with the same cell size in ascii format. To prepare the input data from 2000 to 2019, ArcGIS 10.3.1 was used, and the raster cell size was maintained at 100 m × 100 m. Each parameter has a total of 72,713 grids. WetSpass is based upon the principle of the water balance. As a result, the average annual recharge from the rainfall was between 99.10 and 177.2 mm, with a mean value of 119.4 mm per year, as is displayed in (Figure 4a). The long-term annual evapotranspiration values for the years from 2000 to 2019 ranged between 134.4 and 287.2 mm/year, with a mean value of 231 mm/year, as is depicted in (Figure 4b).

Groundwater Flow Model Setup
A fully distributed three-dimensional groundwater flow model, MODFLOW-2005 [41], was developed with the help of ModelMuse [42] as a graphical user interface. The partial differential equation that was used in MODFLOW for the groundwater flow is from [43]:  , with each layer having a thickness of 8, 32, and 67 m, respectively. In this study, the aquifer layers were defined on the basis of the lithology [32]. According to the lithology of the aquifer model, Layer-1 was defined as "unconfined", Layer-2 was defined as "unconfined to convertible", and Layer-3 of the third aquifer was defined as "confined". To make sure, the precision of the model grid was discretized by 100 m by 100 m, subsequent in a purview of 113 rows and 643 columns and having a total of 72,713 cell numbers, and three layers calculated by the WetSpass model. In boundary conditions, the average groundwater level over twenty years was assigned as the initial hydraulic head. The eastern boundary was taken at Chenab River (Stream-1), and it was assigned river packages to simulate the recharge and river stage impacts at the aquifer. The northern boundary, which is represented by the escape canal (Stream-2), was also assigned river packages. All of the other channels in the study area were also assigned river packages. The western and southern boundary was assigned as a general head boundary on the basis of the records of the observation wells, as is shown in ( Figure 5). Well packages (WEL) (used in MODFLOW to simulate the GW flux) were applied to the pumping wells, and recharge packages (RCH) (used to simulate the recharge at the top of the model in MODFLOW) were applied to the groundwater data.

Groundwater Flow Model Setup
A fully distributed three-dimensional groundwater flow model, MODFLOW-2005 [41], was developed with the help of ModelMuse [42] as a graphical user interface. The partial differential equation that was used in MODFLOW for the groundwater flow is from [43]: where x, y, and z are the three-dimensional major axes of the hydraulic conductivity coordinates along the K xx , K yy , and K zz [LT −1 ]; h is the head, which is known to be potentiometric with each layer having a thickness of 8, 32, and 67 m, respectively. In this study, the aquifer layers were defined on the basis of the lithology [32]. According to the lithology of the aquifer model, Layer-1 was defined as "unconfined", Layer-2 was defined as "unconfined to convertible", and Layer-3 of the third aquifer was defined as "confined". To make sure, the precision of the model grid was discretized by 100 m by 100 m, subsequent in a purview of 113 rows and 643 columns and having a total of 72,713 cell numbers, and three layers calculated by the WetSpass model. In boundary conditions, the average groundwater level over twenty years was assigned as the initial hydraulic head. The eastern boundary was taken at Chenab River (Stream-1), and it was assigned river packages to simulate the recharge and river stage impacts at the aquifer. The northern boundary, which is represented by the escape canal (Stream-2), was also assigned river packages. All of the other channels in the study area were also assigned river packages. The western and southern boundary was assigned as a general head boundary on the basis of the records of the observation wells, as is shown in ( Figure 5). Well packages (WEL) (used in MODFLOW to simulate the GW flux) were applied to the pumping wells, and recharge packages (RCH) (used to simulate the recharge at the top of the model in MODFLOW) were applied to the groundwater data. The hydraulic properties were defined on the basis of the available bore log data. These included the horizontal and vertical hydraulic conductivities, as well as the aquifer storage properties. The hydraulic conductivity that was given to the model varied between 1 and 119 m/day as the maximum and minimum values, respectively, with an average value of 66 m/day. A similar hydraulic conductivity range within the province of the Punjab realm was used by the authors of [44][45][46]. Only 5% of Layer-1 had a horizontal hydraulic conductivity (HHC) value in the range of 1-30 m/day, 80% had an HHC value in the range of 30-80 m/day, and the remaining 15% had an HHC value in the range of 80-108 m/day, as is shown in (Figure 6a). From the spatial distributions of the HHCs for Layer-2 and Layer-3, it was found that only 3% of Layer-2 and Layer-3 had an HHC value in the range of 1-30 m/day, 80% had an HHC value in the range of 30-80 m/day, and 17% had an HHC value in the range of 80-119 m/day, as is shown in (Figure 6b).  The hydraulic properties were defined on the basis of the available bore log data. These included the horizontal and vertical hydraulic conductivities, as well as the aquifer storage properties. The hydraulic conductivity that was given to the model varied between 1 and 119 m/day as the maximum and minimum values, respectively, with an average value of 66 m/day. A similar hydraulic conductivity range within the province of the Punjab realm was used by the authors of [44][45][46]. Only 5% of Layer-1 had a horizontal hydraulic conductivity (HHC) value in the range of 1-30 m/day, 80% had an HHC value in the range of 30-80 m/day, and the remaining 15% had an HHC value in the range of 80-108 m/day, as is shown in (Figure 6a). From the spatial distributions of the HHCs for Layer-2 and Layer-3, it was found that only 3% of Layer-2 and Layer-3 had an HHC value in the range of 1-30 m/day, 80% had an HHC value in the range of 30-80 m/day, and 17% had an HHC value in the range of 80-119 m/day, as is shown in (Figure 6b). The hydraulic properties were defined on the basis of the available bore log data. These included the horizontal and vertical hydraulic conductivities, as well as the aquifer storage properties. The hydraulic conductivity that was given to the model varied between 1 and 119 m/day as the maximum and minimum values, respectively, with an average value of 66 m/day. A similar hydraulic conductivity range within the province of the Punjab realm was used by the authors of [44][45][46]. Only 5% of Layer-1 had a horizontal hydraulic conductivity (HHC) value in the range of 1-30 m/day, 80% had an HHC value in the range of 30-80 m/day, and the remaining 15% had an HHC value in the range of 80-108 m/day, as is shown in (Figure 6a). From the spatial distributions of the HHCs for Layer-2 and Layer-3, it was found that only 3% of Layer-2 and Layer-3 had an HHC value in the range of 1-30 m/day, 80% had an HHC value in the range of 30-80 m/day, and 17% had an HHC value in the range of 80-119 m/day, as is shown in (Figure 6b).  The vertical hydraulic conductivity of Layer-1 was 1 ≤ over 75% of the total area, with 5% in a range of 1-1.5 m/day, and the remaining 20% in a range of 1.5-2 m/day. For the distribution of the spatial vertical hydraulic conductivities for Layer-2 and Layer-3, it was found that 45% of the aquifer was in the range of 0-3 m/day, 35% percent was in the range of 3-4 m/day, and 20% was in the range of 5-6 m/day. The specific yield (Sy) ranged between 0.05 and 0.25 as the minimum and maximum values, respectively. The average daily recharge ranged between 0.000271 and 0.000597 m/day as the minimum and maximum values, respectively, with an average value of 0.000387 mm/day.

Coupling of Surface Water Model and Groundwater Flow Model
Because of its lumped nature, the WetSpass model is essentially restricted in terms of dealing with the groundwater flow. MODFLOW, on the other hand, has problems identifying the distributed evapotranspiration and the groundwater recharge, which are the primary inputs to the groundwater model. As a result, it appears that, if WetSpass-based evapotranspiration and groundwater recharge are used as the input data in MODFLOW, and the groundwater level is estimated and transferred to WetSpass, the spatiotemporal characteristics of the study region will be adequately represented. The coupled WetSpass-MODFLOW framework in this study integrated calibrated hydrological modeling (WetSpass-M) with a calibrated groundwater flow model (MODFLOW), as is illustrated in Figure 7. WetSpass simulates the water balance components. MODFLOW, on the other hand, models the threedimensional groundwater flow and its related sources and sinks (e.g., recharge, discharge to drains, interaction with lakes, and interaction with stream networks). The outputs of the WetSpass-M model were utilized as the inputs to MODFLOW. MODFLOW calculates the groundwater hydraulic head and the groundwater-surface water interactions, which are then fed into WetSpass. A coupled model of this type was used to analyze the current stage and the future conditions of the hydrology (water balance components, system water budget, and groundwater level) of the Khadir canal subdivision.   (Figure 8c). For the calibration of the groundwater head, 19 observation wells were used over a period of 12 years to record the groundwater levels. The value of R 2 was found to be 0.99, and it is a significant fit between the simulated and observed heads. Calibrated groundwater model parameters are shown in Table 1. The mean errors ranged from 0.10 to 0.32 m, and the root mean square error was 0.0027 min total. Table 2 shows the statistical analysis. Data from 2015 to 2019 were used for the model validation, as is shown in (Figure 8b).

Calibration, Validation, and Sensitivity Analyses
Average long-term parameters, such as the average value of the precipitation, the PET mean, the mean level of the groundwater, and the mean values of the groundwater recharge, were used to create the steady-state model as a starting step for the transient-state model. A total of 19 average observation-well groundwater levels were used to calibrate the steady-state model. The value, R 2 = 0.98, was found to be a significant fit between the simulated and observed heads. The root mean square error was 0.0022 m, while the mean error was −0.018 m. The mean absolute error was 0.51 m. The steady-state model was  (Figure 8c). For the calibration of the groundwater head, 19 observation wells were used over a period of 12 years to record the groundwater levels. The value of R 2 was found to be 0.99, and it is a significant fit between the simulated and observed heads. Calibrated groundwater model parameters are shown in Table 1. The mean errors ranged from 0.10 to 0.32 m, and the root mean square error was 0.0027 min total. Table 2 shows the statistical analysis. Data from 2015 to 2019 were used for the model validation, as is shown in (Figure 8b).  (Figure 8c). For the calibration of the groundwater head, 19 observation wells were used over a period of 12 years to record the groundwater levels. The value of R 2 was found to be 0.99, and it is a significant fit between the simulated and observed heads. Calibrated groundwater model parameters are shown in Table 1. The mean errors ranged from 0.10 to 0.32 m, and the root mean square error was 0.0027 min total. Table 2 shows the statistical analysis. Data from 2015 to 2019 were used for the model validation, as is shown in (Figure 8b).    A sensitivity analysis was performed by using the computer program known as UCODE (universal sensitivity analysis, calibration, and uncertainty evaluation) [47], and by using ModelMate [48]. The program was used for studying the sensitivity of the hydraulic conductivity (HKZ), the groundwater recharge (RCH), the discharge (Q), the evapotranspiration (EVT), the general head boundary (GHB), and the river conductance (RIVC). The parameters of the discharge (Q), the recharge (RCH), and the hydraulic conductivity (HKZ4) were found to be more sensitive, and the parameters of the hydraulic conductivity (HKZ2) and the general head boundary represent the least sensitive parameters. Moreover, the pumping rate (Q) and the groundwater recharge were found to be highly sensitive groundwater model parameters, as is shown in (Figure 9).  Table 2. Statistical analysis of field and model data.

Calibration Validation
Mean Error conductivity (HKZ2) and the general head boundary represent the least sensitive parameters. Moreover, the pumping rate (Q) and the groundwater recharge were found to be highly sensitive groundwater model parameters, as is shown in (Figure 9).

Groundwater Flow Model Budget
The groundwater budget at the transient-state condition of the calibrated MODFLOW model at the last time step is displayed in Table 3. The results show that the inflow from the river seepage was 1,216,020.25 m 3 d −1 , and that it contributed 54.84% of the total inflow (2,217,328.50 m 3 d −1 ) to the aquifer. The seepage from the general head boundary represented 29.45% of the total inflow to the aquifer. The groundwater recharge contributed 15.71% to the inflow to the aquifer. The outflow quantification of the aquifer showed that 87.28% (1,935,357.59 m 3 d −1 ) of the total outflow was due to the pumping of the wells from the aquifer, and that the remaining 3.35% of the outflow represented the

Groundwater Flow Model Budget
The groundwater budget at the transient-state condition of the calibrated MODFLOW model at the last time step is displayed in Table 3. The results show that the inflow from the river seepage was 1,216,020.25 m 3 d −1 , and that it contributed 54.84% of the total inflow (2,217,328.50 m 3 d −1 ) to the aquifer. The seepage from the general head boundary represented 29.45% of the total inflow to the aquifer. The groundwater recharge contributed 15.71% to the inflow to the aquifer. The outflow quantification of the aquifer showed that 87.28% (1,935,357.59 m 3 d −1 ) of the total outflow was due to the pumping of the wells from the aquifer, and that the remaining 3.35% of the outflow represented the general head boundary (GBH), with 6.96% due to river recharge, and 2.24% due to evapotranspiration.

Scenarios
The calibrated model was run for different scenarios of the investigated area. These scenarios were developed on the basis of possible situations that may cause a shortage of water resources that may stress the aquifer in the coming days. The developing scenarios are helpful for the development of proper measures and management strategies to counter future problems.
Scenario-I was simulated to assess the impact of pumping if the prevailing conditions of the years from 2003 to 2019 were to continue until 2035, as there are no regulations to control the extensive tube well installation in that area. The scenario was further subdivided into I-A, I-B, I-C, and I-D, according to the simulation times of 2020, 2025, 2030, and 2035, respectively. In this scenario trend of increasing pumping, the volume was also analyzed, as a massive increase in tube wells in the study area was observed. Scenario-II was used to simulate the impacts of pumping on the aquifer by increasing the pumping capacity by 25, 50, 75, and 100% for the coming 10 years. This means that, in this scenario, the aquifer status was analyzed on the basis of increasing the pumping rate with the aforementioned percentages over a period of ten years. It was further subdivided into the scenarios, II-A, II-B, II-C, and II-D, according to the increase in the pumping capacity. Scenario-III was subdivided into the subscenarios, III-A and III-B, under the conditions that were simulated by decreasing, on average, the groundwater recharge from the rainfall by 50%, and from the river by 50%, following the same pumping that is shown in Table 4.
In this study, three main scenarios (which were divided into ten subscenarios) were simulated for the future prediction of the groundwater. A summary of these scenarios is depicted in Table 4.
The results of Scenario I-A show that there was a 5.0 m decrease in the groundwater table in 2020 from the reference point of the calibrated steady-state model (2000-2019), with a pumping rate of 2,915,029 m 3 /day. In the 1-B scenario, the model predicted that the groundwater table would decrease by up to 9.8 m in the year 2025, which is due to the pumping rate of 3,441,109 m 3 /day from the level. For Scenario I-C, the groundwater was predicted to decline by 14.6 m from the reference because of the increase in the pumping rate to 3,904,630 m 3 /day in 2030, and Scenario I-D simulated a decrease of 18.1 m in the year 2035, and the pumping rate was predicted to be 4,262,693 m 3 /day in 2035, as is shown in Table 4. The simulated results of Scenario II-A show that, by increasing the discharge from pumping by 25%, the groundwater table would decrease by up to 2.0 m by the end of 2029, from the reference. The Scenario II-B results show a decrease of 5.5 m at 50% increased pumping. The Scenario II-C results show a decrease of 9.8 m at 75% increased pumping, and the Scenario II-D results show a decrease of 14.3 m at 100% increased pumping in the year 2029, from the reference. The Scenario III-A results predicted that, by decreasing the recharge from the river and other water bodies by 50%, because of uncertain climate change, drought, etc., the aquifer recharge from the river was 1,069,024 m 3 /day instead of 2,134,048 m 3 /day in 2019, which caused a 2.4 m decrease in the water table from the reference. The Scenario III-B results predicted that, by decreasing the recharge from the rainfall by 50%, because of uncertain climate change, drought, etc., the aquifer recharge from the rainfall was 182,482 m 3 /day instead of 364,890 m 3 /day in 2019, which caused a 0.7 m decrease in the water table from the reference. A graphical representation of the recharge to, and the discharge from, the aquifer under different scenarios is shown in (Figure 10).  (Figure 10). Because of the imbalance between the groundwater abstraction and the recharge at different locations, many mixed drawdown levels were seen at every observation well. It was simulated that, at most of the locations, there was lowering in the water table, and at only very few observation wells, there was lowering in the groundwater level. The model predicted that the groundwater level in the vicinity of the monitoring in Pir Panja, and towards the riverside of the Khadir distributary, would be little affected by the discharge until 2035. Similarly predicted results show that the groundwater levels in the observation wells in the vicinity of Vijhalke Minor, Kalri Minor, and Walla Minor were decreasing less than the observation wells that were installed near Khadir Minor. For Scenario-I, the contour intervals of the groundwater table were kept constant for each year. The ranges of the minimum and maximum contour lines of the groundwater table varied from 159 to 189 m, and from 138 to 186 m, for the years 2020 and 2035, respectively, as is shown in Figure 11a,b, respectively. Because of the imbalance between the groundwater abstraction and the recharge at different locations, many mixed drawdown levels were seen at every observation well. It was simulated that, at most of the locations, there was lowering in the water table, and at only very few observation wells, there was lowering in the groundwater level. The model predicted that the groundwater level in the vicinity of the monitoring in Pir Panja, and towards the riverside of the Khadir distributary, would be little affected by the discharge until 2035. Similarly predicted results show that the groundwater levels in the observation wells in the vicinity of Vijhalke Minor, Kalri Minor, and Walla Minor were decreasing less than the observation wells that were installed near Khadir Minor. For Scenario-I, the contour intervals of the groundwater table were kept constant for each year. The ranges of the minimum and maximum contour lines of the groundwater table varied from 159 to 189 m, and from 138 to 186 m, for the years 2020 and 2035, respectively, as is shown in Figure 11a,b, respectively. In the higher side of the study area, a substantial drawdown was observed, while, on the lower side, or riverside, of the study area, there was less of a decrease in the water table until 2035. The greater decrease in the upper side was due to less recharge from the river, as it lies away from the river. This lowering in the groundwater level in the future may cause a very drastic situation for crop production as well as for domestic and industrial uses. Decreases in the groundwater level may have drastic impacts on pumping, such as an increase in the pumping cost, and a decrease in the pumping efficiency. The results of a study that was based on the groundwater economy, which was conducted by the authors of [4] in the Indus basin, Pakistan, show that the decline in the groundwater table could have significant impacts on the pumping cost. A shallow electronic turbine (operating < 6 m) could cost up to USD 1000, while a deep electronic turbine/tube well (operating > 20 m) could cost up to USD 5000 per annum. Similarly, the pumping cost of deep tube wells was also more (i.e., USD 8.61 per acre inch, compared to USD 18.78 for a 31 and 91 m lift in the water [49]). In the higher side of the study area, a substantial drawdown was observed, while, on the lower side, or riverside, of the study area, there was less of a decrease in the water table until 2035. The greater decrease in the upper side was due to less recharge from the river, as it lies away from the river. This lowering in the groundwater level in the future may cause a very drastic situation for crop production as well as for domestic and industrial uses. Decreases in the groundwater level may have drastic impacts on pumping, such as an increase in the pumping cost, and a decrease in the pumping efficiency. The results of a study that was based on the groundwater economy, which was conducted by the authors of [4] in the Indus basin, Pakistan, show that the decline in the groundwater table could have significant impacts on the pumping cost. A shallow electronic turbine (operating < 6 m) could cost up to USD 1000, while a deep electronic turbine/tube well (operating > 20 m) could cost up to USD 5000 per annum. Similarly, the pumping cost of deep tube wells was also more (i.e., USD 8.61 per acre inch, compared to USD 18.78 for a 31 and 91 m lift in the water [49]). In the higher side of the study area, a substantial drawdown was observed, while, on the lower side, or riverside, of the study area, there was less of a decrease in the water table until 2035. The greater decrease in the upper side was due to less recharge from the river, as it lies away from the river. This lowering in the groundwater level in the future may cause a very drastic situation for crop production as well as for domestic and industrial uses. Decreases in the groundwater level may have drastic impacts on pumping, such as an increase in the pumping cost, and a decrease in the pumping efficiency. The results of a study that was based on the groundwater economy, which was conducted by the authors of [4] in the Indus basin, Pakistan, show that the decline in the groundwater table could have significant impacts on the pumping cost. A shallow electronic turbine (operating < 6 m) could cost up to USD 1000, while a deep electronic turbine/tube well (operating > 20 m) could cost up to USD 5000 per annum. Similarly, the pumping cost of deep tube wells was also more (i.e., USD 8.61 per acre inch, compared to USD 18.78 for a 31 and 91 m lift in the water [49]).
For sustainable management, and for aquifer replenishment, the groundwater recharge from rainfall must be increased through the promotion of rainwater harvesting in order to assess the continuous decrease of using it for groundwater recharge. In areas that are suffering a continuous decrease in the groundwater table, the techniques of artificial recharge must be applied by using the rainfall water and extra surface irrigation water for the groundwater replenishment. There should be the formation and implementation of legal boundaries to evaluate the unjudicial abstraction, and the unnecessary use of the groundwater. This study may be extended to the regional scale by incorporating further management and climate change scenarios.

Model Uncertainty
The current model produces a good interaction between the surface water and the groundwater for the proper management of groundwater resources. The data required a comprehensive quantitative analysis for their uncertainty. Most of the major uncertainties in the model are due to the available data constraints, such as the measured groundwater levels and observation wells, as well as poor meteorological, soil, physical, and hydraulic data availability. Recording the river's recharge over time and a comparison of the historical data with the forecasts for current and future conditions, as well as the refinement of the groundwater recharge, would help the model to become more accurate in terms of the short-term future uncertainty. Using the current conceptual model as a starting point, it will be possible to refine the model and reduce its uncertainty as the network for monitoring the groundwater levels increases.

Conclusions
The management of groundwater supplies is critical to the achievement of agricultural, domestic, and industrial growth. The future potential of the Chaj doab aquifer in Punjab, Pakistan, was studied by using a distributive surface water model, which is known as WetSpass, and a groundwater flow model, which is known as MODFLOW. With twelve years of data (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014), and five years of validation data, the groundwater flow model, MODFLOW, was calibrated and validated manually (2015-2019). The R 2 was found to be 0.99, which is an acceptable fit between the simulated and observed heads. In this work, three primary scenarios (which were separated into ten subscenarios) were simulated for the future groundwater predictions, and the root mean square error was 0.0027 m. The predicted results of Scenario-I indicate that the groundwater in the study area may decrease from 1.5 to 18.1 m over a period of 16 years (2019-2035). A greater decrease was observed in the area where the tube well density was high. Similarly, a smaller decrease in the groundwater table (from normal to medium) was simulated in the lower part of the study area of concern.
In Scenario-II, the results were simulated by increasing the pumping rate by 25, 50, 75, and 100% over 10 years (2020-2029). The results predict that, after 10 years, the groundwater table would decrease by up to 2.0, 5.3, 9.8, and 14.3 m through the increases in the pumping rate of 25, 50, 75, and 100%, respectively. The Scenario-III results reveal that, with the decreases in the groundwater recharge from the rainfall and the river by 50% due to climate change or to a decrease in the river recharge, the mean decreases in the water table would be 0.7 m and 2.4 m, respectively. A greater decrease in the water table through the decrease in the groundwater table was recorded by the river recharge, compared to the recharge from the rainfall. On the higher side of the study area, a substantial drawdown was observed, while, on the lower side, or the riverside, of the study area, there was less of a decrease in the water table until 2035. The greater decrease in the upper side was due to less recharge from the river, as it lies away from the river. Spatially, greater decreases in the water table were observed in the areas away from the river. This study would be helpful for the development of future strategies for the proper management of aquifers.