Assessing climate change impact on root-zone water balance and groundwater levels

India is the largest groundwater consumer in the world (230 km³ year⁻¹), with more than a quarter of the global groundwater (World Bank 2010). In rural areas, groundwater is important as more than 60% of irrigated agriculture and 85% of drinking water supplies depend on it. A report from The National Aeronautics and Space Administration (NASA) showed 109 km³ of groundwater depletion, during 2002–2008 in the Indian States of Rajasthan, Punjab, Haryana and Delhi (Rodell et al. 2009). In Punjab state, there is a demand-supply gap of 1.63 Mm³ of water, which is met by groundwater withdrawal. At the global level, the climate is changing and the surface temperature is predicted to rise within 0.4–2.6 °C in 2046–2065 and 0.3–4.8 °C in 2081–2100 compared with the reference period of 1986–2005 (IPCC 2014). Climate change may bring about many environmental problems, and uncertainties in water resources could be the most severe. The effect of global warming on groundwater is complicated and is dictated by recharge and draft, geographical location, land use and soil type. In addition, changes in rainfall (RF) frequency, temperature, the replacement of either forest or agricultural land with concrete and bitumen expansion of urban built-up areas may also vary the groundwater balance in future. Droughts may result in declining water levels not only because of less RF, but increased soil water evaporation. Paradoxically, extreme RF events may lead to less recharge to groundwater in upland areas because of more runoff. In recent years, groundwater modeling has become a standard tool for understanding groundwater behavior. Several groundwater models such as MODFLOW (McDonald & Harbaugh 1988), FEFLOW (Diersch 2005) and PMWIN (Chiang 2005) are available. Although some location specific studies related to the impact of climate change on groundwater exist in literature (Chen et al. 2004; Green et al. 2011), there are none available for this region. Keeping this in mind, we undertook the present study to assess the impact of regional climate model projections on root-zone water balance and groundwater levels using crop and groundwater models.

The study was conducted at Ludhiana district (75°52′E longitude and 30°56′N latitude) of Punjab State (73°53′–76°55′E longitude and 29°33′–32°31′N latitude), India (Figure 1). The total area of the district is 3,767 km². The altitude varies from 221 m in the west to 273 m in the east above mean sea level. The land slopes from northeast to southwest with an average gradient of 0.38 mm per 1,000 m. Agriculture is the predominant land use in the district (85% of the district area). About 2% of district area is under water bodies and other land uses (marshy land/forest land etc.). The soil texture of the profile varies from sand to silt loam. The lithology of the area is heterogeneous with the presence of many sand beds forming the principal aquifers that are separated by clay beds. The granular material becomes coarser in nature with increasing depth. Groundwater occurs under phreatic as well as confined conditions in these alluvial aquifers.

For the observed climate data, past 40 years (1971–2010) daily data on RF, maximum temperature (T_max) and minimum temperature (T_min) recorded at meteorological observatory of PAU, Ludhiana was collected. The A1B climate scenario, representing an intermediate scenario of two extremes, one for economic development (A2) and other for environment development (B2), was chosen for future climate. For this scenario, climate data on T_max, T_min and RF was obtained from regional climate model, HRM3 at daily intervals for baseline (1961–1990), mid-century (MC) (2021–2050) and end-century (EC) (2071–2098). HRM3 data were available at two longitudes (75°47′E and 76°14′E) for the study area. The average of these two longitudes represented the Ludhiana observatory location. The observed and modeled data on T_max, T_min and RF showed bias, which was evaluated from time trends and statistical parameters like mean (μ), standard deviation (σ) and variance (σ²). To minimize the bias, correction functions were developed from the observed and modeled data of 1971–1985 and validated on 1986–1990. Four types of correction functions, i.e. difference, modified difference, statistical bias correction and modified statistical bias correction at daily, monthly and annual time scales were tested (Jalota et al. 2013).

A stepwise procedure to estimate the impact of climate change on groundwater resources of Ludhiana district is provided in Figure 2. Variation for soil type, land use and drainage in the area was accounted for by digitizing thematic maps, viz. soil texture, land use, canal network and drainage in geographical information system (GIS) software. The digitized canal and drainage map (which was available as a line shape file, Figure 3) were buffered in accordance to their width to form a polygon map. A combined map of soil texture (Figure 4(a)), land use (Figure 4(b)), canal and drainage network was obtained by performing overlay operations. This process resulted in the generation of 24 unique simulation units for the entire study area (Figure 5).

Root-zone water balance was simulated with a crop model (Stockle et al. 1994). The model was set for a rice-wheat cropping system which is the predominant system for central Punjab. This model had already been intensively parameterized, using the experimental data collected at the Research Farm, Punjab Agricultural University, Ludhiana (Jalota et al. 2011). Crop management decisions, soil and weather data were provided as input to the crop model. In the model, weather file for individual year was set from the daily observed data of the present time slice (PTS) (2000–2010) and corrected modeled data for mid-MC and EC. Solar radiation was calculated using the Hargreaves & Samani (1982) method. Data for relative humidity and wind speed were generated by ClimGen model (Stockle & Nelson 1999). The soil files were set using the data presented in Table 1.

In the model, irrigation was scheduled when the ratio of depth of irrigation water (IW) to cumulative potential evapotranspiration (PET) minus RF reached the pre-set values (Prihar et al. 1974). Irrigation water at the spatial scale was scheduled based on the soil texture, where the aforementioned ratios differed. The ratio was increased for coarse textured soil to simulate increased frequency of irrigation and decreased to simulate the effect of decreased irrigation requirement in fine textured soil. The ratios used for rice and wheat crops were 3.0 and 1.0 in loamy sand, 2.5 and 0.9 in sandy loam and 2.0 and 0.8 in loam and silt loam soils, respectively. Depending on the availability, preference was given to canal water. Daily canal water releases for the study area were accounted for and ground water draft was calculated.

In rice, irrigation was applied frequently from transplanting to the first 15 days for proper establishment of the transplanted rice and thereafter the irrigation was scheduled as per the ratio. In wheat, IW/PET ratio was started from the 21st day after sowing in loamy sand and the 30th day in sandy loam, loam and silt loam soils. The nitrogen fertilizer applied to each rice and wheat crops was 120 kg ha⁻¹. It was applied in four splits to rice (a quarter 1 day after transplanting, a quarter at Day 21, a quarter at Day 42 and a quarter at initiation of flowering). Nitrogen fertilizer to wheat was applied in three splits in loamy sand (the first at sowing, the second after first irrigation and the third after second irrigation) and in two splits (first at sowing and second after first irrigation) in sandy loam, loam and silt loam soils. Adjustments in time of fertilizer application were made with respect to changed crop duration and RF in every rotation under the PTS, MC and EC. For each rotation, irrigation schedule based on PET and RF were primed separately and added in the respective management file of the model.

Water balance components for rice-wheat system were simulated using daily weather data of PTS, MC and EC, and CO₂ levels for a given year and A1B scenario as per Bern climate change model (IPCC 2013: www.ipcc-data.org). In the model, the effect of CO₂ on crop biomass production was taken care of by G ratio (a ratio of potential growth at specified CO₂ concentration to reference level of 380 ppm), associated to daily crop radiation use efficiency (Tubiello et al. 2000). Simulations were run for normal planting dates, i.e. June 20 for rice (variety PR 111) and November 5 for wheat (variety PBW 343). Crop duration, irrigation requirement, evapotranspiration (ET) and soil water evaporation during fallow periods were estimated by the crop model. ET from marshy lands was estimated using the models of Kingra & Hundal (2005), and for forest land the models of Minhas et al. (2010) were used.

The study area was discretized with a constant grid spacing of 2.5 × 2.5 km. Hydrologically, the study area is bounded by Satluj River in the north and there is no hydrological boundary on the other three sides. Thus, a constant head boundary was applied to the cells representing the Sirhind canal and Satluj river. Here the storage terms are not used, but the other parameters of Equation (3) are considered. The remaining three boundaries were simulated as flux boundaries using the appropriate recharge/well package. Darcy's law was used to compute the flux from/to an individual boundary cell with the help of Arc Hydro toolkit GIS and Microsoft Excel. The cells lying inside the study area were simulated as active cells for which hydraulic heads were computed throughout all time steps of the simulation. The spatial hydraulic parameters used were initial hydraulic head, hydraulic conductivity and specific yield. The initial hydraulic heads used for model simulation were specified as for pre-monsoon 2000, i.e. on 1 June 2000 (Figure 6(a)). Hydraulic parameters such as hydraulic conductivity and specific yield were estimated indirectly (Todd 1980) from the 43 well logs scattered all across the district (Figure 6(b)).

Twenty years (1971–1990) monthly averages of the observed and HRM3 modeled T_max, T_min and RF at Ludhiana showed that the modeled values of T_max (T_max(modeled)) were higher than that of the observed (T_max(observed)) from February to May; and lower from July to December (Jalota et al. 2013). T_min(modeled) also followed the trend similar to that of T_max(modeled), however the modeled values remained higher than those of the observed from March to October. In the case of RF, the modeled RF values were lower than those of the observed RF from January to May and the trend reversed thereafter up to December.

Statistical parameters, i.e. annual μ, σ and σ² of T_max showed that μ of the modeled and observed were comparable, but σ was 2.3 °C higher in the modeled data. In T_min, μ and σ modeled values were higher by 1 and 2.0 °C, respectively, than those of the observed. In RF, μ of modeled RF was 0.3 mm day⁻¹ more than that of the observed and σ for the same was 4.1 mm day⁻¹ less. In the present study bias correction method was applied and the best fitted bias correction functions for T_max, T_min and RF were used, which are given elsewhere (Jalota et al. 2013). These functions made the time trends similar to those of the observed in addition to lowering of root mean square error in corrected modeled data. RMSE of the corrected modeled T_max and T_min was 7 and 12%, respectively, which fall in the excellent to good category. However, RMSE of model corrected RF remained high (67%) even after correction.

Average annual T_max of 30.1 ± 0.5 °C in PTS would increase to 32.1 ± 1.0 °C in MC and 34.8 ± 1.2 °C in EC. Similarly, the T_min of 17.5 ± 0.2 °C of in PTS would increase to 19.0 ± 0.7 °C (14%) in MC and 21.6 ± 0.7 °C (30%) in EC. Thus, for A1B scenario, the mean annual temperature would increase by 1.8 °C in MC and 4.4 °C in EC, compared to that of the PTS. The RF showed an increasing trend with prominent inter-decadal variability. The RF of 728.0 ± 199.4 mm in PTS would increase to 1055.7 ± 144.7 mm (45%) in MC and 1213.8 ± 239.7 mm (67%) in EC.

The area under agricultural and built-up lands would change with increase in population. The present area of 10% of Ludhiana district which is built up would increase to 12, 13, 16 and 20% by 2021, 2031, 2041 and 2051, respectively. Thereafter, the increase in area which was built up was not considered, due to stabilization of population as characterized in the A1B scenario. With increased temperature in MC and EC, the crop model simulations showed shortening of crop duration (Jalota et al. 2013), which would alter the root-zone water balance components in different land uses (Table 2).

In volumetric terms, averaged across years and land uses, the RF over the area in the PTS is 2,663 Mm³, which would increase by 658 Mm³ (25%) in MC and 1,211 Mm³ (45%) in EC. The IW requirement of 4,199 Mm³ in PTS would be reduced by 1,744 Mm³ (42%) in MC and 2,045 Mm³ (49%) in EC. Similarly, evaporative loss (from agricultural, forest and marshy lands) of 3,631 Mm³ in PTS would reduce by 203 Mm³ (6%) in MC and 773 Mm³ (21%) in EC annually. The reduction in evaporative loss is primarily due to reduction in transpiration from agricultural lands with a shortening of crop duration, although loss from soil and water bodies was more. Drainage beyond root-zone in PTS is 3,394 Mm³, which would reduce by 685 Mm³ (20%) in MC and 260 Mm³ (8%) in EC. Runoff of 302 Mm³ in PTS would increase by 84% in MC and 190% in EC. In future, increased RF, decreased crop land and more release of water from canals to avoid floods are the factors responsible to stimulate run off. The magnitude of water balance components in different land uses would also change in future. Volumetric runoff would decrease from agricultural lands and increase from built-up lands in response to the projected land shrinkage under agriculture and expansion in built-up land. Evaporative loss in the future would decrease from agricultural lands and increase from water bodies. Runoff in the future would increase from built-up lands.

The potential groundwater recharge of 3,394 Mm³ in PTS would change to 2,810 Mm³ in MC and 3,519 Mm³ in EC, thus indicating 17.2% decrease in MC and 3.7% increase in EC. Decreased transpiration losses and crop irrigation requirements due to changed land use in MC and EC as compared to PTS, reduced the groundwater draft from the aquifers.

Groundwater draft of 4,122 Mm³ in PTS would decrease by 1,698 Mm³ (41.2%) in MC and 1,851 Mm³ (44.9%) in EC. The changes in recharge and draft results in change in the hydraulic head values, simulated using MODFLOW. The model showed reasonable matching with the observed values during calibration and validation (Figure 7). Simulation error calculated as NRMSE value ranged from 0.011 to 0.016 m for the calibration period and from 0.012 to 0.018 m for the validation period, which is within the accepted calibration accuracy of 10%. Averaged over 5 years, groundwater level showed a fall of 0.63 m year⁻¹ during 2000–2005, which decreased to 0.43 m year⁻¹ during 2005–2010 in PTS. In MC, groundwater levels would rise and the maximum annual rise of 0.41 m year⁻¹ was predicted for the period 2035–2040, followed by 0.40 0.24, 0.11 and 0.09 m year⁻¹ for the periods 2030–2035, 2021–2025, 2025–2030, 2040–2045 and 2045–2049, respectively (Figure 8). The annual rate of change in water levels showed fluctuating trends viz. 7 years fall, 2 years of no change, with rest experiencing rise. The average rise of water level would be 0.23 m year⁻¹ in MC.

The present study gives a methodology to put together the outputs of climate (HRM3), cropping system (CropSyst) and groundwater (MODFLOW) models using GIS for projecting groundwater in future at local/regional scale with varying land uses. In previous studies (Mall et al. 2006; Kumar et al. 2011), future climate was predicted in relation to the modeled climate data of base line (1961–1990) without minimizing bias. Many studies in the literature have demonstrated that biases exist in the regional climate models data and require correction (Leander & Buishand 2007; Sennikovs & Bethers 2009; Haerter et al. 2010; Piani et al. 2010a, 2005b). Likewise, for estimating groundwater recharge in response to climate change, physical based models have been used by various researchers (Jyrkama & Sykes 2007; Ficklin et al. 2010), which were incapable of accounting for boundary conditions, changed CO₂ concentration as a result of climate change, crop phenology, etc. Moreover, in previous studies a steady-state groundwater model has been used to simulate the effects of climate change on groundwater (Woldeamlak et al. 2007). The present study has an edge over the previous studies as in this case the crop duration, yield and root-zone water balance were assessed in different time periods of the 21st century with intensively calibrated and validated cropping system (CropSyst) model using the corrected future climate data, actual soil profile data, management intervention and crop data (Jalota et al. 2009, 2011, 2013). In this model CO₂ concentration is also accounted for (Tubiello & Ewert 2002) in addition to other weather parameters like temperature, relative humidity, RF, solar radiation and wind speed. For groundwater, transient simulations were run through a customized groundwater model (MODFLOW).

The projected climate data showed that T_max, T_min and RF in the future would increase as compared with that of the PTS. Earlier research (Kumar et al. 2011; Chaturvedi et al. 2012) also made such projections. In future, though RF and temperature would increase, increased temperature would shorten crop duration (Jalota et al. 2012, 2013), which would directly affect the components of root-zone water balance, i.e. irrigation applied, evapotranspiration, runoff and drainage; and groundwater components, i.e. recharge and draft. In addition to these, it is also important that RF events do not necessarily correspond to an increase in recharge as anticipated, but is dependent on many other factors, for example, soil type, land use, timing of RF event, local geology and topography. The greater magnitude of RF and increase in urbanized area would contribute significantly to surface runoff. Decreased groundwater pumping due to shortening of their duration and more RF in different years would show fluctuating trends of water levels in future, i.e. rising trends in more numbers of years. There are uncertainties in climate, crop and groundwater models, yet the authors hope that this study gives a sound base for further research on understanding and refining the impact of climate change on water resources estimation and management in the future.

The authors acknowledge the Indian Council of Agricultural research (ICAR) for financing the study under All India Co-ordinated Research Project on Groundwater Utilization and Emeritus Scientist Scheme. Thanks are also due to the Hadley Centre for Climate Prediction and Research, UK Meteorological Office and Indian Institute of Tropical Meteorology, Pune, for providing HRM3 data set, and to Dr Anil Sharma, Additional Director, Centre for Communication & International Linkages, PAU, for his English editing services.