Studying the impact of climate change on spatiotemporal variability of blue and green water resources

Water is an essential natural resource for human beings, which has significant effects on the sustainable development of the environment and the stability of manufacturing activities. Nowadays, water resources all around the world are becoming increasingly vulnerable due to the increasing levels of societal demand. On the other hand, given the global warming and the rising frequency and intensity of extreme weather, the effects of climate change on water resources are serious. Various river-basins all over the world are expected to experience undesirable effects of climate change on water resources and freshwater ecosystems (Abbaspour et al. 2009). Therefore, depicting the effects of climate change on various components of the water cycle is of great significance in the management of water as an indispensable resource.

Most hydrological systems consist of very complicated and heterogeneous processes which are not easily understood (Teshager et al. 2016). Furthermore, it is extremely difficult to measure all inputs, outputs, and states of hydrological processes and/or parameters in spatial and temporal scales due to spatial heterogeneity, temporal dynamics of hydrological forcing functions and processes, deficiencies in measurement methods, and limitations in times and costs (Mengistu et al. 2019). Knowledge about the internal renewable water resources provides valuable information in long term management and planning. Falkenmark (1995) introduced two main types of water: blue water (BW) and green water (GW). BW includes surface and groundwater runoff, and GW refers to precipitation that is stored in the root zone of the soil, which evaporates, transpires, or gets incorporated by plants.

Distributed process-based hydrological models are suitable tools for determining the spatial and temporal information of the watershed area; they help to understand the hydrological processes and/or interactions between the watershed features and the hydrological responses, and support sustainable water resource planning, management, and decision-making procedures (Teshager et al. 2016; Franco & Bonumá 2017; Mengistu et al. 2019). Over the recent decades, extensive researches have been conducted to study the spatial and temporal variations of the hydrological components, BW, and GW indicators, particularly in future climate change scenarios. The impacts of future climate on Iran's water resources were investigated using the calibrated Soil and Water Assessment Tool (SWAT) model for 1980–2002 time horizon at a sub-basin level by Abbaspour et al. (2009). The future climate scenarios, A1B, B1, and A2, were generated by the Canadian Global Coupled Model (CGCM 3.1) and downscaled for 37 climate stations. The effects of future climate scenarios on BW, GW, and wheat yield across the country were analyzed with the SWAT model. Rodrigues et al. (2014) illustrated analysis of the availability and use of BW and GW to represent water scarcity and vulnerability indicators at the basin scale in the Cantareira water supply system in Brazil. The SWAT model was applied to depict the hydrological processes and derive the BW and GW footprint indicators against various water access levels used for human activities during 23 years. The contrasting status of BW indicators was studied in their research based on different hydrological based methodologies to determine the monthly environmental flow requirements (EFRs), and the risk of natural EFR violation. Farsani et al. (2019) evaluated the effects of climate change on the spatiotemporal distribution of water resources in order to analyze the water supply demand in the Bazoft watershed, Iran. The SWAT model was applied to depict the changes of BW flow, GW flow, and GW storage for a future period (2010–2099) in comparison to a historical record (1992–2008). Veettil & Mishra (2018) studied the effects of both anthropogenic and climatic features on the spatiotemporal variability of water security indicators, such as BW scarcity, GW scarcity, Falkenmark index, and freshwater provision indicators in the Savannah River basin. Their results implied that the study area experienced a decline in BW due to the climate variability whereas GW was significantly affected by land use alternation. Zhu et al. (2018) examined the spatiotemporal distribution of GW in the Hai River basin based on the impacts of land use types. SWAT model was employed in their research for the determination of the relation between the soil and land use type with certain indices, such as the maximum possible storage of GW, the GW footprint, and the available GW. The impacts of climate change on runoff, aquifer infiltration, renewable water resources, and drought intensity in the Salt Lake sub-basin, Iran, were assessed by Khalilian & Shahvari (2019). The calibrated and validated SWAT model projected the hydrological responses in watershed scale in a spatiotemporal framework based on various climatic scenarios. Huang et al. (2019) conducted a research to study the effects of climate and land use changes on crop GW and BW consumptions. A crop water use module, developed based on the global change assessment model and its hydrology module (Xanthos), was used to illustrate the effects of climate and land use changes on BW and GW footprints. The modeling results demonstrated the increase in global crop green water footprint and dominance of climate change over land use change on the GW footprint. Furthermore, the global crop BW footprint would increase, specifically in areas with considerable irrigated land development.

Liang et al. (2020) determined the spatiotemporal variations of GW and BW at the rapidly developing Xiangjiang River basin in China according to climate change and anthropogenic activities. SWAT, as the semi-distributed and process-based hydrological model, demonstrated the influence of (1) BW scarcity from precipitation and population growth and (2) GW scarcity from agriculture and urban lands. In other words, the climatic features were found to affect the BW components in the river basin scale. Furthermore, the spatial shortages of BW and GW in the river basin scale showed higher values in downstream reaches of the watershed. Mao et al. (2020) studied the interactions between GW and BW, focusing on the anthropogenic activities in an arid endorheic river basin in China. The knowledge about the interactions between GW and BW in the hydrological system was applied in integrated surface water and groundwater resource management and support of basin scale water resource management, addressing human nature water conflicts. A regional climate model (RCM), the COSMO Climate Limited Area Model (CCLM), linked to the SWAT model, was applied by Mengistu et al. (2021) to assess the effects of climate change on the Upper Blue Nile (Abay) water system. They illustrated an increase in potential evapotranspiration (PET) by 27% under RCP 8.5, an increase in surface runoff by 14%, and a decrease in base flow by 30% in comparison to the baseline scenario. Despite the increase in the surface runoff, the total water yield in the study area was estimated to decrease by −1.7 to −6.5% and −10.7 to −22.7% in RCP 4.5 and RCP 8.5 scenarios, respectively. Serur (2020) simulated BW and GW resources availability at the basin and sub-basin levels in the Weyb River basin in Ethiopia using the SWAT model. According to their results, the mean annual BW flow, GW flow, and GW storage exhibited an increase in the entire basin and in all the sub-basins under representative concentration pathway (RCP) 8.5/4.5/2.6 scenarios.

To our knowledge, the main attempt in modern water governance in the watershed scale is to assess the effects of future climate change on water cycle. The scientific numerical hydrological models can help to define policy for the stakeholders and policy makers, which would result in more effective use of soil and water resources. In view of the changing pattern in climate features, it is desirable to investigate the spatiotemporal variability of water footprint indicators at river basin scales for water resource management and ecological conservation, especially in fragile watersheds. In this research, the integrated hydrological SWAT model (Arnold et al. 1998) was calibrated and validated (time horizon 2003–2016) to study the effect of climate change at a basin level for the Seymareh River Basin (SRB), Karkheh, Iran. We investigated the changes of various components of water balance, including precipitation and evapotranspiration (ET) distribution, stream flow, soil moisture, and aquifer recharges. The spatiotemporal variability in water resources concerning BW and GW indicators was derived under RCP 8.5/6.0/4.5 scenarios. Limited researches have been conducted in identifying areas with high BW/GW potentials to plan the rain-fed or irrigated agriculture in the SRB scale. In addition, the current study would be effective in enhancing socio, economic, and environmental capital and sustainable development aspects in the SRB.

The subsequent sections of this paper are organized as follows: Section 2 introduces the study area, the required data, hydrological model characteristics, setup, and calibration, future climate data and modeling scenarios, and green and BW delineation. Section 3 begins with the calibration and validation results of the SWAT model. These results are demonstrated within the downscaling climate variables and accompanied by an analysis of the climate change impacts. The manuscript concludes with Section 4, by summarizing key findings and suggesting directions for future research.

Seymareh River leaves the Gamasiab sub-basin, flows through Gharasou and then the Seymareh sub-basin, through a series of deep valleys and canyons and joins the Seymareh Reservoir. Gamasiab, Gharasou, and Seymareh sub-basins also contribute about 37, 23, and 43% of the annual flow of Seymareh River. Roughly 9,500 km² of the SRB's total surface area is arable land and is devoted to horticulture (apple, grape and walnut) and field crop (wheat, barley, and sugar beet) production. About 1,700 km² of this land is irrigated using traditional and modern techniques, and 7,800 km² of it is rain fed. The remaining land use of the SRB includes forests, pastures, rocks (totally 64%) and barren lands, wetlands, and urban areas (less than 2%) (MGCE 2014). For modeling of the SRB, most (about 90%) horticulture and field crop productions are considered, such as apple, grape and walnut, wheat, barley, and sugar beet.

To address the current project objectives, the required datasets were collected from multiple sources. The fundamental component of the datasets used for model development included (1) the digital elevation model (DEM), obtained from the global U.S. geological survey's (USGS) public geographic domain database set at a spatial resolution of 90 m. The DEM was utilized to delineate the study area's topographic characteristics; (2) the land use data sets were provided from USGS EROS Archive – Land Cover Products – Global Land Cover Characterization webpage (USGS EROS 2019); (3) the soil data were extracted from FAO/UNESCO (2019) global soil maps of the world; (4) the daily meteorological (precipitation and the minimum and maximum temperatures) data from 2003 to 2016 were downloaded from National Centers For Environmental Prediction (NCEP) climatic SWAT data website (https://globalweather.tamu.edu/); (5) the recorded stream flow data were obtained from the Basic Study Office of Iran water resources management company (IWRMC 2017) for seven hydrometric stations from 2003 to 2016; and (6) historical records on the annual yields of horticulture and field crop productions were collected from 2000 to 2018 from the agricultural statistics and the information center of Ministry of Jahade-Agriculture (Ahmadi et al. 2019).

The SWAT is a continuous, long term, physically-based, and semi-distributed hydrological model. SWAT model is able to assess the effects of climate and land management on the hydrological processes, sediment loading, pollution transport, and crop growth in watershed scales (Arnold et al. 1998). The SWAT model is useful for estimating BW and GW available at a basin scale (Schuol et al. 2008; Abbaspour 2015; Veettil & Mishra 2016).

In the SWAT model, a watershed is partitioned into sub-basins that are further divided into a series of hydrological response units (HRUs). HRUs are uniform units that share unique combinations of soil and land use and are employed as the basis of water balance calculation. The water components, sediment yield, and nutrient cycles are estimated in each HRU and then aggregated for the whole sub-basins in the watershed.

The SRB was divided into sub-basins, and then partitioned into unique HRUs. Five classes of slopes were defined for HRU delineation: 0 –5%, 5–10%, 10–20%, 20–40%, and >40%. The number of HRUs was determined by adjusting the threshold of land use (5%), soil (10%), and slope (10%), which resulted in 1,723 HRUs under land use features of the study area. In this research, the surface runoff was illustrated with soil conservation service curve number (SCS-CN) using daily meteorological data and soil hydrologic group, land use and land cover features, and antecedent soil moisture. The PET and actual evapotranspiration (AET) were simulated based on the Hargreaves Method and the Ritchie Method, respectively. The leaf area index (LAI) and root development were simulated using the “crop growth” component of SWAT, which is the simplified version of the erosion productivity impact calculator (EPIC) crop model (Williams et al. 1984). Plant growth was determined according to leaf area development, light interception, and conversion of intercepted light into biomass, assuming plant species specific radiation use efficiency. Phenological plant development was based on daily accumulated heat units, potential biomass, and harvest index. Plant growth could be inhibited by user specified temperature, water, nitrogen, and phosphorus stress factors (Neitsch et al. 2002).

Simultaneous past and future events of climate change are affected by a combination of external forcing, unforced internal fluctuations, and the response characteristics of the climate system. Due to human induced greenhouse gas (GHG) emissions, the IPCC expects the average global temperature to change in future years (IPCC 2014). It is crucial to understand how much change to the earth's climate system will affect future hydrological cycles and water components. To study the impacts of future climate changes on hydrological cycles, scenarios are created based on global circulation models (GCMs), considering certain boundary conditions (such as the solar constant) or physical parameters (such as the GHG concentration). To evaluate the effects of climatic scenarios on future (2020–2040) hydrologic events, GCM data have been used as an input to hydrologic models. However, such data cannot be applied directly since hydrologic models need local scale daily meteorological data as input (IPCC 2014). Therefore, the GCM outputs should be converted to the appropriate spatial and temporal resolutions using statistical downscaling tools, the LARS-WG model, for instance.

According to the IPCC-V assessment report, cumulative emissions of CO₂ mainly affect global mean surface warming by the late 21st century and beyond. Projections of GHG emissions alter over a wide range, depending on both socio-economic development and climate policy. In this regard, four different 21st century pathways of GHG emissions and atmospheric concentrations, air pollutant emissions, and land use have been considered in the representative concentration pathways (RCPs) (IPCC 2014). In this study, we analyzed the effects of two intermediate scenarios (RCP 4.5 and RCP 6.0) and one scenario with very high GHG emissions (RCP 8.5) on the hydrological cycles.

LARS-WG, as a credible stochastic weather generator tool, could be used for the projection of weather data at a local site, under both current and future climate conditions. This tool provides a means to extend the weather time series simulation of unobserved locations through the interpolation of the weather generator parameters obtained from running the models at neighboring sites (Semenov 2021). This model can serve as a computationally inexpensive tool to generate multiple year climate change scenarios at the daily time scale, which considers changes in both mean climate and climate variability. The LARS-WG 6.0 model incorporates climate projections from the CMIP5 ensemble described in the IPCC Fifth assessment report. The model incorporates the process of generating synthetic weather data into three distinct steps: model calibration, model validation, and generation of synthetic weather data (Semenov 2021).

While using LARS-WG, model calibration consists of calculating the relevant statistical parameters for each meteorological variable from the observed historical data. These parameters or the once modified ones based on future climate change scenarios are then used to stochastically generate realistic climate data corresponding to the present or future climate scenario, respectively. For the first set of experiment, the mean of observed daily precipitation as well as daily maximum and minimum temperatures are used to extract the statistical parameters of the current climate. For precipitation, these parameters consist of monthly histogram intervals and frequency of events in each interval for dry and wet spell lengths, as well as precipitation amounts. On the other hand, temperature is modeled in LARS-WG by using Fourier series which can be constructed with parameters such as the mean value, amplitude of the sine and cosine curves, and phase angle. Both maximum and minimum temperatures are modeled more accurately by considering the wet and dry days separately; therefore, the temperature parameters for wet and dry days are derived separately. The weather generator also uses parameters corresponding to the average autocorrelation values for minimum and maximum temperatures derived from the observed weather data. After the observed weather data are analyzed in this way, the derived statistical parameters are used to generate synthetic weather data representing the current climate (Dibike & Coulibaly 2005).

In this research, the ArcView GIS interface for SWAT 2012 (Winchell et al. 2013) was configured and parameterized to depict the hydrological processes in the SRB. In the SWAT model set up, the watershed was delineated into 33 sub-basins with the main outlet in Seymareh Reach inflow. The daily precipitation data and the minimum and maximum temperature data were obtained from NCEP climatic SWAT data website. The hydrologic SWAT model was calibrated and validated at the sub-basin level based on the monthly observed discharges at seven stations across the river basin, and annual crop and garden production yields. The combinations of river discharge and crop and garden yields in the calibration processes result in a more accurate approximation of both runoff and ET and, therefore, soil moisture and deep aquifer recharge. The SWAT calibration and uncertainty program (SWAT-CUP) was applied to calibrate the SWAT model parameters. Adjusting the model parameters would result in model accuracy to better depict the hydrological process of the study area. In this research, the SWAT simulation of SRB consists of the warming up period (2003–2004), the calibration period (2005–2012), and the validation period (2013–2016). As the SWAT model involves extensive parameters, a sensitivity analysis was accomplished to identify the main parameters across various hydrologic regions.

The SUFI-2 algorithm in the SWAT-CUP program was used in the SWAT model parameter adjustion. The whole uncertainties (parameter, conceptual model, input, and so forth) of the simulation were included in the parameter ranges as the algorithm attempts to involve the 95% prediction uncertainty of the measured data. Two indices quantify the goodness of calibration/uncertainty performance, the P-factor, and the R-factor. In order to compare the observed and simulated monthly discharges and annual crop yields, the coefficient of determination (R²), Nash–Sutcliffe (NS) efficiency coefficients, and percent bias (BIAS) were used as the criteria to evaluate the accuracy of SWAT model performance.

‘Blue water’ is generally defined as ‘the sum of river discharge and deep groundwater recharge’ (Abbaspour et al. 2009). Based on the modeling framework, BW is estimated by combining both water yield and groundwater storage. The amount of water leaving the HRU and inflowing to the main channel is illustrated as water yield. Groundwater storage is defined as the difference between the total amount of water recharging to aquifers (GW_RCHG) and the amount of water from the aquifer entering the main channel flow (GW_W). ‘Green water’ consists of resource and flow as soil moisture and actual evaporation–transpiration, respectively (Abbaspour et al. 2009; Veettil & Mishra 2016). In this study, these two different sources of water were determined based on the SWAT model results and the effects of climate changes on each source were evaluated.

Due to extensive tunable parameters in the SWAT model, a sensitivity analysis was accomplished to identify the main parameters across various hydrologic regions. An automatic calibration with the SUFI-2 algorithm involved in the SWAT-CUP program was applied to calibrate the model and find out how sensitive is the calibration to the tunable parameters. The SUFI-2 algorithm analyzed and optimized more significant parameters, representing the processes of surface and subsurface runoff generation, as well as the annual yields of horticulture and field crop production.

The SWAT model is based on a multitude of physical parameters of mathematical equations. The sensitivity analysis method helps in the (a) avoidance of over-parameterization, and (b) determination of the most influential parameters on the hydrological process. The t-stat and p-value, as statistical measures, determine the sensitive rank of each parameter. The lower the p-value and the higher the absolute value of t-stat, the more sensitive the parameter. The parameter selection scheme for the setup of physical-processes-based distributed-parameter hydrological model is determined according to the sensitivity analysis technique. The model performance achieved the satisfactory level after the main parameter's calibration. Previous studies indicate that suitable results, depicting the real world observed data, are obtained based on a large number sets instead of one best parameter set. The output of the assessment of the sensitivity procedure on the surface runoff indicates that SOL_BD, ESCO, and CN2 are the most influential parameters. In addition, based on previous studies, other important parameters were identified as the significant parameters for SWAT simulation in the SRB (Noori et al. 2017). Table 1 represents the p-value, as the statistical measure, in the sensitivity analysis procedure for streamflow simulation. Having accomplished the sensitivity analysis, SWAT model calibration was then conducted according to the selected sensitive parameters. The adjusted values for each parameter considered in the calibration process are presented in Table 1.

The current study indicates that the sensitivity rank of CN2 is lower than SOL_BD and ESCO on streamflow, representing the occurrence of more complex hydrological processes in the SRB. The dominant performance of SOL_BD is due to the combination of lateral and vertical water movements in loess soils. The lateral and vertical redistribution of water movement in the soil layers effectively affect the lateral flow, groundwater, and ET of the hydrological processes.

A deeper comprehension of both hydrological cycle and model simulation results is achieved by uncertainty analysis. The larger the P-factor value, the greater the contribution of parameter uncertainty to the simulation uncertainty. These results implied that adjusting appropriate parameters makes it possible to represent model simulation uncertainty and attenuate the input and/or structure uncertainty. Table 2 demonstrates the statistical analysis on SWAT model simulation results in the calibration and validation periods at seven hydrometric stations, located in the lower, middle, and upper parts of the SRB. The comparison of the simulation results and the observed streamflow data at the hydrometric stations, as presented by R², NS (Nash–Sutcliffe), P-factor, and R-factor, may be taken into consideration as an indicator of good calibration and validation of the model.

Given the spatial distribution of these hydrometric stations, the results (Table 2) illustrated that the R² value in the downstream station was lower than that of the upstream station. The model produced results with an average P-factor of 0.71 and an R-factor of 1.23, during calibration and validation, 0.70 and 1.24, respectively. Although, most of the observed recorded data were within the 95PPU boundary, the P-factor measures were within no appropriate bound of P-factor > 70% (Abbaspour 2015) in the limited hydrometric stations (Table 2) in the calibration and validation periods. These statistical measures imply that the SWAT model can be applied as a suitable decision support tool for water management and planning in the SRB.

The downscaled meteorological data from LARS-WG agreed relatively well with the recorded historical data in all the 11 control stations. The downscaled temperature in all the meteorological stations had R² values in the range of 0.96–0.99. The R² value and the NS coefficient, calculated between the downscaled rainfall data (simulation with the LARS-WG model) and the measured data, indicated entirely satisfactory results for all the regions in the study area.

The differences between the historical and projected precipitation with different emission scenarios were calculated and compared for the period of 2020–2040 and 2003–2016, respectively. According to the obtained findings, all climate change scenarios indicated an increase in precipitation in the study area in the range of 29–33%.

An analysis of the hydrological status under climate change scenarios of the IPCC-V assessment report revealed intensification of surface runoff, especially in the wet seasons of the year in the SRB. However, the surface runoff intensification in the upstream parts of the SRB seems more evident than that in the downstream parts (Figure 3).

The long term average spatial changes of GW contents in various sub-basins in RCP 4.5 and 8.5 scenarios were −7.4% and −1.9%, respectively. Furthermore, the standard deviations of the spatial changes of GW contents in various sub-basins in RCP 4.5 and 8.5 scenarios were 0.84 and 0.18%, respectively. According to the results (Figure 8) and the spatial analysis, the reductions of GW contents in all the sub-basins of the SRB occurred in both RCP 4.5 and 8.5 scenarios.

In this research, the spatiotemporal impacts of climate change on the hydrological components, precipitation, GW, and BW, were studied utilizing the SWAT model. The knowledge on the hydrological components, as an essential subject in the river basin water resources management and planning, was derived based on the physically-based distributed hydrological model. The SWAT model in the SRB was calibrated and validated based on the measured streamflow in seven hydrological stations and the yields of agricultural and horticultural products in the period 2003–2016. The spatial analysis on precipitation, GW, and BW demonstrated further water availability in the western parts of the SRB.

Future climate scenarios for the period of 2020–2040 were generated based on the LARS-WG 6.0 model, which projected the future climate status from the CMIP5 ensemble described in the IPCC Fifth assessment report. In the climate change scenarios, higher flow discharges are experienced in the hydrological stations, especially in the wet seasons, whereas in the dry seasons, they are extended. According to the results, the climate change scenarios influenced precipitation more substantially compared to the other meteorological components in the SRB. An analysis on the meteorological components demonstrated that in the climate change scenarios, the projected annual precipitation and ET would increase and decrease, respectively, relative to 2003–2016 as the baseline period. The significant spatiotemporal changes of BW were principally anticipated due to the climate change scenarios whereas the effects of the climate change scenarios on GW were insignificant. In other words, BW was mainly affected by precipitation and positively correlated with it.

This study would help regional watershed managers and national policy makers to make more accurate and reasonable decisions concerning water resources planning and protection in the SRB. Future research is recommended to include more water quantity and quality indicators as water security indicators. Moreover, studying further effects of future land use and the newer IPCC report scenarios on the hydrological components in the corresponding HRUs of the SRB would be interesting; this might help to get insights into the spatial variations of each variable.

Calibration and validation of the SRB SWAT model based on remote sensing production (i.e. LAI, soil moisture contents, and snow cover fraction), in addition to the current field data, could depict a more accurate projection of the hydrological and agronomical processes in the SRB. Limited studies have been accomplished in the projection of spatial rain-fed or irrigated agricultural distributions in the SRB scale. Furthermore, the spatiotemporal mapping of GW and BW resources will help the planners and authorities in persuading farmers to adopt cropping patterns in accordance with the territorial capacities.

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

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

The authors declare there is no conflict.