Impact of Crop Diversification on Groundwater Sustainability in the Barind Area of Bangladesh: A Case Study on Replacing Boro Rice with Wheat Using Mathematical Modelling ()
1. Introduction
Groundwater in Bangladesh transpires at a very shallow depth, where the recent river-borne sediments form prolific aquifers in the floodplains. The Pliocene Tipam sands in the hilly areas serve as aquifers. In the higher terraces, the Barind and Madhupur tracts, the Pleistocene DupiTila sands act as aquifers [1]. The groundwater level is at or very close to the surface during the monsoon, whereas it is at maximum depth during the months of April and May. This trend is common over most of Bangladesh except Dhaka City and the Barind Tract [2].
Barind Tract, the largest Pleistocene physiographic unit of the Bengal basin, covering an area of about 7727 sq km [1] [2], can be divided into high, medium and low based on their elevation [3]. Elevation of the area varies from 9 m to 47 m PWD (Public Works Datum) [4]. Because of the elevation of high Barind, Rajshahi is one of the most drought prone districts of Bangladesh [5]. The impact of drought can be much higher and can cause greater loss than flood, cyclone and storm surge [6]-[8]. Drought is related to groundwater recharge.
Groundwater recharge in Bangladesh mainly occurs through monsoon rainfall and flooding. Due to its higher elevation (topography varies from 20.0 m PWD to 47.0 m PWD) [3], most of the areas of Barind area are flood free zone. So, the main source of groundwater recharge in this area is rainfall [9]. Apart from Rajshahi, where the annual rainfall is about 1600 mm, most parts of the country receive at least 2000 mm of annual rainfall [10]. Moreover, the thick sticky clay surface of Barind Tract acts as aquitard which impedes groundwater recharge and increases surface run-off [11]. As a result, groundwater level in this part is successively falling by years with increasing withdrawal of water for irrigation [11]. The irrigation in this area is mainly dependent on groundwater due to absence or presence of a very limited amount of surface water. As a result of excess abstraction lower than that of recharge, groundwater level is declining at an alarming rate in this area. Over abstraction of groundwater, lack of surface water bodies, low rainfall, high elevation, thick clay layer, etc., are the major hindrances in the study area to sustain groundwater resources.
A recent study shows that groundwater levels in some areas falls between 5 - 10 m in the dry season, and most of the tube wells fail to lift sufficient water [12]. The GW dependent irrigation scheme in the area has reached a critical stage where GW level has dropped under the depth of the shallow tube wells in several places [13]. Rice governs the cropping pattern of Barind, which suffer from scarcity of groundwater in the dry season. Only Aman paddy in wet season was cultivated in Barind [14]. Due to the expansion of GW irrigation after the 1980s, High Yielding Variety (HYV) paddies are familiarized in this area. Now Barind Tract produces three crops in one agricultural season with the blessing of groundwater irrigation [11]. Researchers and policymakers are advocating sustainable development as the best approach to today’s and future water problems [15] [16]. But sustainability of groundwater resources is at risk in terms of quantity in the northwest region [17]. So, it is urgent to take necessary initiatives to reduce the declining trend of groundwater and its sustainability in this area.
The primary objective of this research is to advance the sustainable management of groundwater resources in Bangladesh’s most drought-prone region. Specifically, this study addresses the decision question: Which future groundwater scenario represents the most critical threat to sustainability, and can crop substitution replace high water-demand Boro rice with more water-efficient wheat serve as an effective intervention to mitigate groundwater decline? To answer this, we employ a coupled MIKE SHE-MIKE 11 modeling framework to (i) screen ten future scenarios to identify the most extreme case, and (ii) quantify the hydrogeological impact of replacing Boro rice with wheat as a demand-side management strategy.
2. Description of the Study Area
For research purposes, Rajshahi District has been chosen as the study area which is around 2407.01 sq km and is located in between 24˚07' and 24˚43' north latitudes and in between 88˚17' and 88˚58' east longitudes. It is bounded by Naogaon district on the north, West Bengal state of India, Kushtia district and Ganges River on the south, Natore district on the east, Nawabganj on the west. The study area is shown in Figure 1.
Figure 1. Study area (Rajshahi district).
The topography of the study area varies from 9.72 mPWD to 46.85 mPWD. Almost 13% of study areas have elevation greater than 25.00 mPWD. Most of the area remains flood free because of its high elevation. Several Upazilas of this district were identified as resource constraint area [3].
2.1. Climatological Conditions
The average annual temperature in Rajshahi is 25.8˚C. The variation in temperatures throughout the year is 10.9˚C. The annual average precipitation is around 1419 mm. Precipitation is the lowest in December, with an average of 2 mm and highest in July, with an average of 301 mm. The histogram of average monthly rainfall of the study area (Rajshahi District) is shown in Figure 2.
Figure 2. Average monthly rainfall in Rajshahi district.
2.2. River System
There are ten rivers in this district, totaling around 146 km in length. Major rivers of the study area are Ganges, Mahananda, Baral and Sib-Barnai. Moreover, there are several minor rivers in this area. Most of the rivers of this region flow from very steep to flat ground. The river systems are included in the model described earlier. The river model is coupled with the groundwater model. The study area hasn’t enough major rivers but there are some major rivers surrounding its area. Some of the name of rivers/khals are given below in a list and shown in Figure 3 and listed in Table 1.
3. Methodology
The methodology of the study includes selection of the study area, data collection and processing, development of surface water (SW) model, development of groundwater (GW) model, development of SW-GW interaction model, calibration and validation of the model, selection of design year, formulation of different scenarios, simulation of the developed scenarios, result analysis of different scenarios and concluding remarks. The schematic diagram of the methodologies of the study is shown in Figure 4.
Table 1. Name of rivers in the study area.
S.N. |
River Name |
S.N. |
River Name |
1 |
Atrai |
15 |
Little Jamuna |
2 |
Atrai-Sib Link |
16 |
Little Nagor-Atrai Link |
3 |
Baral |
17 |
Little Nagor |
4 |
Baranai Branch |
18 |
Mohanada |
5 |
Borokhari |
19 |
Monohorganga Khal |
6 |
Durgapur Khal |
20 |
Musakhan |
7 |
Fakirni |
21 |
Nandakuja |
8 |
Ganges |
22 |
Narod |
9 |
Ghaelgurhari Khal |
23 |
Pagla |
10 |
Godai |
24 |
Rasulpur Khal |
11 |
Hoja Khal |
25 |
Sabaibeel |
12 |
Joai Khari |
26 |
Sib-Barnai |
13 |
Kamargaon Khari |
27 |
Zia Khal |
14 |
Khalishadingi |
28 |
Pangal |
Figure 3. River network of the study area and surroundings.
Figure 4. Schematic diagram of methodology of the study.
3.1. Selection of the Study Area
Rajshahi district is one of the most drought prone districts of Bangladesh [5]. Groundwater level in this area is successively falling at an alarming rate [18]. Groundwater recharge conditions are very poor in Tanore, Godagari, Mohanpur and Baghmara upazilas and vulnerable for Boro rice [18]. Groundwater dependent irrigation system in the area has reached a critical phase [13] and GW level has dropped below the depth of the shallow tube wells in many places [13]. The government of Bangladesh is about to announce a “state of water emergency” in the drought prone Barind tract in the northwest of the country, where over-extraction of groundwater for rice farming and dwindling rainfall caused by climate change have combined to create a crisis (thethirdpole.net). [3] has also identified several upazilas of Rajshahi district as resource constraint areas. This is why this district has been chosen as study area where sustainable GW is a major concern.
3.2. Data Collection and Data Processing
According to the model requirements, significant amount of data have been collected from different relevant secondary sources, mainly from Institute of Water Modelling (IWM), Bangladesh Water Development Board (BWDB), Bangladesh Agricultural Development Corporation (BADC), Barind Multipurpose Development Authority (BMDA) and Soil Resource Development Institute (SRDI). Historical time series data for river water level and groundwater level available from BWDB stations have been considered for developing surface water model using MIKE 11 and groundwater model using MIKE SHE respectively, in current study. All the data collected has to be used in this study after checking its quality & consistency and then processed as per required format for the model setup. In addition to data quality checking, data analysis has also been carried out for estimation of different model parameters.
3.3. Base Model Set Up
The hydrodynamic model of the project area has been set up using Hydrodynamic Module (HD) of the MIKE 11 tool. The study area is approximately 2400 square km and the model area has been considered larger than this for minimizing the boundary effect and to use observed water level and discharge stations of BWDB and IWM. The base model has been set up for 5 years (January 2012 to December 2016). Expected major outputs from this model are water level and discharge at every grid point within model area.
The groundwater model has been set up using MIKE SHE model of DHI. This model has been set up for 5 years (January 2012 to December 2016). Expected major outputs from this model are groundwater levels at every grid point within model area and water balance of the model area. The MIKE 11(HD) and MIKE SHE models have been linked interactively and then they can produce water balance and change in storage in the form of groundwater recharge/discharge and show fluctuations in the water table.
Model Domain and Discretization
Spatial Discretization:
Model domain: ~2900 km2 (extended beyond the 2407 km2 study area to minimize boundary effects).
Grid resolution: 500 m × 500 m regular grid (approximately 11,600 active cells for the groundwater model).
Vertical discretization: 3 layers representing i) upper clayey aquitard (0 - 10 m), ii) main Dupi Tila aquifer (10 - 80 m), and iii) lower semi-confined layer (80 - 120 m).
Hydrogeological Configuration:
Aquifer type: Unconfined to semi-confined Pleistocene Dupi Tila sands.
Horizontal hydraulic conductivity (Kₓ): Calibrated range 1.0 × 10−5 to 8.5 × 10−5 m/s.
Vertical hydraulic conductivity (Kγ): 1.0 × 10−6 to 5.0 × 10−6 m/s.
Specific yield (Sγ): 0.15 - 0.25 (upper layer), 0.05–0.15 (deeper layers).
Storage coefficient: 1.0 × 10−4 to 5.0 × 10−4.
Boundary Conditions:
Lateral boundaries: Eight boundary wells with time-varying specified head, interpolated from observed hydrographs.
River boundaries: 28 river reaches (146 km total length) including Ganges, Mahananda, Baral, and Sib-Barnai, modeled as head-dependent flux boundaries using MIKE 11 HD coupled with MIKE SHE.
Upper boundary: Specified recharge and evapotranspiration; irrigation applied as negative recharge (pumping).
Lower boundary: No-flow boundary at 120 m depth (impermeable basement).
Observation Network:
Groundwater: 24 monitoring wells (BWDB and IWM networks) distributed across all nine Upazilas, with continuous water level records (2012-2016).
Surface water: 6 gauging stations (Mohananda at Chapainawabganj, Atrai at Bagmara, Baral at Charghat, plus 3 minor stations).
Irrigation inputs: 15,847 DTW (Deep Tube Well) and STW (Shallow Tube Well) locations from BADC database, with monthly abstraction rates.
Model Coupling:
MIKE 11 (1D hydrodynamic) and MIKE SHE (3D distributed) dynamically coupled at 1-day time steps.
River-aquifer exchange calculated using Darcy’s law with streambed conductance: 1.0 × 10−4 to 5.0 × 10−4 m/s.
3.4. Calibration and Validation of Base Models
The surface water and groundwater models have been calibrated individually for three years (January 2012 to December 2014) and have been validated for the next 2 years (January, 2015 to December, 2016). Then the coupled model (MIKE 11 and MIKE SHE) has also been calibrated and validated with observed groundwater level data. Main calibration parameter for MIKE-11(HD) model is Bed Resistance of stream channel (Manning’s n) and for MIKE SHE model, calibration parameters are horizontal hydraulic conductivity, vertical hydraulic conductivity, specific yield and storage co-efficient etc.
Calibration and Validation of SW Model: The MIKE 11 model has been calibrated and validated with observed surface water level. This model has been calibrated for 3 years (January to December of 2012 to 2014) and has been validated for next 2 years (January to December of 2015 to 2016). Main calibration parameter for MIKE-11(HD) model is Bed Resistance of stream channel (Manning’s n). The comparison plot of model simulated and observed data for surface water level on Mohananda and Atrai river is shown in Figure 5 and Figure 6 respectively.
From Figure 5 and Figure 6, it is observed that there is a quite good matching of model simulated water level with the observed water level on Mohananda River, both in dry as well as in wet season, though there is a slight gap during the monsoon or wet season, in the case of Atri. So, it may be considered that the model is well calibrated.
Calibration and Validation of SW-GW Integrated Model: SW model has been coupled with GW model and this integrated model has been calibrated for 3 years (January to December of 2012 to 2014) and has been validated for next 2 years (January to December of 2015 to 2016). The comparison of model simulated and observed data of groundwater level at Puthia and Durgapur is shown in Figure 7 and Figure 8 respectively. From the Figure, it is observed that the matching is quite satisfactory.
Figure 5. Comparison observed & model simulated water level on the Mohananda River.
Figure 6. Comparison observed & model simulated water level on Atri River.
Figure 7. Comparison of observed and model simulated groundwater level at Puthia Upazila.
Figure 8. Comparison of observed and model simulated groundwater level at Durgapur Upazila.
Model Performance Metrics
Model performance metrics are tabulated in Table 2.
Table 2. Model performance metrics.
Model Component |
Period |
Metric |
Value |
Target/Interpretation |
Surface Water (MIKE 11) |
Mohananda River |
Calibration (2012-2014) |
NSE |
0.84 |
>0.75 (good) |
RMSE (m) |
0.42 |
<0.5 m (excellent) |
Bias (m) |
−0.08 |
Near zero |
Validation (2015-2016) |
NSE |
0.81 |
>0.75 (good) |
RMSE (m) |
0.51 |
<0.6 m (acceptable) |
Atrai River |
Calibration (2012-2014) |
NSE |
0.79 |
>0.75 (good) |
RMSE (m) |
0.58 |
<0.6 m (acceptable) |
Bias (m) |
0.12 |
Slight overestimation |
Validation (2015-2016) |
NSE |
0.76 |
>0.75 (good) |
RMSE (m) |
0.63 |
<0.7 m (acceptable) |
Groundwater (MIKE SHE) |
Puthia Upazila |
Calibration (2012-2014) |
NSE |
0.71 |
>0.65 (satisfactory) |
RMSE (m) |
1.84 |
<2.0 m (good for regional model) |
MAE (m) |
1.42 |
<1.5 m |
Bias (m) |
0.35 |
Slight overprediction |
Validation (2015-2016) |
NSE |
0.68 |
>0.65 (satisfactory) |
RMSE (m) |
2.12 |
<2.5 m (acceptable) |
Durgapur Upazila |
Calibration (2012-2014) |
NSE |
0.74 |
>0.65 (satisfactory) |
RMSE (m) |
1.67 |
<2.0 m |
MAE (m) |
1.28 |
<1.5 m |
Bias (m) |
−0.22 |
Slight underprediction |
Validation (2015-2016) |
NSE |
0.69 |
>0.65 (satisfactory) |
RMSE (m) |
1.95 |
<2.5 m |
3.5. Selection of Design Year
For the development of future scenarios, design year has been selected based on statistical analysis. In the present study, design year has been selected based on return period of mean annual rainfall of the study area. The mean annual rainfall has been obtained from the average of 10 stations falls in the study area. Observed annual rainfall for a period of 47 years (1970-2016) has been considered for statistical analysis. According to the recommendation of FAP25 study, data has been fitted to 3-parameter Log Normal distribution to find out the average and extreme dry year. The statistical software HYMOS 4.0 has been used for this purpose. From this analysis 1990 year has been selected as a design year of average hydrological condition and 2014 and 1994 have been selected for design year of dry and extreme dry conations respectively.
3.6. Development Future Scenarios
In order to sustain groundwater resources up to year 2030 we have to foresee future condition of groundwater resources under different Scenarios. For this reason, there are ten (10) Scenario s have been chosen to understand future groundwater level. It also helps to find out the most extreme scenario. Scenario 10 in which rainfall and evaporation data generated from GCM, boundary GWL of extreme dry event of 10 years return period, population demand from BBS and crop water demand of the base year 2016 has been increased annually by 2.5% up to year 2030 has been found the most extreme scenario from analysis.
Option for replacement of Boro with Wheat: Production of rice requires huge amount of groundwater for irrigation purposes during non-monsoon period (December to May), which lead to progressive lowering of groundwater tables. Around 90% of the irrigation water comes from groundwater sources in Bangladesh. To achieve food security, Bangladesh has increased Boro rice production substantially at the cost of unsustainable groundwater depletion that will hinder water and food security in the long run. So, in this option, the area of Boro rice field has been replaced by wheat which requires around half of water required by rice.
3.7. Prediction of Data up to Year 2030
For simulating the model for this future scenario all the required predicted data have been used. The major data that has been predicted for use in the model is described below:
Precipitation Data: Monthly precipitation data for 2017 to 2030 have been collected from Global Circulation Model HADCM3 using Climate Editor Tool of Mike 11 software. Emission scenario has been chosen for this analysis, SRA2, which assumes that a very heterogenous world with continuously increasing global population and regionally oriented economic growth that is more fragmented and slower than in another storyline.
Evapotranspiration Data: Monthly evapotranspiration data for 2017 to 2030 have been collected from Global Circulation Model HADCM3 using Climate Editor Tool of Mike 11 software. Emission scenario has been chosen for this analysis, SRA2, which assumes that a very heterogenous world with continuously increasing global population and regionally oriented economic growth that is more fragmented and slower than in another storyline.
River Water Level Data: For projection of groundwater levels up to 2030, it is considered in this study that river water level data follows the same time series plotted for year 2016 that has to be extended up to 2030.
Groundwater Level at Model Boundary: Projected groundwater level data for eight (8) boundary wells are very important variables of the MIKE SHE hydrologic model. For projected groundwater level of the model boundary wells, four (4) boundary conditions named as “Boundary Condition 01”, “Boundary Condition 02”, “Boundary Condition 03” and “Boundary Condition 04” have been considered.
Boundary Condition 01: Assumption is, same groundwater level of hydrological condition of the year 2016 (base condition) individual model boundary wells will be continued over the years up to 2030.
Boundary Condition 02: Assumption is, same groundwater level of hydrological condition of the year 1990 (average condition of 2.33 years return period) individual model boundary wells will be continued over the years up to 2030.
Boundary Condition 03: Assumption is, same groundwater level of hydrological condition of the year 2014 (dry condition of 5 years return period) individual model boundary wells will be continued over the years up to 2030.
Boundary Condition 04: Assumption is, same groundwater level of hydrological condition of the year 1994 (extreme dry condition of 10 years return period) individual model boundary wells will be continued over the years up to 2030.
Irrigation Water Demand: Upazila wise crop water requirement for different cropping patterns has been received from the study report conducted by IWM survey [3]. These data have used in this study for base irrigation demand in the year 2016. According to National Water Management Plan (NWMP), irrigation demands are expected to increase potentially by at least a quarter over the next 25 years [19]. From 2016 to 2030, 3 sets of crop demand increment were taken that exponentially increased from 1.0 in 2016 to values of i) no increment; ii) 1% yearly increment of present demand; iii) 1.5% yearly increment demand; iii) 2% yearly increment demand and v) 2.5% yearly increment of present demand respectively in 2030 is shown in Figure 9. 2016 is considered as the base year (present demand year). To attain 1%, 1.5%, 2% and 2.5% yearly incremental demand, demand for the base year has been multiplied by 1.010, 1.015, 1.020 and 1.025 respectively.
Domestic and Municipal Water Demand: According to the NWMP report, per capita gross water demand for municipal town and rural areas are 166 lpcd and 30 lpcd respectively [19]. The gross water demand of municipal town includes 119 lpcd net domestic water demand, 20% of it as a system loss, 10% as gross commercial demand and 15% as industrial demand. On the other hand, it has 50% returned flow from commercial demand and 75% return flow from industrial demand, thus the net water demand for municipal town becomes 76 lpcd. The gross water demand for rural areas doesn’t include any loss and commercial and industrial demand. Thus, the net water demand for rural areas is same as the gross water demand. In Bangladesh, the population census is carried out in ten years interval. Last population was enumerated in the year 2011 and published by the BBS (2012) which is the only sources of population data. In this study, a government report titled “Population Projection of Bangladesh, Dynamics and Trends (2011-2061)” by BBS under Ministry of Planning has been taken as an authentic source for projection of population of Rajshahi district up to 2030.
![]()
Figure 9. Irrigation demands considered in the models (for 2017-2030).
Domestic and municipal water requirements for future conditions have been assessed based on the population projected to the year 2030. From calculation it has found that average domestic and municipal water requirement from 2012 to 2030 is 0.05 mm/day for all Upazilas except Rajshahi City Corporation where demand is 0.85 mm/day. For simplicity and for low demand compared to crop demand, a constant demand of 0.05 mm/day and 0.85 mm/day have been applied in groundwater model for all Upazilas and Rajshahi city corporation respectively to address population and municipal demand.
Boundary Condition Assumptions and Limitations
Justification for Fixed Boundary Conditions:
1) River Water Levels: The assumption that river water levels follow 2016 patterns through 2030 is based on i) the absence of reliable future river flow projections for the Ganges and tributaries under climate change, and ii) the dominance of upstream water withdrawal and Farakka Barrage operations, which are policy-driven rather than climate-driven. The 2016 base year represents average recent hydrological conditions (neither extreme wet nor dry).
2) Boundary Groundwater Levels: The four boundary conditions (BC-01 to BC-04) represent a gradient of stress scenarios from base (2016) to extreme dry (1994, 10-year return period). This approach allows systematic testing of boundary uncertainty while maintaining computational tractability. Fixed-head boundaries are appropriate because the model domain extends 10 - 15 km beyond the study area, placing boundaries in regions with less intensive pumping stress.
Sensitivity Analysis: To assess the impact of boundary assumptions on scenario ranking, we conducted sensitivity tests varying boundary heads by ±10% from the base case. Results indicate:
Scenario ranking (identification of S-10 as most extreme) remained unchanged across all boundary perturbations.
Absolute groundwater levels shifted by 0.3 - 0.8 m, but the relative difference between scenarios (2.5 - 4.2 m) was preserved.
The Boro-to-wheat intervention effect (improvement magnitude) showed <5% variation across boundary conditions.
Limitations: We acknowledge that fixed historical boundaries do not account for potential future changes in regional groundwater flow patterns due to i) climate change impacts on upgradient recharge, ii) transboundary groundwater flow from India, or iii) large-scale pumping trends in adjacent districts. These limitations suggest our absolute groundwater level projections should be interpreted as scenario-based rather than predictive, with emphasis on relative differences between scenarios and interventions rather than absolute future states.
3.8. Determination of the Most Extreme Scenario from Different Scenarios
Due to randomness of rainfall distribution and cropping pattern, extreme scenario has been selected by analyzing data of the Upazila instead of district. Yearly minimum and maximum groundwater levels have been analyzed for all Upazilas and most emphasis has been given to identify the worst and less vulnerable on minimum groundwater level. At first Upazila wise analysis has been done for all scenarios to find out which scenario is most extreme, and which is good among these ten scenarios. It is expected that extreme and less vulnerable of one Upazila may not be extreme or good for every Upazila of Rajshahi. Therefore, extreme scenario for the maximum number of Upazilas has been selected for the whole Rajshahi district.
4. Result and Analysis
After generation of possible future scenarios, extreme scenario has been identified. So, in this study, interventions have been applied to this extreme scenario to sustain future groundwater resources up to 2030. Attempts have been made to stop the rate of declining groundwater levels and to increase the trends of groundwater levels of the study area.
4.1. Impact on Groundwater Level for Replacement of Boro with Wheat
Before application of intervention, it is necessary to explore the extent of Boro rice field. Total area of Boro rice field within Rajshahi district is approximately 724 km2, which is about one-third of the study area (Figure 10). From Table 3, it is clear that Tanore, Godagari, Baghmara and Mohanpur are the Boro dominant Upazilas of Rajshahi district.
Figure 10. Extent of Boro rice field within rajshahi district.
Table 3. Upazila-wise distribution of Boro rice field.
Upazila Name |
Upazila Area (km2) |
Boro Rice Field (km2) |
% Coverage of Each Upazila |
Tanore |
298 |
183 |
61 |
Godagari |
494 |
223 |
45 |
Baghmara |
368 |
137 |
37 |
Mohanpur |
163 |
57 |
35 |
Paba |
303 |
57 |
19 |
Durgapur |
200 |
37 |
18 |
Puthia |
193 |
18 |
9 |
Bagha |
194 |
11 |
6 |
Charghat |
175 |
2 |
1 |
Replacement of Boro rice by Wheat has significantly improved the water stressed situation within this study area. Almost all the wells of the Upazilas have been responded positively for this intervention. Wise use of groundwater in the dry season could be a better solution towards sustainability of this scarce resource. Graphical representation of extreme scenario vs intervention 01 at Godagari, Paba, Baghmara and Mohanpur Upazila is shown in Figures 11-14. In these graphs, extreme scenario is represented by red line and intervention 01 is represented by green line. Groundwater level simulated under extreme scenario shows significant recovery due to application of intervention.
Figure 11. GWL of extreme scenario vs. intervention 01 at Godagari Upazila.
Figure 12. GWL of extreme scenario vs. intervention 01 at Paba Upazila.
Figure 13. GWL of extreme scenario vs. intervention 01 at Baghmara Upazila.
Figure 14. GWL of extreme scenario vs. intervention 01 at Mohanpur Upazila.
4.2. Depth of Phreatic Surface under Different Scenarios and Interventions
The phreatic surface or water table is the surface where the water pressure head is equal to the atmospheric pressure. It may be visualized as the “surface” of the subsurface materials that are saturated with groundwater in a given vicinity. Depth of this surface is measured from existing ground level. From the statistical analysis, it has been noticed that the 27 April 2028 will be the driest event for this whole study area. That’s why, depth to phreatic surface of different Scenarios and interventions have been shown below for this particular date to visualize specially the improvement of groundwater table after application of different interventions. It helps to understand the future availability of GW due to lowering of GWL. To counter extreme scenario S-10, the intervention has been applied. In this intervention, high water consuming Boro rice has been replaced by low water consuming Wheat (Intervention 01). Due to this intervention, a significant improvement in groundwater level is observed and is shown in Figure 15 and Figure 16.
![]()
Figure 15. Depth of phreatic surface of extreme scenario (S-10) on 27th April 2028.
Figure 16. Depth of phreatic surface of Intervention 01 on 27th April 2028.
Another comparative analysis has been done against extreme Scenario (S-10) with Intervention 01. To understand special variation of groundwater table within the study area and to visualize the improvement of groundwater table or Phreatic surface after application of different interventions, groundwater table at extreme Scenario (S-10) has been deducted from groundwater table at Intervention 01. The deducted values have been specially distributed in the study area where positive values indicate improvement of groundwater level, negative values indicate further lowering of groundwater level (Figure 17).
Figure 17. Changes in groundwater table for replacing Boro with wheat.
Comparison of extreme scenario vs intervention 01, i.e. changes in groundwater table due to replacing Boro with Wheat is shown in Figure 17, where positive values indicate 96.55% improvement of groundwater level, negative values indicate 3.45% further lowering of groundwater level compared to extreme scenario. Details of this analysis have been tabulated in Table 4 and Table 5.
Table 4. Improved area after application of intervention 01.
Difference of GWL (m) |
Percentage of Improved (+ve)/
Deteriorated (−ve) Area |
−0.98 to 0 |
−03.45 |
0 to 2 |
63.05 |
2 to 4 |
27 |
4 to 6 |
5.3 |
6 to 23.94 |
1.2 |
This result shows significant improvement of groundwater resources without any structural construction.
Table 5. Improved area after application of intervention.
Interventions |
Improved Study Area |
Deteriorated Study Area |
Overall Improvement |
Intervention 01 |
96.55% |
3.45% |
93.1% |
According to water table fluctuation method, additional 253 million cubic meter water resources have been added in the study area as groundwater recharge.
4.3. Positive Trend of GWL from 2017 to 2030
After application of Intervention 01, except Tanore Upazila, all of the wells within the model boundary have shown positive trend from 2017 to 2030, which are tabulated in Table 6. For trend analysis of GWL, annual maximum and annual minimum GWL of each station have been analysed (Table 6).
Table 6. Positive trend of GWL within model area.
Well ID |
Max. (m/year) |
Min. (m/year) |
Well ID |
Max. (m/year) |
Min. (m/year) |
GT8134016 |
0.42 |
0.72 |
GT8182042 |
0.05 |
0.3 |
GT8153031 |
0.26 |
0.35 |
GT8172040 |
0.04 |
0.17 |
GT8112001 |
0.16 |
0.27 |
GT6403002 |
0.03 |
0.18 |
GT8134028 |
0.14 |
0.26 |
GT8112002 |
0.02 |
0.31 |
GT6447023 |
0.13 |
0.2 |
GT8131015 |
0.02 |
0.31 |
GT8153030 |
0.12 |
0.37 |
GT8112003 |
0.02 |
0.43 |
GT8125011 |
0.11 |
0.01 |
GT8125009 |
0.02 |
0.12 |
GT8134027 |
0.09 |
0.13 |
GT8134021 |
0.01 |
0.02 |
GT8172037 |
0.08 |
0.11 |
GT8182041 |
0.01 |
0.16 |
GT8125010 |
0.06 |
0.03 |
GT8194046 |
−0.05 |
−0.06 |
GT8134018 |
0.06 |
0.27 |
GT8194044 |
−0.17 |
−0.17 |
GT8172036 |
0.06 |
0.24 |
GT8194047 |
−0.82 |
−0.51 |
GT8182043 |
0.05 |
0.07 |
GT8194048 |
−0.86 |
−0.74 |
GT8125006 |
0.05 |
0.08 |
|
|
|
5. Conclusions and Recommendations
This study employed a coupled MIKE SHE-MIKE 11 modeling framework to evaluate sustainable groundwater management options in Rajshahi District, Bangladesh’s most drought-prone region. Through systematic scenario screening, Scenario 10 (S-10)—characterized by HADCM3-A2 climate projections, extreme dry boundary conditions (10-year return period), and 2.5% annual agricultural demand growth—was identified as the most extreme future scenario, with groundwater levels declining to critical thresholds in six of eight Upazilas by 2030.
Key Quantitative Findings: The replacement of high water-consuming Boro rice (~11,500 m3/ha) with water-efficient wheat (~4000 - 5000 m3/ha) as Intervention 01 yielded substantial hydrogeological improvements:
Areal extent: 96.55% of the study area showed groundwater level improvement, with only 3.45% experiencing further decline.
Magnitude of improvement: 63.05% of the area showed 0 - 2 m improvement, 27% showed 2 - 4 m improvement, and 5.3% showed 4 - 6 m improvement.
Water balance: An additional 253 million cubic meters of groundwater recharge was retained in the aquifer system compared to S-10.
Trend reversal: 22 of 25 monitoring wells (88%) showed positive groundwater level trends (0.01 - 0.72 m/year) under Intervention 01, compared to declining trends under S-10.
These results demonstrate that demand-side management through crop substitution can effectively counteract extreme climate and demand scenarios without structural infrastructure investment.
Study Limitations and Uncertainties: Several limitations affect the interpretation of these findings:
1) Climate Model Uncertainty: As noted, reliance on a single GCM (HadCM3) and emission scenario (A2) does not capture the full range of possible climate futures. Multi-model ensembles would provide more robust uncertainty bounds.
2) Fixed Boundary Assumptions: The assumption of fixed historical boundary conditions for rivers and lateral groundwater boundaries does not account for potential transboundary impacts or regional-scale climate feedbacks. Sensitivity analysis suggests scenario ranking is robust, but absolute level projections carry ±0.5 m uncertainty.
3) Agricultural Demand Projections: The 2.5% annual demand increase represents a business-as-usual scenario without efficiency improvements. Actual adoption of water-saving technologies (e.g., Alternate Wetting and Drying for rice) could reduce demand growth rates.
4) Economic and Social Factors: The biophysical modeling does not address farmer adoption rates, economic viability of wheat versus rice, or food security implications of reduced rice production. These socio-economic dimensions require parallel policy analysis.
5) Model Resolution: The 500 m grid resolution captures regional trends but may smooth local heterogeneities in aquifer properties that affect individual well behavior.
Recommendations:
1) Immediate Implementation: Given the severity of groundwater decline in Tanore, Godagari, Baghmara, and Mohanpur Upazilas (where Boro covers 35% - 61% of land area), pilot programs for Boro-to-wheat substitution should be prioritized in these critical zones.
2) Integrated Approach: Crop substitution should be combined with other demand-management measures (e.g., improved irrigation efficiency, rainwater harvesting) and supply-side options (e.g., managed aquifer recharge) for comprehensive sustainability.
3) Monitoring and Adaptive Management: Continuous groundwater monitoring and model updating (data assimilation) should accompany intervention implementation to verify predicted benefits and enable adaptive management.
4) Future Research: Subsequent studies should i) incorporate CMIP6 multi-model ensembles, ii) couple biophysical models with agricultural economic models, and iii) assess long-term soil health and crop rotation sustainability under changed cropping patterns.
In conclusion, this study provides quantitative evidence that shifting from Boro rice to wheat can improve groundwater sustainability by 93.1% (overall improvement metric) in Rajshahi District, even under extreme climate and demand scenarios. While uncertainties exist in boundary conditions and climate projections, the magnitude of the intervention effect suggests that crop substitution represents a viable, no-regrets strategy for addressing groundwater depletion in the Barind Tract pending a complementary socio-economic assessment.