Assessing Long-Term Hydrological Impact of Climate Change Using an Ensemble Approach and Comparison with Global Gridded Model-A Case Study on Goodwater Creek Experimental Watershed

Abstract

Potential impacts of climate change on the hydrological components of the Goodwater Creek Experimental Watershed were assessed using climate datasets from the Coupled Model Intercomparison Project Phase 5 and Soil and Water Assessment Tool (SWAT). Historical and future ensembles of downscaled precipitation and temperature, and modeled water yield, surface runoff, and evapotranspiration, were compared. Ensemble SWAT results indicate increased springtime precipitation, water yield, surface runoff and a shift in evapotranspiration peak one month earlier in the future. To evaluate the performance of model spatial resolution, gridded surface runoff estimated by Lund–Potsdam–Jena managed Land (LPJmL) and Jena Diversity-Dynamic Global Vegetation model (JeDi-DGVM) were compared to SWAT. Long-term comparison shows a 6–8% higher average annual runoff prediction for LPJmL, and a 5–30% lower prediction for JeDi-DGVM, compared to SWAT. Although annual runoff showed little change for LPJmL, monthly runoff projection under-predicted peak runoff and over-predicted low runoff for LPJmL compared to SWAT. The reasons for these differences include differences in spatial resolution of model inputs and mathematical representation of the physical processes. Results indicate benefits of impact assessments at local scales with heterogeneous sets of parameters to adequately represent extreme conditions that are muted in global gridded model studies by spatial averaging over large study domains.

1. Introduction

Climatic shift due to anthropogenic activities has been linked to water supply shortages [1,2], declining biodiversity [3,4,5], ecosystem damage [6], and economic impact [7]. Increased greenhouse gas (GHG) concentrations alter the radiative balance of the Earth’s atmosphere, causing an increase in average temperature (T) and changes in precipitation (P) patterns [8,9]. The Intergovernmental Panel on Climate Change (IPCC) reports increases in mean surface temperature of between 0.3 and 4.8 °C by 2081–2100 relative to 1986–2005, and variable changes in precipitation across the globe for the same period of comparison [10], dependent upon GHG concentrations in the atmosphere. The atmospheric carbon dioxide (CO₂) concentration is estimated to increase from the present day concentration of approximately 400 ppm to approximately 850 ppm by year 2075 in the highest emission scenario (RCP 8.5) [10]. The projected changes in meteorological variables have been reported to impact hydrological components, and ultimately agroecosystem functioning [11,12,13]. The majority of modeling exercises have been conducted at larger scales than this study, despite the fact that in reality, management decisions are taken at a local scale by individual farmers. Therefore, it is important to develop approaches to increase the ability to predict how future climate changes could impact hydrological components, and correspondingly, agricultural productivity, through the representation of localized systems.

Multiple Global Circulation Models (GCMs) have been developed to project future Earth climate given multiple plausible futures [10]. GCMs are three-dimensional numerical models that represent the physical dynamics of the atmosphere, ocean, cryosphere, and land surface, and are the best available methods for predicting how these systems behave given increasing GHG concentrations in the atmosphere [14]. In the Fifth Assessment Report, the climate projections are organized into four Representative Concentration Pathways (RCPs) (RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5) based on radiative forcing levels associated with socioeconomic and technical assumptions by the end of the century [15]. Radiative forcing is defined as the difference between incoming solar radiation and outgoing infrared radiation caused by increased concentration of GHGs expressed in Watts per square meter (W/m²) [16]. Note that climate change impact assessment studies using hydrological models include different sources of uncertainty: uncertainty in future climate projections (GCM), downscaling techniques, and hydrological modeling [17,18].

Hydrological models have been used to simulate the long-term impact of projected changes in climate to assist with the quantification of future risk and, if necessary, to assist with developing adaptive management strategies [19,20]. Hydrological simulation models have been applied to assess the impacts of climate change on nutrient loading [21], crop productivity [22], and streamflow [23,24]. The majority of climate impact assessment studies have been conducted at larger scales using gridded model structures matching with the resolution of Coupled Model Intercomparison Project (CMIP3 and CMIP5) [22,24,25,26]. Studies in the Midwestern United States have focused on quantifying the impact of future climate on hydrological processes at a regional scale using gridded climate data (50 × 50 km) [18,27,28] and with some studies conducted in mid-sized watersheds [29,30,31]. There have been relatively few studies at smaller scales [32,33,34], in part due to lack of fine resolution climate model output and corresponding, data-intensive requirements to further downscale these data to a specific location to use in hydrological models for local scale applications [19]. Watershed sizes were defined using United States Geological Survey (USGS) classification scheme based on hydrological unit code (HUC) which consists of six levels of classification [35]. Specifically, subwatershed to watershed (~100 km² to 585), subbasin to basin (>585 to 27,445 km²), and sub region to region (>27,445 km² and above) were considered as local scale, mid-size and regional scale, respectively.

Hydrological components can be simulated at various scales. The spatial and temporal scale of a hydrological model also defines the resolution of key output variables. For regional reservoir managers, a coarser resolution of both space and time may be sufficient for anticipating and planning for predicting water levels. For example, the Variable Infiltration Capacity (VIC) model [36] operates at a 10-kmgrid cell size and a daily time-step, and is useful for regional-scale water management planning and identification of general trends. In contrast, the Agricultural Policy Environmental eXtender (APEX) model [37] operates at field scale (a few to few hundreds of hectares), and is useful for field-scale management planning. Thus, hydrological models can produce distinct results, because of the model structure, parametrization and resolution of inputs [38,39]. Comparison of model results at multiple scales can help to evaluate which resolution is necessary and appropriate for particular risk and adaptation questions. For example, early work on hydrological impacts associated with future climate change have focused on determining which regions of countries and even the globe may be at risk of changes in water availability, a question for which a coarse resolution may be adequate. However, in order to develop strategic responses in a specific location, and particularly the need to link hydrology to agricultural productivity, a finer resolution may be required in order to evaluate low probability, high risk events (e.g., flooding, drought, etc.). Similarly, a finer resolution for key parameters representing physical processes is likely to result in an increased ability to predict responses to future climate conditions.

The geographic focus of this study is the Goodwater Creek Experimental Watershed (GCEW), which is part of a Long Term Agroecosystem Research (LTAR) site, and has been a part of a broader USDA-ARS watershed network since 1971; thus, it has ample, high quality data [40]. The Soil Water Assessment Tool (SWAT) [41] has been used to model this watershed for many years in conjunction with this LTAR research program. The SWAT model [41] is a widely used watershed-scale process-based distributed parameter hydrological simulation model designed to simulate hydrological components, and enables detailed representation of agricultural land cover and management practices [42,43]. SWAT includes algorithms for predicting how P, T, and CO₂ concentration affect plant growth and hydrological components.

The GCEW is a headwater watershed in the Salt River Basin, and the geophysical context of the study area represents the Central Claypan region, which includes Northeast Missouri, Southeast Iowa and Southern Illinois [40,44]. GCEW includes a restrictive clay layer in the subsurface soil (B-Horizon), which results in lower hydraulic conductivity, and ultimately, in considerable surface water runoff despite shallow slopes [45,46]. The presence of this restrictive clay layer may result in more significant responses in hydrological components as a result of future changes in T and P. For example, intense precipitation events could result in higher runoff, whereas short-term dry periods may result in short-term drought, both of which can impact agricultural productivity. Hydrological models with a coarse spatial resolution may not be able to adequately capture hydrological responses to climate in areas with unique physical features, such as the claypan in the GCEW.

The objectives of this study are to contribute to the body of knowledge developed in conjunction with the LTAR project in the GCEW, to characterize the potential hydrological impacts in relation to climate change in the GCEW, and to compare hydrological output from models of different spatial resolution, as well as with and without localized downscaling of weather data. This is an important inquiry for developing modeling techniques to assess risks to hydrological components and agricultural activity. To support these objectives, 12 climate model datasets obtained from CMIP5 were statistically downscaled using historic T and P data from GCEW [47] and modeled hydrological responses were analyzed using an ensemble approach. To support the last objective, SWAT simulated runoff was compared with simulated runoff from the Lund–Potsdam–Jena managed Land (LPJmL) and Jena Diversity-Dynamic Global Vegetation (JeDi-DGVM) models forced with precipitation data bias corrected using a nonlinear regression with a spatial resolution of 0.5 degree [48]. These data were made available via the Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP) [49]. Runoff results were compared to evaluate the performance of these varying spatial resolutions, to generate more insight on why there is a need for detailed simulation with higher spatial resolution in an agricultural setup.

2. Materials and Methods

2.1. Watershed Description and Data Sources

The study was conducted in the GCEW, located in Boone and Audrain counties, Missouri (Figure 1), with a drainage area of approximately 73 km² [40]. The land uses in the watershed include row crops (72.8% of total land area), pasture (14.3%), forest (6.0%), and urban land (6.9%) (

Table 1. Land use land cover for the Goodwater Creek Experimental Watershed based on NASS 2010.

Land Use	Area (ha)	%Watershed Area
Corn	1702.6	23.1
Forest	444.5	6.0
Hay	540.7	7.3
Pasture	1049.5	14.3
Soybean	3052.1	41.4
Switchgrass	73.5	1.0
Urban	501.5	6.9

). The study area received an average annual precipitation of 1027 mm between 1993 and 2010, with maximum and minimum annual precipitation of 1614 and 776 mm for 2008 and 2007 respectively [47]. The historical average annual estimated snowfall depth for the study watershed is around 380 mm; given the difficulty in obtaining snow measurements, snow water equivalent is measured in the study area [47]. The daily average maximum and minimum temperatures during this same period were 17.5 °C and 6.5 °C respectively. The average annual discharge of the watershed (expressed as a depth over the drainage area) is 310 mm, proportioned of surface (80%) and base flow (20%) [50].The watershed has flat topography, with 69% of the area having a slope of 0–2%, 15.2% having a slope of 2–3%, and the remaining 15.8% having a slope greater than 3%. The average elevation of the watershed is 256 m above sea level, ranging from 223 to 281 m.

Data required for the SWAT model included a Digital Elevation Model (DEM), soil database, Land Use/Land Cover (LULC), and weather datasets. The DEM (10-m resolution), soil (30-m resolution), and land use/land cover (30-m resolution) datasets were downloaded from Missouri Spatial Data Information Service (MSDIS) (http://msdis.missouri.edu/data/dem/) [51], Soil Survey Geographic Data (SSURGO) (http://websoilsurvey.sc.egov.usda.gov/App/WebSoilSurvey.aspx) data set [52], and National Agricultural Statistics Service (NASS) Cropland Data Layer (CDL) (http://nassgeodata.gmu.edu/CropScape/) of the United States Department of Agriculture (USDA) [53], respectively. The historical weather datasets required for the model, including daily P, daily maximum and minimum air T, wind speed, and solar radiation, were obtained from the weather station and rain gauge located in the watershed (Figure 1). The precipitation datasets from 5 different rain gauges (Figure 1) were used to better represent the spatial variability of precipitation across the watershed.

2.2. Description of SWAT Model and Model Setup for GCEW

The SWAT 2012 version 635 was used by incorporating the changes included in an earlier modeling effort to simulate the percolation through and saturation above the claypan [44] which controls the hydrological components of Central Claypan region. The watershed was delineated using Arc Hydro tools built in the Arc SWAT interface using DEM data for the area (MSDIS). The watershed was delineated using the threshold drainage area of 200 ha, with delineation resulting in 7 sub-basins. The delineation scheme was consistent with previous modeling efforts for GCEW [44]. LULC data were assumed to be constant throughout the simulation. A two-year rotation of corn and soybean (conventional tillage corn and no-till soybean rotation) was used to represent the dominant year-to-year cropping patterns in the watershed. Pasture land was further divided into pasture1, pasture2, hay and switchgrass, to include more detail for SWAT simulation and to represent proper stocking density by applying grazing in alternate pasture [44]. The watershed was divided using five slope classes (<0.5%, 0.5–1%, 1–2%, 2–3%, and >3%). The overlay of soil, land use layer, and slope, resulted in 93 hydrological response units (HRUs). Subbasin area thresholds were used to define minor land uses, minor soils for each land use, and minor slope ranges for each land use/soil combination, using values of 5% for land use, 20% for soils, and 25% for slopes. Minor land uses, soils, or slope ranges were not simulated and the corresponding area was redistributed among the simulated ones. The historic observed data values for all five weather variables (T, P, wind, relative humidity, and solar radiation) were input for historic model calibration and validation. The wind, relative humidity and solar radiation were generated using the SWAT weather generator for future simulations, and T and P were based on statistically downscaled projections from CMIP5. The Natural Resources Conservation Service (NRCS) curve number (CN) method was used to calculate the surface runoff, and potential evapotranspiration (PET) was estimated using the Penman-Monteith (PM) method. Note that the PM method must be used for the climate change scenarios that account for CO₂ changes. This method accounts for the impact of CO₂ change on plant growth.

2.3. LPJmL and JeDi Model Overview

LPJmL and JeDi-DGVMare dynamic global vegetation models represent vegetation dynamics in addition to land surface processes with a spatial resolution of 0.5 degree [54,55]. Runoff simulation for both models was conducted at a daily time step using daily meteorological variables including P, daily mean T, long wave net radiation, and long wave downwelling radiation gridded at 0.5 degree. The change in leaf area index and stomatal conductance due to changes in CO₂ concentration was applied in the simulation to represent dynamic vegetation growth in both LPJmL and JeDi-DGVM. The saturation excess runoff scheme was used for simulation of runoff, i.e., excess water above field capacity was considered runoff and infiltration processes were not simulated. Evapotranspiration was calculated using the Priestley-Taylor (PT) method. Soil data for the LPJmL model setup was taken from the Harmonized World Soil Database (HWSD). Soil data of HWSD is a 1-km-resolution raster database. HWSD soil data are available for two layers, top- and subsoil [56]. The land use land cover database from gridded land use (HYDE 3.2) [57] with a horizontal resolution of 0.5-km wasused for LPJmL, which means land use cannot be subdivided at a resolution of less than 0.25 km². The JeDi-DGVM model uses a parameter-based land surface model to simulate soil and land use. The parameters used in the JeDi-DGVM model include maximum plant available water storage in the root zone, transpiration rate, leaf area index, fractional vegetative cover, fractional forest cover, snow-free surface albedo, and canopy storage. Details about the parameter-based land surface model used in JeDi-DGVM can be found in Diimenil, et al. [58]. The major differences between LPJmL and JeDi-DGVM include the definition of land use land cover, soil, and dynamic vegetation representation. Both models represent vegetation dynamics using semi-empirical plant functional types (PFTs). Vegetation in LPJmL is represented using only 12 PFTs, compared to the large number of PFTs generated based on 15 trait parameters in JeDi-DGVM [55]. LPJmL and JeDi-DGVM surface runoff output, generated based on weather data from the CMIP5 climate data, were retrieved via the ISI-MIP data repository [https://esg.pik-potsdam.de/search/isimip-ft/] [59] for the grid cell overlapping the GCEW for comparison to SWAT output. The monthly average runoff projections for near and far future for LPJmL, JeDi-DGVM, and SWAT were compared using the paired-t test. Two Climate models, MIROC-ESM-CHEM and IPSL-CM5A-LR, were used for this study to represent both kinds of plausible futures (higher temperature and precipitation change and smaller temperature and precipitation change). Note that the runoff simulation data are available using two impact models (LPJmL and JeDi-DGVM) in the ISI-MIP project.

2.4. SWAT Model Calibration, Validation and Evaluation Criteria

The model was calibrated for daily streamflow during the period 1993–2001 using data collected at the watershed outlet (Weir 1), and the model was validated using 2002–2010 streamflow data at Weir 1 and 1993–2001 streamflow data at Weir 11 [50] (Figure 1). An initial three–year (1990–1992) warm-up period was used to initialize the SWAT model. The model sensitive parameters for streamflow were determined using the SWAT-CUP (SUFI-2) [60], and calibration was conducted manually by adjusting a single parameter at a time to fit the observed and simulated streamflow. Manual calibration was preferred over auto-calibration, as SWAT parameter ranges were well understood based on previous studies [44]. The model was calibrated to get the best fit value of the streamflow-sensitive parameters, and their values were changed within the range defined by Neitsch, et al. [61]. Calibration was conducted through visual assessment of superposition of observed and simulated flow duration curves for whole range of flow, and further model performance statistics were calculated to confirm. The model performance was tested using coefficient of determination (r²) and Nash Sutcliffe Efficiency (NSE). The r² value above 0.60 and NSE values above 0.50 and p-bias of ±10% during the calibration and validation periods were set as a criterion for satisfactory model performance during the model simulation of daily streamflow as defined by Moriasi, et al. [62]. The calibrated SWAT model was used for the future projection of hydrological condition using the statistically downscaled temperature and precipitation data for the near (2016–2045) and far future (2046–2075) periods. The simplified flow chart shown in Figure 2 describes the steps followed in this study.

2.5. Climate Data and Bias Correction Methods

Climate projection data were obtained from the CMIP5 project [63], which provides Bias Corrected Constructed Analog (BCCA) downscaled to one-eighth degree (12 km²) from many GCM models (available at: http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/) [64]. Data from twelve GCM models for four emission scenarios (RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5) for a total of 48 unique datasets were selected for this study (

Table 2. List of the General Circulation Models used in this study.

Model	CMIP5 Climate Modeling Group	Model ID	Reference
1	Beijing Climate Center, China Meteorological Administration	bcc-csm1-1	[65]
2	National Center for Atmospheric Research	ccsm4.1	[66]
3	National Center for Atmospheric Research	ccsm4.2	[66]
4	NOAA Geophysical Fluid Dynamics Laboratory	gfdl-esm2g	[67]
5	NOAA Geophysical Fluid Dynamics Laboratory	gfdl-esm2m	[68]
6	Institut Pierre-Simon Laplace	ipsl-cm5a-lr	[68]
7	Institut Pierre-Simon Laplace	ipsl-cm5a-mr	[68]
8	Japan Agency for Marine-Earth Science and Technology, Atmosphere and Ocean Research Institute (The University of Tokyo), and National Institute for Environmental Studies	miroc-esm	[69]
9	Japan Agency for Marine-Earth Science and Technology, Atmosphere and Ocean Research Institute (The University of Tokyo), and National Institute for Environmental Studies	miroc-esm-chem	[69]
10	Atmosphere and Ocean Research Institute (The University of Tokyo), National Institute for Environmental Studies, and Japan Agency for Marine-Earth Science and Technology	miroc5	[70]
11	Meteorological Research Institute	mri-cgcm3	[71]
12	Norwegian Climate centre	noresm1-m	[72]

[63] and Table S1). Daily P and air T (min and max) were extracted for the four grid cells spanning the GCEW. These data were further downscaled using observed P and T data over the years 1971–2010 from the weather station located in the watershed. The bias in temperature was corrected using the delta method, and the bias in precipitation datasets was corrected using modified quantile mapping.

2.5.1. Statistical Downscaling: Temperature-Delta Method

The delta method aims to match the monthly mean of the corrected temperature values with observed temperature [73]. The correction of the temperature involves linear shifting of future temperatures using a correction factor calculated based on the difference between the historic GCM simulated and historic observed temperature data for each month [29]. Temperature data extracted from CMIP5 were corrected to correspond to the study area using the delta method. The monthly correction factor was calculated by subtracting average monthly baseline CMIP5 simulated T data and baseline observed T data (1971–2010). Daily T over the analysis period was corrected with an additive term (b_j) computed for each month as Equation (1):

$$b_j=\frac{{\displaystyle \sum _{C=1}^{C_i}{T}_{ij}-{\displaystyle \sum _{C=1}^{C_i}\widehat{T_{ij}}}}}{C_{ij}}$$

(1)

where, C_i = is the total number of days for month j over the entire baseline historic period, $T_{ij}$ is the daily observed T in °C (maximum or minimum), and $\widehat{T_{ij}}$ is the daily model temperature (maximum or minimum) on day i in month j. Finally, the additive factor for each respective month was added to raw data for bias correction (Equation (2)).

$$T_{ij}=\widehat{T_{ij}}+b_j$$

(2)

2.5.2. Statistical Downscaling: Precipitation-Quantile Mapping

Quantile mapping adjusts the daily CMIP5 precipitation values to match the statistical distribution of the observed precipitation. This method is known by many other names, for example as ‘distribution mapping’ [74,75] and ‘probability mapping’ [76,77]. One major assumption of quantile mapping is that the precipitation distribution does not change over time, and variance and skewness of the distribution remain the same with only a change in mean [78]. Cumulative Distribution Functions (CDFs) were constructed for both CMIP5 and observed daily P (1970–2010) on a monthly basis. A transfer function was constructed to transform CMIP5 data values to probabilities based on the CDF of the model distribution, and these values were transformed back to data values using the inverse CDF or quantile function of the observed distribution. This procedure is expressed mathematically as (Equation (3)):

$$P_{i,j}^{corrgcm}=CD{F}_{obs,j}^{-1}\left[CD{F}_{gcm,j}^{his}\left(P_{i,j}^{raw}\right)\right](P_{i,j}^{raw}\le P_{max,j}^{hisgcm})$$

(3)

where, $P_{i,j}^{raw}$, $P_{i,j}^{corrgcm}$ and $P_{max,j}^{hisgcm}$ are the projection’s raw CMIP5, corrected CMIP5, and maximum historical CMIP5 precipitation for i-th day of j-th month, respectively. $\mathrm{The}CD{F}_{obs,j}^{-1}$ is the quantile function of the observed precipitation and $CD{F}_{gcm,j}^{his}$ is the CDF of an individual GCM.

Quantile mapping was modified to correct two major limitations in CMIP5 P data: (1) the P frequency or dry bias, i.e., the data contain a large number of very low P days, i.e., “drizzle days”, [79] and (2) an underestimation of extreme P compared to historical observation. The number of P days was more than 300 days/year for projected gridded data compared to 94 days/year in the observed dataset (1971–2010). To correct for these excess “drizzle days”, a P threshold was determined and the P values below the threshold were removed (Equation (4)). A P threshold for each month was determined for each GCM model to ensure the same number of P days for GCM baseline and observed, as was done by Grillakis, et al. [80]. Monthly thresholds ranged from 0.60 to 4.5 mm/day; similar to ranges of P thresholds reported by Themeßl, et al. [81]. Quantile mapping was performed after removing “drizzle day” corrections.

$$P_{i,j}^{raw}\left(d\right)=\{\begin{array}{c}0,if{P}_{i,j}^{raw}<P_{threshold,j}\\ P_{i,j}^{raw},otherwise\end{array}$$

(4)

The second limitation of quantile mapping is that the method limits the upper value of future daily precipitation to the highest value observed in the past. Given that climate projections indicate that precipitation events are likely to increase in magnitude and this sort of extreme event poses potential risks a second adjustment was made to the CMIP5 P data to allow a daily P to exceed the historic maximum using Equations (5) and (6). P was scaled on a monthly basis to represent the monthly extreme. A separate scaling factor was calculated for each month, where S_j is the scaling factor for a given month and $P_{i,j}^{hisgcm}$ is the CMIP5 daily P data for month j during the historic period. This approach retained the relationship between the historic period and future periods within a CMIP5 P dataset for a given model without constraining future P data to the maximum historical observation.

$$P_{i,j}^{corrgcm}=S_j\ast P_{i,j}^{raw}(P_{i,j}^{raw}>P_{i,j}^{hisgcm})$$

(5)

$$S_j=\frac{Max\left(P_{i,j}^{obs}\right)}{Max\left(P_{i,j}^{hisgcm}\right)}$$

(6)

2.6. Simulation Scenarios. Impacts of Climate Change on Water Yield, Evapotranspiration, and Surface Runoff

The validated SWAT model was run at a daily time step starting in 1981 and continuing until 2075. For analysis, the baseline represents 1981–2010, the near future is defined as 2016–2045, and the far future is defined as 2046–2075. Different atmospheric CO₂ concentrations were applied within SWAT for each RCP for different time periods; these values were calculated by averaging the projected CO₂ concentrations over each 30-year period. The average was used because it is not possible to change CO₂ concentration annually in SWAT. The details of CO₂ concentration used for each model for each time frame and RCP are presented in

Table 3. Variable initial carbon dioxide concentration used for SWAT simulation under different emission scenarios [84].

Scenarios	2016–2045	2046–2075
	CO2 Concentration in Part per Million (ppm)
RCP 2.6	428	440
RCP 4.5	437	507
RCP 6.0	431	515
RCP 8.5	454	611

. The CO₂ concentration of 330 ppm was used for all the historic simulations. The ensemble median and quartile values (Q1 and Q3) of monthly water yield, evapotranspiration, and surface runoff were calculated for the three analyses periods. Reporting results as ensembles allows for greater confidence in simulated hydrological outputs [82]. Pierce, et al. [83] found that the average of a multi-model ensemble is a better predictor than any individual GCM for global studies examining the mean climate, as the ensemble of results cancels out the bias and errors.

The modified PM method was used to calculate ET. This equation allows for the incorporation of variability in radiation-use efficiency, plant growth, and transpiration induced by change in atmospheric CO₂ concentration [43]. SWAT simulates the impact of change in CO₂ concentration on stomatal conductance by linear interpolation between current conditions and a 40% decrease of stomatal conductance for doubling of CO₂ up to 660 ppm [43]. However, SWAT does not simulate effects of CO₂ on leaf area index; thus, this study does not incorporate the impact of increased CO₂ concentration on leaf area.

3. Results and Discussion

3.1. Model Calibration and Validation

The default and adjusted model parameter values are presented in

Table 4. Sensitivity results of key SWAT parameters for stream discharge in the Goodwater Creek Experimental Watershed including the default values.

Parameter	Definition	Default Value	Adjusted Value
ESCO	Soil evaporation compensation factor	0.95	0.88
GWQMN	Shallow aquifer depth of water required for return	0	55
	flow to occur (mm)
GW_DELAY	Groundwater delay time (days)	31	55
SOL_AWC	Available water capacity of the soil layer (mm H2O/mm soil)	-	0.04 †
CH_N	Manning’s “n” value for the main channel	0.014	0.019
CH_K	Effective hydraulic conductivity in main channel	0	0.08
	Alluvium (mm/hr.)
SMTMP	Snow melt temperature (°C)	0.5	−2.5
SMFMN	Snow melt min rate (mm H2O/°C-day)	4.5	1.5
EVRCH	Reach evaporation adjustment factor	1	0.5
MSK_X	Weighting factor that controls the relative importance	0.2	0.1
	of inflow and outflow in determining the storage in the reach
MSK_CO2	Calibration coefficient used to control impact of the	0.25	3.5
	storage time constant for low flow upon the time constant value calculated for the reach
ALPHA_BF	Baseflow alpha factor (1/days)	0.0048	0.1
SHALLST	Initial depth of water in the shallow aquifer (mm)	1000	600

Note: † denote the relative change in the parameter value.

. The model parameters controlling surface runoff, groundwater, and snowmelt were important parameters considered during calibration. The calibrated model parameters are similar to those reported by previous researchers for GCEW [44]. Daily streamflow simulation at weir1 (watershed outlet) resulted in r² of 0.66, NSE of 0.66, and p-bias of 5.6% during calibration (1993–2001), and r² of 0.60, NSE of 0.59, and p-bias of −6.9% during the validation. The model was further validated using the streamflow data at weir11 resulting in r² of 0.73, NSE of 0.68, and p-bias of 18%. SWAT simulated dry (2007) and wet (2008) years accurately, indicating that the model was suitable for evaluating climate change impacts. The flow duration curve is a description of the frequency at which a given discharge is reached or exceeded over the simulation period. The flow duration curves for observed and simulated flow for calibration and validation periods are presented in Figure 3. The calibration objective was to superimpose observed and simulated flow duration curves to match low, intermediate, and high flow.

3.2. Downscaled GCM Projection

3.2.1. Temperature

Minimum temperatures were under-predicted and maximum temperatures, particularly in the summer months, were over-predicted in CMIP5 output, compared to historical observations. The scaling factors used for correction of the maximum and minimum daily T are presented in supplement Tables S2 and S3. The maximum bias was for the month of March for both the maximum and minimum T. Average monthly T (observed, CMIP5, bias-corrected) for the period 1981–2010 for one of the GCM models included in the study is presented in Figure 4 as an example. These results are consistent with previous reports of over-prediction of annual mean T for CMIP5 GCMs for Northern Eurasia [85,86], Northern hemisphere [87], and globally [88].

Figure 5 shows GCM projected changes in average annual T, between the baseline (1981–2010) and near (2016–2045) (Figure 5a) and far future (2046–2075) (Figure 5b) for the 12 sets of downscaled T data from CMIP5 for four RCP scenarios. Ensemble average T change in the near future period shows minimum T change for RCP 6.0 (1.13 °C) and maximum T change for RCP 8.5 (1.63 °C). The far future period shows minimum T change for RCP 2.6 (0.33 °C) and maximum T change for RCP 8.5 (1.92 °C). The most significant increase in T occurred for models under RCP 8.5 scenarios for both near and far future (Figure 5a,b). This is consistent with the assumptions associated with technological and socioeconomic conditions leading to the concentrations of GHGs associated with each RCP. That is, RCP 2.6 includes a decrease in the concentration of GHGs in the atmosphere to below currently observed levels, while RCP 8.5 represents the largest increases of GHG concentrations in the atmosphere. For RCP 4.5, atmospheric CO₂ peaks around 2040, then declines compared to RCP 6.0, where CO₂ peaks in 2080 and declines [89]. Therefore, in the near future, RCP 6.0 has the lowest CO₂ concentration compared to RCP 4.5 scenarios resulting in the least T change (Figure 5a). The spread of predicted changes in T for models in RCP 4.5 (0.7 °C–2.2 °C) and RCP 8.5 (0.9 °C–2.5 °C) demonstrates the value of a multiple model ensemble approach for future climate change impact assessment studies.

3.2.2. Precipitation

Figure 5 also shows GCM projected percentage changes in average annual P, for baseline (1981–2010), near future (2016–2045), and far future (2046–2075) respectively, for the 12 different GCMs used in this study. An ensemble P change for the far future scenario was the least for the RCP 6.0 scenario (1.0% increase) and most for RCP 8.5 (6.9% increase). The highest simulated increase in average annual precipitation was for RCP 8.5 with a change in P ranging from −7% to 16% across 12 GCMs in the near future and −10% to 28% in the far future, followed by the RCP 4.5. Note that in near future (2016–2045), RCP 4.5 had a higher CO₂ concentration than RCP 6.0 [90], which may explain the elevated T and P for RCP 4.5 compared to RCP 6.0 over this period. Use of variable monthly thresholds for daily precipitation enabled exact matching of wet days between historic modeled P to observed P, which is critical for quantile mapping. The monthly P thresholds and scaling factors used to correct these biases are presented in supplement Tables S4 and S5. Quantile mapping of the CMIP5 data to the historical data for one model under the RCP 8.5 scenario (micro-esm.1) is presented in Figure 6.

Monthly ensemble median P results indicate increased spring P compared to the baseline for both near and far future for all the emission scenarios (Figure 7). A previous study also reported an increase in precipitation during the spring months based on historic data for Central US [91]. Decrease in ensemble monthly P compared to the baseline was simulated for near-future for the months from July-November for RCP 2.6 (Figure 7). For RCP 4.5, the decrease in P for July was 7% and 14%, and August was 6% and 14%, for near and far future respectively. The decrease in precipitation for July and August was 12% and 14%, and 9% and 21% for far future of RCP 6.0 and RCP 8.5 respectively. There was minimal change in the monthly precipitation during the winter months. Occurrence of more P during spring may delay agricultural operations, and dry summers with very low P may result in short-term drought and possible crop losses. Previous regional modeling efforts using the Regional Climate Models (RCM) to assess the climate change impact on hydrology of the Mississippi River Basin indicated increased future P, as high as 21% on a mean annual basis [28]. Downscaling of P based on point scale (weather station) data is more relevant at the small watershed scale and may help to better understand how changes in P and T impact the hydrological components in the watershed.

3.3. Climate Change Impact on Hydrological Output

3.3.1. Water Yield

The ensemble median annual water yield was predicted to increase in both the near and far future compared to the baseline for RCP 2.6, RCP 4.5, and RCP 8.5 climate scenarios. Climate scenario RCP 6.0 showed 1% reduction in median water yield for far future scenario, due to minimum P change projected for far future for models under these climate scenarios. The highest change of 29% was observed for the far future scenario of RCP 8.5, which also resulted in the greatest increase for the near future (19%). The changes in water yield for RCP 2.6 and 4.5 were similar for near future, with the water yield increasing by around 10%. The water yield for RCP 2.6 and RCP 4.5 increased by 12% and 6% respectively, for the far future (

Table 5. Ensemble median and quartiles of annual simulated water yield, surface runoff and ET changes (in percent) for the Goodwater Creek Experimental Watershed.

Water Yield (% Change over the Baseline)
	RCP 2.6		RCP 4.5		RCP 6.0		RCP 8.5
	Near Future	Far Future	Near Future	Far Future	Near Future	Far Future	Near Future	Far Future
Median	9.9	12.0	10.4	6.2	5.4	−1.0	19.0	29.4
1st Quartile	5.6	8.4	2.4	0.1	−2.9	−2.7	7.7	19.9
3rd Quartile	11.1	19.2	11.6	28.3	8.9	16.3	23.1	31.4
Surface Runoff
Median	10.1	12.2	9.7	7.2	5.9	−0.9	20.9	29.9
1st Quartile	6.6	10.0	3.6	1.6	−3.0	−3.0	7.8	20.0
3rd Quartile	12.4	20.0	11.3	28.1	8.2	17.6	24.4	31.0
Evapotranspiration
Median	0.1	3.1	1.6	0.1	0.8	−0.3	2.1	−2.1
1st Quartile	−1.1	1.9	1.8	−0.8	−0.7	−2.5	0.4	−1.4
3rd Quartile	1.0	2.2	0.1	1.3	0.5	0.7	0.4	−2.4

). The models in RCP 6.0 showed very small changes in P (2% for near and 1% for far future). The increase in water yield was due to increased P; the model with the highest increase in P resulted in the highest increase in water yield. For example, the greatest increase in water yield occurred for the month of May where the P was maximum for RCP 8.5 (Figure 7 and Figure 8).

The first and third quartiles for water yield in Figure 8 show the variation in predicted water yield, which reflects the differences in GCM projections for each RCP. Despite this variation across GCM outputs, in all cases the ensemble median for near and far future was above the baseline median (Figure 8, Table 5).

For all the RCPs, the GCM ensemble range indicated a projected increase in the 3rd quartile of water yield. Monthly peak water yield occurs in June, except for the far future of RCP 6.0 and RCP 8.5 scenario where it occurs in May (Figure 8). Elevated water yield during the spring months (March, April, and May) is due to an increase in spring P projected by all four emission pathways. Such change in distribution of water yield compared to baseline may cause increase risk of extreme (flood and drought) events in the watershed. Similar trends were simulated for the region, as shown by elevated future spring and mean annual water yield for the Missouri River Basin simulated using the GCM output and the SWAT model [27]. These results, as well as previous studies [92], document that the runoff processes can be amplified due to fluctuations in P, resulting in significant impacts on streamflow. Our study shows similar results: the highest increase of water yield (83%) (270.5 to 494.9 mm) was simulated using the microc-esm.1 model under RCP 8.5, with a 28% (from 969.2 to 1239 mm) increase in precipitation compared to the baseline.

3.3.2. Surface Runoff

Ensemble mean annual surface runoff was predicted to increase compared to the historic period in all scenarios with the exception of RCP 6.0 in the far future (Figure 9). RCP 8.5 resulted in the greatest increase compared to the historic period in surface runoff for both near future (20.9%) and far future (29.9%) (Table 5). The first and the third quartile ranges for surface runoff in Figure 9 show the range in predicted runoff associated with individual climate models. These variations reflect the difference in GCM projections among RCPs (Figure 9 and Table 5). The largest increases are shown in the spring months as a direct result of increased precipitation. The monthly peak runoff occurred in June, except for in the far future of RCP 8.5, where it occurred in May (Figure 9). The magnitude of the peaks is larger than those observed in the historical dataset for all the RCPs, except RCP 6.0, which shows small change in ensemble water yield (Table 5). Decreased surface runoff in August was simulated for RCP 4.5, RCP 6.0, and RCP 8.5, compared to historic baseline (Figure 9).

The relationship between increased P and surface runoff is a function of physical characteristics of the location, e.g., soil type, land cover, and in a mathematical model is also a function of the spatial and temporal resolution of the model inputs and response. Hydrological models use geomorphologic and topographic parameters to characterize land forms to represent hydrological process [93]. The parameters used to represent the scale lead to differences in the simulated results. For example, SWAT uses detailed parameters to represent small sub-areas within a watershed, in this case 93 HRUs to represent the heterogeneity in GCEW, compared to single, spatially-averaged parameters applied in JeDi and LPJml to represent the grid that encompasses GCEW.

In this analysis, the GCM model microc-esm.1 under RCP 8.5 results in the largest increase in P, with an average annual increase of 28% (969.2 to 1239 mm), results an increase in runoff of 88% (230.6 to 434.9 mm). Jha, et al. [28] found an increased runoff of up to 51% in the Upper Mississippi River Basin (UMRB) for a 21% increase in P, using regional climate model output and SWAT. The SWAT model version used by Jha, et al. [28] was not modified to represent percolation through claypan soil. In contrast, a modified version of SWAT was used in our study to represent claypan hydrology [44]. The GCEW’s relatively larger modeled surface runoff response is likely due to increased heavy precipitation events, differences in the model spatial resolution, and accounting for the claypan condition in the watershed [44]. Lower saturated hydraulic conductivity of the claypan layer results in lower infiltration and higher runoff response. The representation of small-scale details (e.g., claypan) conditions is often overlooked in regional-scale simulations, as these models are constructed to simulate average conditions. Model correction to represent each unique local situation is not feasible in a model with a large spatial resolution. The small watershed-scale studies are important, as they help to better represent heterogeneity within the watershed, which can ultimately help local-scale decision-making.

Monthly surface runoff estimated using SWAT in GCEW was compared to the surface runoff output from the LPJmL and JeDi-DGVM models forced using two sets of CMIP5 data (miroc-esm-chem and ipsl-cm5a-lr), and obtained from ISI-MIP, Figure 10 and supplement Figures S6 and S7 [94]. For peak runoff, the flow duration curve comparison showed good agreement between SWAT and observed data (Figure 3). A lower peak value by the gridded models indicates an under-prediction of peak flow. The comparison indicates that the gridded models were unable to simulate the monthly peaks predicted by the SWAT simulations (Figure 10). Average annual surface runoff during 2016–2075 was over-predicted by 6–8% using LPJmL and under-predicted by 5–30% using JeDi-DGVM compared to the SWAT model. Paired t-test results indicate significantly lower prediction of average monthly runoff by JeDi-DVGM compared to SWAT and greater prediction by LPJmL. The results for paired t-tests are presented in

Table 6. Comparison of average monthly runoff predicted by the hydrological models SWAT, JeDi-DVGM, and LPjmL forced with two different global climate model (GCM) datasets.

GCM	Runoff (mm)-SWAT		Runoff (mm)-LPJmL		Runoff (mm)-JeDi-DVGM
Period	NF	FF	NF	FF	NF	FF
ipsl-cm5a-lr	24.4 aA	20 cC	23.9 a	23.3 d	26.5 A	15.5 D
miroc-esm-chem	20.7 aA	26.4 cC	24 a	26.8 c	18.4 B	14.8 D

Note: Means followed by different letters between the columns for each of the time periods (Near Future (NF) and Far Future (FF)) are significant at p ≤ 0.05. Small letter shows difference in runoff between SWAT vs. LPJmL and capital letter shows difference between SWAT vs. JEDI-DVGM. GCM: General Circulation Model.

.

The causes for these differences include model assumptions, model parameterization, input spatial resolution, and model structure in the respective models. Evapotranspiration was calculated using the PT method in LPJmL, compared to the PM approach in SWAT. Suleiman and Hoogenboom [95] report that the accuracy of the PM method for ET estimation is superior for predicting ET in the U.S., and that the PT method tends to overestimate ET for U.S. conditions. The difference in runoff prediction between the LPJmL and JeDi-DGVM models is also due to differences in approaches for representing dynamic vegetation in the model as presented in Section 2.3. The lower peak monthly runoff and greater low runoff values in the gridded model output may also be due to simulation at a coarse resolution with a homogenous set of parameters compared to the fine resolution input data of SWAT model (as described in Section 2.1 and Section 2.3), as well as differences in accounting for local weather (i.e., statistical downscaling) and local hydrological conditions (e.g., claypan). The comparison shows the importance of impact assessment at the smaller scale to predict hydrological change through representing soil and land use with greater detail. Recommending simulation at local scale for every single point on the globe can be resource intensive and even impossible. However, creating simulations of small watersheds that represent major characteristics of regions can be useful for anticipating risks and corresponding decision-making, especially within the agriculture sector, which can experience quick response to change in precipitation and temperature. Likewise, global gridded models, which are generally better at simulating average conditions, may be sufficient for other applications. For a global resource analyst, general trends noted by a global gridded model may be sufficient. Although LPJmL was able to represent the long-term average monthly conditions, it failed to represent extremes (peak and minimum runoff), which are important from an agricultural perspective.

3.3.3. Evapotranspiration

Median annual ET in the near and far future had little change, −1% to 3% compared to historical observations on an annual basis (Table 5). There was minimal change in annual ET in the near future for RCP 2.6 and RCP 6.0 (Table 5). However, there was an increase in ET in spring months and a decrease in summer months for both the near and far future for all four RCPs (Figure 11). The decrease in summer ET was maximum for RCP 8.5 and minimum for RCP 2.6.

Results indicated shifts in the peak monthly ET from July (baseline) to June for RCP 4.5, RCP 6.0, and RCP 8.5 in the far future (Figure 11). Land use and land cover were fixed during the entire simulation period. Therefore, variation of the simulated ET from the watershed was mostly controlled by meteorological variables and corresponding plant growth responses. For each crop, a base T must be reached for plant growth to begin in SWAT. With increased T, base T values were reached sooner than they have been reached historically, resulting in a shift in the plant growth cycle to earlier in the calendar year. This earlier emergence of biomass resulted in a corresponding shift in evapotranspiration to earlier in the season, with an earlier peak for all far future RCPs except for 2.6. The impact of elevated T on the hydrological cycle due to shift of the crop growth cycle has been documented previously by Ficklin, et al. [96]. However, it remains unclear whether farmers would actually plant crops earlier. In addition, for this study region, increased spring precipitation may preclude early planting.

4. Conclusions

A calibrated and validated SWAT model for the GCEW was used to simulate the effect of climate change using downscaled T and P data from 12 GCMs under four emission scenarios. Downscaled climate model projections suggested an increase in spring precipitation and increased temperatures in the future for the study region, with the magnitude varying with the GCM model and emission scenarios. The change in precipitation for the extreme scenario (RCP 8.5) ranged from −7% to 16% for near future, and from −10% to 28% for far future scenarios. The majority of CMIP5 data showed a reduction of summer precipitation and increase in spring precipitation. SWAT simulated hydrological results indicated a range of possible futures due to increased temperature, precipitation changes, and elevated CO₂. Results show increased water yield, surface runoff during spring months, and a shift in ET for all the RCPs except RCP 2.6 in the far future. The greatest increase in the median water yield and surface runoff (29% and 30%) compared to the baseline was for the far future of the most extreme future climate scenario, i.e., RCP 8.5. Small changes in precipitation amounts can impact the runoff and water yield response, which is especially apparent at the small watershed scale. The shift in the peak ET from July to June for RCP 4.5, RCP 6.0, and RCP 8.5 indicated the probable impact of increased CO₂ concentration and temperature on planting dates and seasonality, which may ultimately impact future crop yields. Of particular concern for this watershed is the increase in precipitation, and the potential for more frequent occurrence of high-intensity rainfall events in the springtime. This could negate potential for earlier planting, noted by some researchers as a potential benefit of increasing temperatures due to climate change. Further, this study has found that cumulative precipitation is likely to decrease in the late summer; this, combined with higher temperatures, could lead to drought conditions in July and August. Depending on planting date, this could overlap with the tasseling stage of corn development. Modeling at this scale allows for a finer understanding of the patterns and processes of the hydrological response, as well as of the potential implications with regard to agricultural activity, to changes in temperature and precipitation.

The comparison with gridded model output indicated that hydrological modeling at larger scales fails to capture peak and minimum runoff, which may be very important for adaptation responses at local scales. Average annual surface runoff during 2016–2075 was over-predicted by 6–8% using LPJmL and under-predicted by 5–30% using JeDi-DGVM, compared to the SWAT model. Global gridded data with coarse spatial resolution does not capture small-scale details that are necessary for agricultural management decisions. Simulation using higher spatial resolution, as well as additional downscaling of weather data provided by CMIP5, helps to adequately represent the hydrological components of small watersheds. This may be particularly important in watersheds with problematic soils, such as claypan soils found in GCEW and the surrounding region, which have a very thin top-soil layer above a claypan to hold and supply water to plants.