# Potential Impacts of Projected Climate Change under CMIP5 RCP Scenarios on Streamflow in the Wabash River Basin

Authors: Jingrui Wang; Litang Hu; Didi Li; Meifang Ren
Journal: Advances in Meteorology
 (2020)

---

## Abstract

Global climate change is becoming an increasingly important issue that threatens the imperiled planet. Quantifying the impact of climate change on the streamflow has been an essential task for the proper management of water resources to mitigate this impact. This study aims to evaluate the skill of an artificial neural network (ANN) method in downscaling precipitation, maximum temperature, and minimum temperature and assess the potential impacts of climate change on the streamflow in the Wabash River Basin of the Midwestern United States (U.S.) using the Soil and Water Assessment Tool (SWAT). A statistical downscaling technique based on an ANN method was employed to estimate precipitation and temperature at a higher resolution. The downscaled climate projections from five general circulation models (GCMs) under the three representative concentration pathway (RCP) scenarios (i.e., RCP2.6, RCP4.5, and RCP8.5) for the periods of 2026–2050 and 2075–2099 as well as the historical period were incorporated into the SWAT model to assess the potential impact of climate change on the Wabash River regime. Calibration and validation of the SWAT model indicated the streamflow simulations matched the observed results very well. The ANN method successfully reproduced the observed maximum/minimum temperature and precipitation; however, bias in precipitation was observed in regard to the frequency distribution. Compared with the simulated streamflow in the historical period, the predicted streamflow based on the RCP scenarios showed an obvious decreasing trend, where the annual streamflows will be decreased by 13.00%, 17.59%, and 6.91% in the midcentury periods and 25.29%, 27.61%, and 15.04% in the late-century periods under the RCP2.6, RCP4.5, and RCP8.5 scenarios, respectively. Climate warming dominated the streamflow decrease under the RCP2.6 and RCP4.5 scenarios. By contrast, under RCP8.5, the streamflow was affected by the joint actions of changes in temperature and precipitation.


---

## Body

## 1. Introduction

Global climate change, which is mainly characterized by climate warming, is an unequivocal reality. Multisource data verified that the annual-mean globally averaged surface temperature has already increased approximately 0.89°C during the period of 1900–2012 [1, 2]. The International Panel on Climate Change (IPCC) suggested a greater than 90% probability that global warming is caused by human activities, such as the build-up of heat-trapping gases in the atmosphere because of fossil fuel burning [3]. Recently, scientists explicitly claimed that some extreme weather events are blamed on anthropogenic factors [4]. Fresh water resources are immensely important to human society and ecosystems; however, these resources are vulnerable and have the potential to be detrimentally affected by climate change [5, 6]. Global hydrological cycles are expected to be accelerated under the conditions of climate warming given that the capacity of the air to hold water vapor increases exponentially in a warming climate and precipitation will increase on average, which will reduce the increase of people who live under conditions of water stress [5, 7]. However, the study also demonstrated that climate change is projected to cause shifts in precipitation patterns and a likely increase in the frequency and distribution of floods and droughts. The United Nations World Water Assessment Programme (WWAP) estimated that approximately 30% of the global population is impacted by either flood or drought events [8]. An increasingly tense contradiction between supply and demand of fresh water resources makes it necessary to predict the hydrological response to climate change and to manage the water resource in a sustainable manner.As the state river of Indiana, USA, historically benefitting from the abundant water resource with trends of increasing stream flows, the Wabash River Watershed with a drainage area of 85,237 km2 is a key contributor in corn and soybean production [9]. A recent report evaluated the degree to which the economy of each state in the U.S. was dependent on water resources, and Indiana ranked first [10]. For instance, the severe and devastating drought in 2012 in the Midwestern United States resulted in significant loss in economic efficiency [11]. Moreover, a pronounced alteration in frequency and amount of precipitation and increase in average temperatures have been observed in recent years [12]. An investigation on the impact of climate change on the annual streamflow over the continental U.S. found that a 1% change in precipitation would result in a 1.5–2.5% change in the streamflow in the Wabash River Watershed. Storage processes play a key role in shaping the nonlinear precipitation-runoff relationship [13]. Mishra et al. [14] assessed the drought caused by the climate change during 1916–2007. The results showed that precipitation, minimum temperature, total column soil moisture, and runoff experienced upward trends, while maximum air temperature, frozen soil moisture, and snow water equivalent experienced downward trends. Chien et al. [15] studied the potential impact of climate change on the streamflow over the Midwestern U.S. under CMIP3 SRES scenarios, and a drastic decrease up to 44.6% was projected. Recently, investigators have examined the uncertainty in modeled and observed climate change impacts on the American Midwest hydrology and concluded that the impacts of climate change on the regional hydrology of this region remain uncertain [16]. By the end of this century, the IPCC predicted that the global surface temperature is likely to be 2°C warmer than the temperature experienced during the period of 1850–1900 [2]. Even small perturbations in temperature can result in significant impacts on the mean annual streamflow, and the local impacts of 2°C warming are beyond what many societies can cope with [17, 18]. The Wabash River runoff is highly sensitive to temperature alteration [19]; thus, it is essential to quantify the potential impacts of future climate changes on the hydrological regime over the Wabash River Watershed. This information may assist policy-makers and water managers in adopting strategies.Typically, hydrologic models are combined with climate scenarios generated from global climate models (GCMs) to produce potential scenarios of climate change effects on water resources over a large range of scales, including the global scale, continental scale, large drainage network, regional flow system, and small-scale catchments. The Coupled Model Intercomparison Project Phase 5 (CMIP5) archive that contributed to the IPCC Fifth Assessment Report (AR5) provides an unprecedented opportunity to analyze the projections of the 21st century climate change [20]. However, resolutions of climate projections from the synoptically scaled GCMs are often too coarse (generally between 150 and 250 km) to reflect the regional microclimate and need to be downscaled to the required finer resolution [21]. Basically, two downscaling methods, including dynamical downscaling (regional climate models (RCMs)) and statistical downscaling (SD), are used to bridge this scale mismatch. RCMs are dynamical models that are forced by using meteorological lateral boundary conditions from GCMs. The main drawbacks of RCMs are their complicated design and high computational cost. In addition, biased values, especially for precipitation, are very common in GCMs, which will be passed on to RCMs. During recent decades, various SD methods have been developed and implemented to overcome these limitations [22]. These SD methods fall into three major categories: regression models, weather typing, and weather generator. All of these methods are computationally cheap and relatively easy to implement, while each has its advantages and disadvantages. Among the wide range of regression models, the artificial neural network (ANN) has been one of the most popular models [23, 24]. Harpham and Wilby [25] forecasted the rate of precipitation by means of ANN models, and the results demonstrated the superiority of the ANN in catching the nonlinear relationship between predictors and predictands. Mendes and Marengo [26] compared the ANN with an autocorrelation statistical downscaling model to downscale five coupled GCMs for a region in Amazon Basin and found that the ANN model significantly outperformed the other statistical models for downscaling of daily precipitation. Successful downscaling studies carried out with ANNs under different climate change scenarios can also be found in the literature [27, 28].High-fidelity simulation of climate change is essential for quantitative prediction of its potential hydrological impacts. However, climate modeling is an imperfect science with uncertainties in simulated climate variables that vary in space and with the forecast time horizon. Knutti and Sedláček [29] evaluated the robustness and uncertainties in the CMIP5 climate model projections. Although the models employed in the AR5 exhibit improved model resolution and complexity compared with that in the AR4, epistemic uncertainties are still difficult to narrow. Several solutions in coping with this issue were proposed, such as generating perturbed physics ensembles or improving the estimation of each model’s forced climate signal by using sufficiently large ensembles from single-physics climate model implementations that differ only in the initial conditions. However, both proved to be computationally impractical [30]. A practical method to obtain reliable future climate change information is to utilize state-of-the-art multimodel ensemble climate simulation, which performs better compared with the individual model [31].In the present study, we will design a SWAT hydrological modeling experiment and focus on future changes to the streamflow of the Wabash River Basin resulting from climate change. In view of the good relation of precipitation with the amount of water resources in the Wabash River Basin, evaluation of potential impacts of climate change on the streamflow is of great importance. The main objectives of the study include (1) downscaling precipitation and temperature from five GCMs to a higher resolution to represent the climate pattern of the study area using an ANN method and (2) evaluating the skill of the ANN method in downscaling the daily precipitation, daily maximum temperature, and daily minimum temperature over the study basin as well as quantifying the hydrological response of the Wabash River Basin to the climate change over the midcentury (2021–2050) and the end of this century (2070–2099) using a hydrological model (SWAT model). Minimal work is available to assess the impacts of future climate change on water resources in the Wabash River Basin using the CMIP5 archive and SWAT model. The findings from this study could be used by water managers to develop a comprehensive water resource management plan in this basin.The rest of this paper is organized as follows: Section2 offers a brief description of the study area and various data sources used for hydrological modeling and statistical downscaling. Section 3 describes the statistical downscaling method and hydrologic modeling framework. Section 4 outlines the performance of the ANN method in downscaling the GCMs’ outputs under a statistical framework and hydrological framework, followed by the evaluation results of their use for hydrological runoff simulations. Concluding remarks are presented in Section 5.


## 2. Study Area and Input Data

### 2.1. Study Basin

Originating in West Central Ohio, USA, the Wabash River flows southwest across Indiana to Southern Illinois, where it forms the Illinois-Indiana border before draining into the Ohio River (Figure1). It is the largest northern tributary of the Ohio River [32]. With a drainage area of 85,237 km2, the Wabash River Basin is a 4-digit HUC (hydrologic unit code) area formed by two 6-digit HUC subbasins, viz., the main Wabash River Basin and the Patoka-White River Basin. Developed by the U.S. Geological Survey (USGS) in cooperation with the U.S. Water Resource Council and the U.S. Department of Agriculture Natural Resources Conservation Service (NRCS), the HUC is a numbering system for watersheds to provide a common coding system for state and federal agencies. Every hydrologic unit is identified by a unique HUC consisting of 2 to 12 digits based on the six levels of classification in the hydrologic unit system [33]. Among them, HUC 4 represents the subregion level, delineating large river basins, and HUC 6 subdivides many of the subregions into basins known as accounting units. The Wabash River Basin lies in two physiographic provinces, the Central Lowland and the Interior Low Plateaus, and the altitude varies from 60 to 373 m a.s.l. (Figure 1). Predominant land use in the basin changes from agriculture to urban to forest moving from north to south. Approximately 65% of the total land is row-crop agriculture (predominantly corn and soybeans), while urban land coverage represents 2.3%. The types of soils in the Wabash River Basin include clay glacial till soils in the upper Wabash, loamy glacial tills to moderately thick loess shed, and shale in the lower watershed. The humid continental climate of the Wabash River Watershed is shaped by the presence of the Great Lakes and its location. Situated in the middle of the North American continent and far from the temperature-moderating effects of the oceans, the Wabash River Basin contributes to large seasonal variations in air temperature from hot, humid summers to cold winters. Precipitation in the Wabash River Basin is fairly well distributed throughout the year. The annual precipitation is approximately 1,141 mm in Indianapolis (IN). Over the entire basin, the values range from 940 mm in northern areas to 1,245 mm in southern areas. Long-term average temperatures range from 9.4°C in the north to 13.8°C in the south. The average monthly temperatures range from approximately 0°C in January to 25°C in July [32] (Table 1).Figure 1
Shaded-relief map of the Wabash River Watershed showing the stream channels delineated with ArcSWAT and the locations of the climate data nodes and selected U.S. Geological Survey (USGS) gauge stations. The numbering and exact locations of USGS gauge stations are shown in Table1.Table 1
Four USGS streamflow gauge stations located in the midstream and downstream of the main stream and major tributaries.
ID
Streamflow gauge stations
Catchment area
Latitude
Longitude
Websites
Station 1
Wabash River at Terre Haute
3.15 × 104 km2
39°27′57″N
87°25′11″W
https://waterdata.usgs.gov/in/nwis/uv?site_no=03341500
Station 2
White River at Petersburg
2.88 × 104 km2
38°30′39″N
87°17′22″W
https://waterdata.usgs.gov/nwis/inventory/?site_no=03374000
Station 3
Wabash River at Mount Carmel
7.39 × 104 km2
38°23′54″N
87°45′23″W
https://waterdata.usgs.gov/nwis/inventory/?site_no=03377500
Station 4
Little Wabash River at Carmel, IL
0.83 × 104 km2
38°30′39″N
87°17′22″W
https://waterdata.usgs.gov/nwis/inventory?agency_code=USGS&site_no=03381500


### 2.2. Data Used

#### 2.2.1. Geospatial Data

Topography, land cover, and soil data are used to represent the catchment heterogeneity for the hydrologic analysis and were compiled from different sources (Table2). The USGS National Elevation Dataset (NED) digital elevation model with a 30 m spatial resolution furnished land surface topography for the simulation domain. Land use data were obtained from National Land Cover Database 2011 (NLCD 2011), while the Soil Survey Geographic Database (SSURGO) that incorporated in ArcSWAT and overlapping the study area was used.Table 2
Data source and description of the SWAT modeling.
Datatype
Description
Source
DEM
U.S. Geological Survey’s (USGS) National Elevation Dataset (NED)
https://lta.cr.usgs.gov/NED
Land use
National Land Cover Database 2011 (NLCD 2011)
https://www.mrlc.gov/
Soil
Soil Survey Geographic Database
The ArcSWAT built-in soil database
Climate
Daily temperature and daily precipitation
https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/
Daily solar radiation, relative humidity, and wind speed
Weather generator data: WGEN_US_COOP_1960_2010The gauge-based climate dataset, including precipitation and maximum and minimum temperatures, was collected from National Oceanic and Atmospheric Administration (NOAA) National Climate Data Center (NCDC) Global Historical Climatology Network-Daily (GHCN-D) Version 3.24. The GHCN-D includes daily climate records from numerous sources that have been integrated and subjected to a common suite of quality assurance reviews (http://www.ncdc.noaa.gov) and thus selected for SWAT model calibration/validation as well as statistical downscaling in this study. The gauge-based climate measurements are unevenly distributed, which is often not ideal for large-scale hydrological simulation. Therefore, gauge-based climate data were gridded to a 0.385° × 0.385° spatial resolution and henceforth referred to as GHCN-0.385. The multiple nearest neighbors interpolating method developed by Liu et al. [34] was used in this study, averaging daily records of all weather stations within a certain distance or range (i.e., 20 km) from each node and assigning mean precipitation or temperature values to that node. This method is very similar to the pixel-based method wherein data are interpolated to a certain resolution grid (e.g., 0.385° by 0.385°) by averaging the available data over each grid area and assigning the mean value to the center of the grid [35]. Taking one grid of GHCN-0.385, for example, there are three weather stations within the distance of 20 km, viz., USC00117636, USC00112344, and US1ILJF0007 (Figure 1). The precipitation records of the three weather stations ranged from 1915 to 2018, from 1972 to 2008, and from 2012 to 2013, respectively. The interpolated records of the grid that ranged from 1971 to 2000 are composed of two parts: the first part that ranged from 1971 to 1972 is the same as the USC00117636 records and the second part that ranged from 1973 to 2000 is the average of USC00117636 and USC00112344.


#### 2.2.2. NCEP Reanalysis Data

Daily reanalysis data during 1971–2000 obtained from the NCEP/NCAR (National Centers for Environmental Prediction/National Center for Atmospheric Research) are usually regarded as a set of observed large-scale atmospheric variables and are used as the predictors in the downscaling model [36]. With a spatial resolution of 2.5° × 2.5°, the NCEP/NCAR includes many large-scale atmosphere variables such as mean sea level pressure, geopotential height, horizontal wind, specific humidity, relative humidity, and others. All of the variables were interpolated into the same spatial resolution as GHCN-0.385 (e.g., 0.385° by 0.385°) for downscaling model construction.


#### 2.2.3. GCM Variables for the Historical Period and Future Climate Scenarios

Based on the easy availability of daily data of the atmospheric variable from CMIP5, five GCMs were considered in this study (Table3): (i) GFDL-ESM2G developed by the NOAA Geophysical Fluid Dynamics Laboratory, USA; (ii) CanESM2 developed by the Canadian Centre for Climate Modeling, Canada; (iii) IPSL-CM5A-LR developed by the Institute Pierre-Simon Laplace, France; (iv) MIROC-ESM-CHEM developed by the National Institute for Environmental Studies, The University of Tokyo, Japan; and (v) NorESM1-M developed by the Norwegian Climate Centre, Norway. We chose three IPCC AR5 emission scenarios, RCP2.6, RCP4.5, and RCP8.5, to define greenhouse gas concentration trajectories in the year 2100 compared to preindustrial values of radiative forcing. RCP2.6 is the mitigation scenario with a forcing peaking in 2035 to approximately 3 W/m2 and decreasing to 2.6 W/m2 by 2100, which will lead to a very low forcing level. In contrast, RCP8.5 represents the highest emission scenario with an increasing radiative forcing to 8.5 W/m2 in 2100. RCP4.5 represents a medium-low emission scenario with a stabilization of radiative forcing at 4.5 W/m2 in 2100 [2]. Furthermore, two 30-year future time periods, viz., a midcentury period with a time frame of 2021–2050 and a late-century period with a time frame of 2070–2099, and the historical period with a time frame of 1971–2000 were considered.Table 3
Information of the five CMIP5 global climate models.
Institute
Nation
Modeling center (or group)
Model information
Model name
Atmosphere resolution
GFDL
USA
NOAA Geophysical Fluid Dynamics Laboratory
GFDL-ESM2G
C45 (∼2° × 2.5°) L24
CCCma
Canada
Canadian Centre for Climate Modeling
CanESM2
T63 (2.8125° × 2.8125°) L35
IPSL
France
Institute Pierre-Simon Laplace
IPSL-CM5A-LR
LMDZ (∼1.9° × 3.75°)
MIROC
Japan
National Institute for Environmental Studies, The University of Tokyo
MIROC-ESM-CHEM
T42 (∼2.8125° × 2.8125°) L80
NCC
Norway
Norwegian Climate Centre
NorESM1-M
F19 (∼1.875° × 2.5°) L26
Note. Expansions of all CMIP model acronyms are available online athttp://www.ametsoc.org/PubsAcronymList.


## 2.1. Study Basin

Originating in West Central Ohio, USA, the Wabash River flows southwest across Indiana to Southern Illinois, where it forms the Illinois-Indiana border before draining into the Ohio River (Figure1). It is the largest northern tributary of the Ohio River [32]. With a drainage area of 85,237 km2, the Wabash River Basin is a 4-digit HUC (hydrologic unit code) area formed by two 6-digit HUC subbasins, viz., the main Wabash River Basin and the Patoka-White River Basin. Developed by the U.S. Geological Survey (USGS) in cooperation with the U.S. Water Resource Council and the U.S. Department of Agriculture Natural Resources Conservation Service (NRCS), the HUC is a numbering system for watersheds to provide a common coding system for state and federal agencies. Every hydrologic unit is identified by a unique HUC consisting of 2 to 12 digits based on the six levels of classification in the hydrologic unit system [33]. Among them, HUC 4 represents the subregion level, delineating large river basins, and HUC 6 subdivides many of the subregions into basins known as accounting units. The Wabash River Basin lies in two physiographic provinces, the Central Lowland and the Interior Low Plateaus, and the altitude varies from 60 to 373 m a.s.l. (Figure 1). Predominant land use in the basin changes from agriculture to urban to forest moving from north to south. Approximately 65% of the total land is row-crop agriculture (predominantly corn and soybeans), while urban land coverage represents 2.3%. The types of soils in the Wabash River Basin include clay glacial till soils in the upper Wabash, loamy glacial tills to moderately thick loess shed, and shale in the lower watershed. The humid continental climate of the Wabash River Watershed is shaped by the presence of the Great Lakes and its location. Situated in the middle of the North American continent and far from the temperature-moderating effects of the oceans, the Wabash River Basin contributes to large seasonal variations in air temperature from hot, humid summers to cold winters. Precipitation in the Wabash River Basin is fairly well distributed throughout the year. The annual precipitation is approximately 1,141 mm in Indianapolis (IN). Over the entire basin, the values range from 940 mm in northern areas to 1,245 mm in southern areas. Long-term average temperatures range from 9.4°C in the north to 13.8°C in the south. The average monthly temperatures range from approximately 0°C in January to 25°C in July [32] (Table 1).Figure 1
Shaded-relief map of the Wabash River Watershed showing the stream channels delineated with ArcSWAT and the locations of the climate data nodes and selected U.S. Geological Survey (USGS) gauge stations. The numbering and exact locations of USGS gauge stations are shown in Table1.Table 1
Four USGS streamflow gauge stations located in the midstream and downstream of the main stream and major tributaries.
ID
Streamflow gauge stations
Catchment area
Latitude
Longitude
Websites
Station 1
Wabash River at Terre Haute
3.15 × 104 km2
39°27′57″N
87°25′11″W
https://waterdata.usgs.gov/in/nwis/uv?site_no=03341500
Station 2
White River at Petersburg
2.88 × 104 km2
38°30′39″N
87°17′22″W
https://waterdata.usgs.gov/nwis/inventory/?site_no=03374000
Station 3
Wabash River at Mount Carmel
7.39 × 104 km2
38°23′54″N
87°45′23″W
https://waterdata.usgs.gov/nwis/inventory/?site_no=03377500
Station 4
Little Wabash River at Carmel, IL
0.83 × 104 km2
38°30′39″N
87°17′22″W
https://waterdata.usgs.gov/nwis/inventory?agency_code=USGS&site_no=03381500


## 2.2. Data Used

### 2.2.1. Geospatial Data

Topography, land cover, and soil data are used to represent the catchment heterogeneity for the hydrologic analysis and were compiled from different sources (Table2). The USGS National Elevation Dataset (NED) digital elevation model with a 30 m spatial resolution furnished land surface topography for the simulation domain. Land use data were obtained from National Land Cover Database 2011 (NLCD 2011), while the Soil Survey Geographic Database (SSURGO) that incorporated in ArcSWAT and overlapping the study area was used.Table 2
Data source and description of the SWAT modeling.
Datatype
Description
Source
DEM
U.S. Geological Survey’s (USGS) National Elevation Dataset (NED)
https://lta.cr.usgs.gov/NED
Land use
National Land Cover Database 2011 (NLCD 2011)
https://www.mrlc.gov/
Soil
Soil Survey Geographic Database
The ArcSWAT built-in soil database
Climate
Daily temperature and daily precipitation
https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/
Daily solar radiation, relative humidity, and wind speed
Weather generator data: WGEN_US_COOP_1960_2010The gauge-based climate dataset, including precipitation and maximum and minimum temperatures, was collected from National Oceanic and Atmospheric Administration (NOAA) National Climate Data Center (NCDC) Global Historical Climatology Network-Daily (GHCN-D) Version 3.24. The GHCN-D includes daily climate records from numerous sources that have been integrated and subjected to a common suite of quality assurance reviews (http://www.ncdc.noaa.gov) and thus selected for SWAT model calibration/validation as well as statistical downscaling in this study. The gauge-based climate measurements are unevenly distributed, which is often not ideal for large-scale hydrological simulation. Therefore, gauge-based climate data were gridded to a 0.385° × 0.385° spatial resolution and henceforth referred to as GHCN-0.385. The multiple nearest neighbors interpolating method developed by Liu et al. [34] was used in this study, averaging daily records of all weather stations within a certain distance or range (i.e., 20 km) from each node and assigning mean precipitation or temperature values to that node. This method is very similar to the pixel-based method wherein data are interpolated to a certain resolution grid (e.g., 0.385° by 0.385°) by averaging the available data over each grid area and assigning the mean value to the center of the grid [35]. Taking one grid of GHCN-0.385, for example, there are three weather stations within the distance of 20 km, viz., USC00117636, USC00112344, and US1ILJF0007 (Figure 1). The precipitation records of the three weather stations ranged from 1915 to 2018, from 1972 to 2008, and from 2012 to 2013, respectively. The interpolated records of the grid that ranged from 1971 to 2000 are composed of two parts: the first part that ranged from 1971 to 1972 is the same as the USC00117636 records and the second part that ranged from 1973 to 2000 is the average of USC00117636 and USC00112344.


### 2.2.2. NCEP Reanalysis Data

Daily reanalysis data during 1971–2000 obtained from the NCEP/NCAR (National Centers for Environmental Prediction/National Center for Atmospheric Research) are usually regarded as a set of observed large-scale atmospheric variables and are used as the predictors in the downscaling model [36]. With a spatial resolution of 2.5° × 2.5°, the NCEP/NCAR includes many large-scale atmosphere variables such as mean sea level pressure, geopotential height, horizontal wind, specific humidity, relative humidity, and others. All of the variables were interpolated into the same spatial resolution as GHCN-0.385 (e.g., 0.385° by 0.385°) for downscaling model construction.


### 2.2.3. GCM Variables for the Historical Period and Future Climate Scenarios

Based on the easy availability of daily data of the atmospheric variable from CMIP5, five GCMs were considered in this study (Table3): (i) GFDL-ESM2G developed by the NOAA Geophysical Fluid Dynamics Laboratory, USA; (ii) CanESM2 developed by the Canadian Centre for Climate Modeling, Canada; (iii) IPSL-CM5A-LR developed by the Institute Pierre-Simon Laplace, France; (iv) MIROC-ESM-CHEM developed by the National Institute for Environmental Studies, The University of Tokyo, Japan; and (v) NorESM1-M developed by the Norwegian Climate Centre, Norway. We chose three IPCC AR5 emission scenarios, RCP2.6, RCP4.5, and RCP8.5, to define greenhouse gas concentration trajectories in the year 2100 compared to preindustrial values of radiative forcing. RCP2.6 is the mitigation scenario with a forcing peaking in 2035 to approximately 3 W/m2 and decreasing to 2.6 W/m2 by 2100, which will lead to a very low forcing level. In contrast, RCP8.5 represents the highest emission scenario with an increasing radiative forcing to 8.5 W/m2 in 2100. RCP4.5 represents a medium-low emission scenario with a stabilization of radiative forcing at 4.5 W/m2 in 2100 [2]. Furthermore, two 30-year future time periods, viz., a midcentury period with a time frame of 2021–2050 and a late-century period with a time frame of 2070–2099, and the historical period with a time frame of 1971–2000 were considered.Table 3
Information of the five CMIP5 global climate models.
Institute
Nation
Modeling center (or group)
Model information
Model name
Atmosphere resolution
GFDL
USA
NOAA Geophysical Fluid Dynamics Laboratory
GFDL-ESM2G
C45 (∼2° × 2.5°) L24
CCCma
Canada
Canadian Centre for Climate Modeling
CanESM2
T63 (2.8125° × 2.8125°) L35
IPSL
France
Institute Pierre-Simon Laplace
IPSL-CM5A-LR
LMDZ (∼1.9° × 3.75°)
MIROC
Japan
National Institute for Environmental Studies, The University of Tokyo
MIROC-ESM-CHEM
T42 (∼2.8125° × 2.8125°) L80
NCC
Norway
Norwegian Climate Centre
NorESM1-M
F19 (∼1.875° × 2.5°) L26
Note. Expansions of all CMIP model acronyms are available online athttp://www.ametsoc.org/PubsAcronymList.


## 2.2.1. Geospatial Data

Topography, land cover, and soil data are used to represent the catchment heterogeneity for the hydrologic analysis and were compiled from different sources (Table2). The USGS National Elevation Dataset (NED) digital elevation model with a 30 m spatial resolution furnished land surface topography for the simulation domain. Land use data were obtained from National Land Cover Database 2011 (NLCD 2011), while the Soil Survey Geographic Database (SSURGO) that incorporated in ArcSWAT and overlapping the study area was used.Table 2
Data source and description of the SWAT modeling.
Datatype
Description
Source
DEM
U.S. Geological Survey’s (USGS) National Elevation Dataset (NED)
https://lta.cr.usgs.gov/NED
Land use
National Land Cover Database 2011 (NLCD 2011)
https://www.mrlc.gov/
Soil
Soil Survey Geographic Database
The ArcSWAT built-in soil database
Climate
Daily temperature and daily precipitation
https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/
Daily solar radiation, relative humidity, and wind speed
Weather generator data: WGEN_US_COOP_1960_2010The gauge-based climate dataset, including precipitation and maximum and minimum temperatures, was collected from National Oceanic and Atmospheric Administration (NOAA) National Climate Data Center (NCDC) Global Historical Climatology Network-Daily (GHCN-D) Version 3.24. The GHCN-D includes daily climate records from numerous sources that have been integrated and subjected to a common suite of quality assurance reviews (http://www.ncdc.noaa.gov) and thus selected for SWAT model calibration/validation as well as statistical downscaling in this study. The gauge-based climate measurements are unevenly distributed, which is often not ideal for large-scale hydrological simulation. Therefore, gauge-based climate data were gridded to a 0.385° × 0.385° spatial resolution and henceforth referred to as GHCN-0.385. The multiple nearest neighbors interpolating method developed by Liu et al. [34] was used in this study, averaging daily records of all weather stations within a certain distance or range (i.e., 20 km) from each node and assigning mean precipitation or temperature values to that node. This method is very similar to the pixel-based method wherein data are interpolated to a certain resolution grid (e.g., 0.385° by 0.385°) by averaging the available data over each grid area and assigning the mean value to the center of the grid [35]. Taking one grid of GHCN-0.385, for example, there are three weather stations within the distance of 20 km, viz., USC00117636, USC00112344, and US1ILJF0007 (Figure 1). The precipitation records of the three weather stations ranged from 1915 to 2018, from 1972 to 2008, and from 2012 to 2013, respectively. The interpolated records of the grid that ranged from 1971 to 2000 are composed of two parts: the first part that ranged from 1971 to 1972 is the same as the USC00117636 records and the second part that ranged from 1973 to 2000 is the average of USC00117636 and USC00112344.


## 2.2.2. NCEP Reanalysis Data

Daily reanalysis data during 1971–2000 obtained from the NCEP/NCAR (National Centers for Environmental Prediction/National Center for Atmospheric Research) are usually regarded as a set of observed large-scale atmospheric variables and are used as the predictors in the downscaling model [36]. With a spatial resolution of 2.5° × 2.5°, the NCEP/NCAR includes many large-scale atmosphere variables such as mean sea level pressure, geopotential height, horizontal wind, specific humidity, relative humidity, and others. All of the variables were interpolated into the same spatial resolution as GHCN-0.385 (e.g., 0.385° by 0.385°) for downscaling model construction.


## 2.2.3. GCM Variables for the Historical Period and Future Climate Scenarios

Based on the easy availability of daily data of the atmospheric variable from CMIP5, five GCMs were considered in this study (Table3): (i) GFDL-ESM2G developed by the NOAA Geophysical Fluid Dynamics Laboratory, USA; (ii) CanESM2 developed by the Canadian Centre for Climate Modeling, Canada; (iii) IPSL-CM5A-LR developed by the Institute Pierre-Simon Laplace, France; (iv) MIROC-ESM-CHEM developed by the National Institute for Environmental Studies, The University of Tokyo, Japan; and (v) NorESM1-M developed by the Norwegian Climate Centre, Norway. We chose three IPCC AR5 emission scenarios, RCP2.6, RCP4.5, and RCP8.5, to define greenhouse gas concentration trajectories in the year 2100 compared to preindustrial values of radiative forcing. RCP2.6 is the mitigation scenario with a forcing peaking in 2035 to approximately 3 W/m2 and decreasing to 2.6 W/m2 by 2100, which will lead to a very low forcing level. In contrast, RCP8.5 represents the highest emission scenario with an increasing radiative forcing to 8.5 W/m2 in 2100. RCP4.5 represents a medium-low emission scenario with a stabilization of radiative forcing at 4.5 W/m2 in 2100 [2]. Furthermore, two 30-year future time periods, viz., a midcentury period with a time frame of 2021–2050 and a late-century period with a time frame of 2070–2099, and the historical period with a time frame of 1971–2000 were considered.Table 3
Information of the five CMIP5 global climate models.
Institute
Nation
Modeling center (or group)
Model information
Model name
Atmosphere resolution
GFDL
USA
NOAA Geophysical Fluid Dynamics Laboratory
GFDL-ESM2G
C45 (∼2° × 2.5°) L24
CCCma
Canada
Canadian Centre for Climate Modeling
CanESM2
T63 (2.8125° × 2.8125°) L35
IPSL
France
Institute Pierre-Simon Laplace
IPSL-CM5A-LR
LMDZ (∼1.9° × 3.75°)
MIROC
Japan
National Institute for Environmental Studies, The University of Tokyo
MIROC-ESM-CHEM
T42 (∼2.8125° × 2.8125°) L80
NCC
Norway
Norwegian Climate Centre
NorESM1-M
F19 (∼1.875° × 2.5°) L26
Note. Expansions of all CMIP model acronyms are available online athttp://www.ametsoc.org/PubsAcronymList.


## 3. Methodology and Model Description

A schematic diagram of statistical downscaling and climate change impact on SWAT streamflow simulation is shown in Figure2, and the process will be described in the following sections.Figure 2
Schematic illustration of the hydrological impact study framework (note: the methods are shown in italics).


### 3.1. ANN Downscaling Modeling

The screening of large-scale variables is the most important process in statistical downscaling. In this study, variables used for downscaling were screened out using the Statistical Downscaling Model (SDSM) [37] and are listed in Table 4. Then, the selected variables were used for rainfall and temperature (maximum temperature and minimum temperature) downscaling that was implemented using the ANN method.Table 4
NCEP and GCM daily predictor variables, codes, and descriptions.
Variable no.
Variable code
Description
1
hurs
Near-surface specific humidity
2
mslp
Mean sea-level pressure
3
ω5
Omega wind at 500 hPa
4
ω8
Omega wind at 850 hPa
5
p5_v
Meridional velocity component at 500 hPa
6
p5_u
Zonal velocity component at 500 hPa
7
p850
850 hPa geopotential height
8
r_700
Relative humidity at 700 hPa height
9
s_700
Specific humidity at 700 hPa height
10
s_850
Specific humidity at 850 hPa height
11
temp
Near-surface air temperatureAn artificial neuron network is a computational model based on the structure and functions of biological neural networks. The ANNs are analogous to multiple regressions but more efficient in finding and representing relationships between predictors and predictands with respect to the fact that they can actually learn from observing datasets. In other words, the ANN is used as a random function approximation tool that estimates the most cost-effective and ideal methods for arriving at solutions while defining computing functions or distributions. During this process, the ANN uses data samples rather than entire datasets to arrive at solutions, which saves time. The ANN has been widely used in fields of hydrology and climatology [38].In this study, the feedforward multilayer perceptron (MLP) network with a backpropagation (BP) error-correction learning algorithm was adopted. The ANN has three layers that are interconnected. The first layer (i.e., input layer) consists of input neurons. These neurons send data to the second layer (i.e., hidden layer), which subsequently sends the output neurons to the third layer (i.e., output layer). In practice, the number of hidden layers is usually determined using a trial-and-error procedure. The standard training algorithm used in the ANN is the backpropagation algorithm. The full mathematical derivation of this algorithm is found in the study of Hagan et al. [39]. The error between the predicted output and the observations is estimated using the root mean square error (RMSE). The calculated error is then backpropagated and will be minimized by iteratively changing the weights of inputs (i.e., predictors) until the error reaches a minimum value via the process known as training.A correlation coefficient is used for measuring the strength of the relationship between the relative movements of downscaled climate variables and their corresponding observations. Given the typical bad performance of downscaled precipitation, its frequency distribution was further evaluated. The frequency distribution of downscaled precipitation from the NCEP/NCAR and its observations were calculated based on the rainfall intensity classification of the World Meteorological Organization (WMO) standard [40]: (a) rain < 1 mm (no/minimal rain), (b) 1 mm ≤ rain < 2 mm (light rain), (c) 2 mm ≤ rain < 5 mm (low moderate rain), (d) 5 mm ≤ rain < 10 mm (high moderate rain), (e) 10 mm ≤ rain < 20 mm (low heavy rain), (f) 20 mm ≤ rain < 50 mm (high heavy rain), and (g) rain ≥ 50 mm (violent rain).Once the downscaling modeling is completed, the GCM atmospheric variables for the historical period and two future periods were first bias-corrected by adopting a quantile mapping method based on observations (e.g., NCEP/NCAR at the 2.5° by 2.5° resolution) [41]. The bias-corrected GCM outputs are then spatially disaggregated to the 0.385° × 0.385° grids (atmospheric horizontal spectral resolution of T382∼38 km) using inverse distance weighting interpolation and inputted into the calibrated and validated ANN model to generate daily precipitation, daily maximum temperature, and daily minimum temperature in the historical period and the two future periods.


### 3.2. The SWAT Hydrological Model

The physics-based SWAT model is a semispatially distributed and continuous hydrological model developed by the U.S. Department of Agriculture-Agriculture Research Service (USDA-ARS) in the early 1990s to assess the impact of management decisions on water, sediment, and agricultural chemical yields in large, complex watersheds with varying soil, land use, and management conditions over long periods [42]. The SWAT is able to operate on daily and subdaily time steps, and its spatial parameterization is performed by dividing the watershed into subbasins on the basis of topography, which are then further subdivided into a series of homogeneous hydrological response units (HRUs) based on soil, land use, and slope characteristics. In the SWAT, subbasins are connected by the stream network, while the HRUs represent certain percentages of the subbasins’ areas and are not spatially identified. Responses of each HRU, including dynamics of water, nutrients, sediments, and other chemicals, are simulated individually, aggregated at the subbasin level, and then routed to the associated reach and catchment outlets through the channel network.An ArcGIS graphical user interface and extension for the SWAT, i.e., ArcSWAT, was used to delineate the Wabash River Basin and write SWAT input files [43]. The whole basin was divided into 516 subbasins based on DEM data and into a total of 2,803 HRUs based on soil, land use, and slope in the subbasins. Climate data, including daily precipitation, daily maximum air temperature, and daily minimum air temperature, were inputted into the SWAT model. The curve number and Penman–Monteith equations were chosen for the simulation of infiltration and potential evapotranspiration, respectively. The baseline simulation was implemented at a daily time step spanning from 1971 to 2000, and the first five years (1971–1975) were used as a spin-up or warm-up period.


### 3.3. SUFI-2 Calibration and Uncertainty Analysis Procedure

In this study, a multisite calibration/validation scheme based on the single-objective optimization algorithm was used. Daily streamflow data from four USGS gauge stations located in the downstream portions of the Wabash River Basin and major tributaries were used for model calibration (January 1976–December 1991) and validation (January 1992–December 2000) (Table1; Figure 1). Similar to the previous study by Liu et al. [34], this study calibrated all key SWAT model parameters, and names and ranges are specified in Table 5.Table 5
Initial and final ranges, best estimations, andP values for SWAT calibration.
Parameters
Initial range
Final range
Best estimations
P values
r_CN2.mgt
−0.50/0.50
−0.24/0.11
0.07
0.00
v_ALPHA_BF.gw
0.00/1.00
0.47/0.82
0.75
0.57
v_GW_DELAY.gw
30.00/450.00
43.82/195.70
178.46
0.00
r_GWQMN.gw
−0.90/1.00
−0.20/0.60
0.21
0.03
v_ALPHA_BNK.rte
0.00/1.00
0.61/1.00
0.88
0.00
v_CH_N2.rte
0.00/0.30
0.05/0.16
0.09
0.77
v_CH_N1.sub
0.00/1.00
0.00/0.34
0.06
0.08
v_EPCO.hru
0.01/1.00
0.52/1.00
0.57
0.99
r_SLSUBBSN.hru
−0.50/0.65
−0.19/0.15
0.13
0.37
v_ESCO.hru
0.00/1.00
0.70/1.00
0.80
0.00
v_GW_REVAP.gw
0.02/0.20
0.05/0.15
0.07
0.00
v_RCHRG_DP.gw
0.00/1.00
0.10/0.53
0.11
0.00
r_OV_N.hru
−0.80/4.00
0.41/2.84
1.40
0.67
v_REVAPMN.gw
0.00/500.00
0.00/158.87
55.44
0.77
v_CH_K2.rte
0.00/500.00
0.00/195.14
156.21
0.02
r_SOL_K(1).sol
−0.90/9.00
−0.9/2.16
0.90
0.06
r_SOL_AWC(1).sol
−0.80/9.00
1.91/4.98
3.38
0.91
v_SFTMP.bsn
−5.00/5.00
−2.18/2.61
0.53
0.83
v_SMTMP.bsn
−5.00/5.00
1.90/5.00
2.92
0.76
r_SOL_BD(1).mgt
−0.50/0.60
0.10/0.57
0.54
0.01
CN2: SCS runoff curve number for moisture condition II; ALPHA_BF: baseflow alpha factor (days); GW_DELAY: groundwater delay time (days); GWQMN: threshold depth of water in the shallow aquifer required for the return flow to occur (mm); ALPHA_BNK: baseflow alpha factor for bank storage (days); CH_N2: Manning’s “n” value for the main channel; CH_N1: Manning’s “n” value for the tributary channel; EPCO: plant uptake compensation factor; SLSUBBSN: average slope length (m); ESCO: soil evaporation compensation factor; GW_REVAP: groundwater “revap” coefficient; RCHRG_DP: deep aquifer percolation fraction; OV_N: Manning’s “n” value for the overland flow; REVAPMN: threshold depth of water in the shallow aquifer for “revap” or percolation to the deep aquifer to occur (mm); CH_K2: effective hydraulic conductivity in the main channel alluvium (mm/hr); SOL_K: saturated hydraulic conductivity (mm/hr); SOL_AWC: available water capacity of the soil layer (mm); SFTMP: snowfall temperature (°C); SMTMP: snowmelt base temperature (°C); SOL_BD: moist bulk density (mg/m3). V_ means the existing parameter value is to be replaced by a given value; r_ means an existing parameter value is multiplied by (1 + a given value).A semiautomatic sequential uncertainty fitting algorithm (SUFI-2) procedure was used for parameter estimation and uncertainty propagation analysis [44]. Uncertainty associated with the driving parameters, conceptual models, and measured data is depicted on parameter ranges. These ranges are often called prior estimates because they are not conditioned on any input-output time-series data but the decisions from the modeler. Three performance metrics, including the Nash–Sutcliffe (NS) coefficient, R2, and PBIAS [44], were employed to evaluate the match or fitness of the simulation results. Based on NS coefficient values, for example, the model performance can be judged as follows: unsatisfactory (NS coefficient ≤ 0.50), satisfactory (0.50 < NS coefficient ≤ 0.65), good (0.65 < NS coefficient ≤ 0.75), or very good (0.75 < NS coefficient ≤ 1.00) [45]. The Latin hypercube (LH) sampling method was used to sample from the initial parameter distribution. Global sensitivity analysis was performed. A t-test was used to identify the relative significance of each parameter. The t-value provides a measure of sensitivity, where larger absolute values indicate higher sensitivity. The P value evaluates the significance of the sensitivity, where a value close to zero indicates greater significance.After each iteration, new and narrow parameter ranges were identified, leading to uncertainty reduction in the model output. Two factors, theP-factor and the R-factor, were used to quantify the performance of uncertainty analysis. The P-factor represents the percentage of observed data bracketed by 95% prediction uncertainty (95PPU). The R-factor calculates the average width of the 95% uncertainty band divided by the standard deviation of the corresponding measured variable. A P-factor (ranging from zero to 100%) of 100% and an R-factor (ranging from zero to infinity) of zero indicate an exact match between the simulations and observations, which would never occur in practice given the inherent uncertainty of the model. A larger P-factor can be achieved at the expense of a larger R-factor, and a compromise between the two factors must be reached. A more detailed description of the calibration protocol can be found in the study of Abbaspour [44]. For streamflows, P-factor >0.7 and R-factor <1.5 are considered acceptable in terms of prediction uncertainty.


### 3.4. Multimodel Ensemble Mean Projection for the Future Scenario

The ensemble mean of simulated streamflows from the hydrologic model obtained by using precipitation and maximum/minimum temperature data from the five GCMs in the historical period was first evaluated by comparing it to that of the observed streamflow. Station 3 covering the majority of the drainage area of the Wabash River Basin was selected for this evaluation. The Q-Q plots that illustrate the distribution of simulated annual and daily streamflows during 1976–2000 with downscaled precipitation and maximum/minimum temperature from each GCM as the driving force to the observed values were adopted for the comparison.Typically, the downscaled climate variables, especially precipitation, are seldom consistent with the observations. Moreover, the hydrological model is highly sensitive to precipitation. Therefore, the ensemble mean of simulated streamflows from the hydrologic model obtained by using precipitation and maximum/minimum temperature data from the five GCMs is used to evaluate the hydrologic impact of climate change in the basin by comparing hydrologic responses in the future time periods with those obtained under the historical scenario. For the future simulation, the first five years (i.e., 2021–2025 in the midcentury period and 2070–2074 in the late-century period) were used as spin-up periods.


## 3.1. ANN Downscaling Modeling

The screening of large-scale variables is the most important process in statistical downscaling. In this study, variables used for downscaling were screened out using the Statistical Downscaling Model (SDSM) [37] and are listed in Table 4. Then, the selected variables were used for rainfall and temperature (maximum temperature and minimum temperature) downscaling that was implemented using the ANN method.Table 4
NCEP and GCM daily predictor variables, codes, and descriptions.
Variable no.
Variable code
Description
1
hurs
Near-surface specific humidity
2
mslp
Mean sea-level pressure
3
ω5
Omega wind at 500 hPa
4
ω8
Omega wind at 850 hPa
5
p5_v
Meridional velocity component at 500 hPa
6
p5_u
Zonal velocity component at 500 hPa
7
p850
850 hPa geopotential height
8
r_700
Relative humidity at 700 hPa height
9
s_700
Specific humidity at 700 hPa height
10
s_850
Specific humidity at 850 hPa height
11
temp
Near-surface air temperatureAn artificial neuron network is a computational model based on the structure and functions of biological neural networks. The ANNs are analogous to multiple regressions but more efficient in finding and representing relationships between predictors and predictands with respect to the fact that they can actually learn from observing datasets. In other words, the ANN is used as a random function approximation tool that estimates the most cost-effective and ideal methods for arriving at solutions while defining computing functions or distributions. During this process, the ANN uses data samples rather than entire datasets to arrive at solutions, which saves time. The ANN has been widely used in fields of hydrology and climatology [38].In this study, the feedforward multilayer perceptron (MLP) network with a backpropagation (BP) error-correction learning algorithm was adopted. The ANN has three layers that are interconnected. The first layer (i.e., input layer) consists of input neurons. These neurons send data to the second layer (i.e., hidden layer), which subsequently sends the output neurons to the third layer (i.e., output layer). In practice, the number of hidden layers is usually determined using a trial-and-error procedure. The standard training algorithm used in the ANN is the backpropagation algorithm. The full mathematical derivation of this algorithm is found in the study of Hagan et al. [39]. The error between the predicted output and the observations is estimated using the root mean square error (RMSE). The calculated error is then backpropagated and will be minimized by iteratively changing the weights of inputs (i.e., predictors) until the error reaches a minimum value via the process known as training.A correlation coefficient is used for measuring the strength of the relationship between the relative movements of downscaled climate variables and their corresponding observations. Given the typical bad performance of downscaled precipitation, its frequency distribution was further evaluated. The frequency distribution of downscaled precipitation from the NCEP/NCAR and its observations were calculated based on the rainfall intensity classification of the World Meteorological Organization (WMO) standard [40]: (a) rain < 1 mm (no/minimal rain), (b) 1 mm ≤ rain < 2 mm (light rain), (c) 2 mm ≤ rain < 5 mm (low moderate rain), (d) 5 mm ≤ rain < 10 mm (high moderate rain), (e) 10 mm ≤ rain < 20 mm (low heavy rain), (f) 20 mm ≤ rain < 50 mm (high heavy rain), and (g) rain ≥ 50 mm (violent rain).Once the downscaling modeling is completed, the GCM atmospheric variables for the historical period and two future periods were first bias-corrected by adopting a quantile mapping method based on observations (e.g., NCEP/NCAR at the 2.5° by 2.5° resolution) [41]. The bias-corrected GCM outputs are then spatially disaggregated to the 0.385° × 0.385° grids (atmospheric horizontal spectral resolution of T382∼38 km) using inverse distance weighting interpolation and inputted into the calibrated and validated ANN model to generate daily precipitation, daily maximum temperature, and daily minimum temperature in the historical period and the two future periods.


## 3.2. The SWAT Hydrological Model

The physics-based SWAT model is a semispatially distributed and continuous hydrological model developed by the U.S. Department of Agriculture-Agriculture Research Service (USDA-ARS) in the early 1990s to assess the impact of management decisions on water, sediment, and agricultural chemical yields in large, complex watersheds with varying soil, land use, and management conditions over long periods [42]. The SWAT is able to operate on daily and subdaily time steps, and its spatial parameterization is performed by dividing the watershed into subbasins on the basis of topography, which are then further subdivided into a series of homogeneous hydrological response units (HRUs) based on soil, land use, and slope characteristics. In the SWAT, subbasins are connected by the stream network, while the HRUs represent certain percentages of the subbasins’ areas and are not spatially identified. Responses of each HRU, including dynamics of water, nutrients, sediments, and other chemicals, are simulated individually, aggregated at the subbasin level, and then routed to the associated reach and catchment outlets through the channel network.An ArcGIS graphical user interface and extension for the SWAT, i.e., ArcSWAT, was used to delineate the Wabash River Basin and write SWAT input files [43]. The whole basin was divided into 516 subbasins based on DEM data and into a total of 2,803 HRUs based on soil, land use, and slope in the subbasins. Climate data, including daily precipitation, daily maximum air temperature, and daily minimum air temperature, were inputted into the SWAT model. The curve number and Penman–Monteith equations were chosen for the simulation of infiltration and potential evapotranspiration, respectively. The baseline simulation was implemented at a daily time step spanning from 1971 to 2000, and the first five years (1971–1975) were used as a spin-up or warm-up period.


## 3.3. SUFI-2 Calibration and Uncertainty Analysis Procedure

In this study, a multisite calibration/validation scheme based on the single-objective optimization algorithm was used. Daily streamflow data from four USGS gauge stations located in the downstream portions of the Wabash River Basin and major tributaries were used for model calibration (January 1976–December 1991) and validation (January 1992–December 2000) (Table1; Figure 1). Similar to the previous study by Liu et al. [34], this study calibrated all key SWAT model parameters, and names and ranges are specified in Table 5.Table 5
Initial and final ranges, best estimations, andP values for SWAT calibration.
Parameters
Initial range
Final range
Best estimations
P values
r_CN2.mgt
−0.50/0.50
−0.24/0.11
0.07
0.00
v_ALPHA_BF.gw
0.00/1.00
0.47/0.82
0.75
0.57
v_GW_DELAY.gw
30.00/450.00
43.82/195.70
178.46
0.00
r_GWQMN.gw
−0.90/1.00
−0.20/0.60
0.21
0.03
v_ALPHA_BNK.rte
0.00/1.00
0.61/1.00
0.88
0.00
v_CH_N2.rte
0.00/0.30
0.05/0.16
0.09
0.77
v_CH_N1.sub
0.00/1.00
0.00/0.34
0.06
0.08
v_EPCO.hru
0.01/1.00
0.52/1.00
0.57
0.99
r_SLSUBBSN.hru
−0.50/0.65
−0.19/0.15
0.13
0.37
v_ESCO.hru
0.00/1.00
0.70/1.00
0.80
0.00
v_GW_REVAP.gw
0.02/0.20
0.05/0.15
0.07
0.00
v_RCHRG_DP.gw
0.00/1.00
0.10/0.53
0.11
0.00
r_OV_N.hru
−0.80/4.00
0.41/2.84
1.40
0.67
v_REVAPMN.gw
0.00/500.00
0.00/158.87
55.44
0.77
v_CH_K2.rte
0.00/500.00
0.00/195.14
156.21
0.02
r_SOL_K(1).sol
−0.90/9.00
−0.9/2.16
0.90
0.06
r_SOL_AWC(1).sol
−0.80/9.00
1.91/4.98
3.38
0.91
v_SFTMP.bsn
−5.00/5.00
−2.18/2.61
0.53
0.83
v_SMTMP.bsn
−5.00/5.00
1.90/5.00
2.92
0.76
r_SOL_BD(1).mgt
−0.50/0.60
0.10/0.57
0.54
0.01
CN2: SCS runoff curve number for moisture condition II; ALPHA_BF: baseflow alpha factor (days); GW_DELAY: groundwater delay time (days); GWQMN: threshold depth of water in the shallow aquifer required for the return flow to occur (mm); ALPHA_BNK: baseflow alpha factor for bank storage (days); CH_N2: Manning’s “n” value for the main channel; CH_N1: Manning’s “n” value for the tributary channel; EPCO: plant uptake compensation factor; SLSUBBSN: average slope length (m); ESCO: soil evaporation compensation factor; GW_REVAP: groundwater “revap” coefficient; RCHRG_DP: deep aquifer percolation fraction; OV_N: Manning’s “n” value for the overland flow; REVAPMN: threshold depth of water in the shallow aquifer for “revap” or percolation to the deep aquifer to occur (mm); CH_K2: effective hydraulic conductivity in the main channel alluvium (mm/hr); SOL_K: saturated hydraulic conductivity (mm/hr); SOL_AWC: available water capacity of the soil layer (mm); SFTMP: snowfall temperature (°C); SMTMP: snowmelt base temperature (°C); SOL_BD: moist bulk density (mg/m3). V_ means the existing parameter value is to be replaced by a given value; r_ means an existing parameter value is multiplied by (1 + a given value).A semiautomatic sequential uncertainty fitting algorithm (SUFI-2) procedure was used for parameter estimation and uncertainty propagation analysis [44]. Uncertainty associated with the driving parameters, conceptual models, and measured data is depicted on parameter ranges. These ranges are often called prior estimates because they are not conditioned on any input-output time-series data but the decisions from the modeler. Three performance metrics, including the Nash–Sutcliffe (NS) coefficient, R2, and PBIAS [44], were employed to evaluate the match or fitness of the simulation results. Based on NS coefficient values, for example, the model performance can be judged as follows: unsatisfactory (NS coefficient ≤ 0.50), satisfactory (0.50 < NS coefficient ≤ 0.65), good (0.65 < NS coefficient ≤ 0.75), or very good (0.75 < NS coefficient ≤ 1.00) [45]. The Latin hypercube (LH) sampling method was used to sample from the initial parameter distribution. Global sensitivity analysis was performed. A t-test was used to identify the relative significance of each parameter. The t-value provides a measure of sensitivity, where larger absolute values indicate higher sensitivity. The P value evaluates the significance of the sensitivity, where a value close to zero indicates greater significance.After each iteration, new and narrow parameter ranges were identified, leading to uncertainty reduction in the model output. Two factors, theP-factor and the R-factor, were used to quantify the performance of uncertainty analysis. The P-factor represents the percentage of observed data bracketed by 95% prediction uncertainty (95PPU). The R-factor calculates the average width of the 95% uncertainty band divided by the standard deviation of the corresponding measured variable. A P-factor (ranging from zero to 100%) of 100% and an R-factor (ranging from zero to infinity) of zero indicate an exact match between the simulations and observations, which would never occur in practice given the inherent uncertainty of the model. A larger P-factor can be achieved at the expense of a larger R-factor, and a compromise between the two factors must be reached. A more detailed description of the calibration protocol can be found in the study of Abbaspour [44]. For streamflows, P-factor >0.7 and R-factor <1.5 are considered acceptable in terms of prediction uncertainty.


## 3.4. Multimodel Ensemble Mean Projection for the Future Scenario

The ensemble mean of simulated streamflows from the hydrologic model obtained by using precipitation and maximum/minimum temperature data from the five GCMs in the historical period was first evaluated by comparing it to that of the observed streamflow. Station 3 covering the majority of the drainage area of the Wabash River Basin was selected for this evaluation. The Q-Q plots that illustrate the distribution of simulated annual and daily streamflows during 1976–2000 with downscaled precipitation and maximum/minimum temperature from each GCM as the driving force to the observed values were adopted for the comparison.Typically, the downscaled climate variables, especially precipitation, are seldom consistent with the observations. Moreover, the hydrological model is highly sensitive to precipitation. Therefore, the ensemble mean of simulated streamflows from the hydrologic model obtained by using precipitation and maximum/minimum temperature data from the five GCMs is used to evaluate the hydrologic impact of climate change in the basin by comparing hydrologic responses in the future time periods with those obtained under the historical scenario. For the future simulation, the first five years (i.e., 2021–2025 in the midcentury period and 2070–2074 in the late-century period) were used as spin-up periods.


## 4. Results and Discussion

### 4.1. Calibration and Validation of Hydrological Model Using Observed Climate Datasets

The parameter uncertainty definitely leads to uncertainties in model simulation. The best-fit parameter values and the sensitivities of model parameters are listed in Table3. The main sources of streamflow uncertainty were the SCS runoff curve number for moisture condition II (CN2), soil evaporation compensation factor (ESCO), groundwater “revap” coefficient (GW_REVAP), deep aquifer percolation fraction (RCHRG_DP), baseflow alpha factor for bank storage (ALPHA_BNK), groundwater delay time (GW_DELAY), moist bulk density (SOL_BD), effective hydraulic conductivity in the main channel alluvium (CH_K2), and threshold depth of water in the shallow aquifer required for the return flow to occur (GWQMN). The P values of these factors are less than 0.05. These results were basically consistent with the study [15].As shown in Figure3, the model reproduced observed river streamflow records fairly well albeit with slight underestimation of the peak flow values. The evaluation metrics of the simulated daily streamflow are listed in Table 6. The NS coefficient values ranged from 0.74 to 0.79 and 0.74 to 0.80 in the calibration period and validation period, respectively, indicating that good performance or very good performance was achieved. With respect to PBIAS values, the streamflows were slightly underestimated at Stations 1, 2, and 3, with PBIAS values ranging from 6.7 to 8.8 in calibration and 0.94 to 10.24 in validation. In contrast, system errors at Station 4 suggested that streamflows were highly overestimated in both calibration and validation periods. Similar results were found in P-factor values. Only Station 4 was below the acceptable level of 0.70, which indicated a relatively large uncertainty in this station, blaming the parameter-adjusting scheme provided by SWAT-CUP. The scheme calibrates the parameters using a single global modification term that scales the initial estimates by values, multiplication or additive terms and thus is not sufficient to cope with the heterogeneity of the precipitation input.Figure 3
Comparisons of daily streamflow simulations (black lines) with observations (gray areas) during calibration (1976–1991) and validation (1992–2000) periods.Table 6
Calibration (1976–1991) and validation (1992–2000) performance metrics for daily runoff at the four monitoring points.
ID
Calibration (1976–1991)
Validation (1992–2000)
R
2
NS coefficient
PBIAS
P-factor
R-factor
R
2
NS coefficient
PBIAS
Station 1
0.75
0.74
8.8
0.89
0.92
0.77
0.74
10.24
Station 2
0.79
0.77
7.6
0.87
0.82
0.82
0.80
0.94
Station 3
0.81
0.79
6.7
0.88
0.94
0.82
0.80
3.77
Station 4
0.80
0.78
−18.3
0.51
0.88
0.84
0.83
−12.50


### 4.2. Calibration and Validation of ANN Downscaling Model

The Pearson correlation coefficient between the observed climate and downscaled variables, including daily maximum temperature, daily minimum temperature, and daily precipitation, from the NCEP/NCAR for the calibration period and validation period was used to evaluate the performance of downscaled climate variables. For daily maximum temperature and daily minimum temperature, the correlation coefficients were all greater than 0.96 in the calibration period and greater than 0.94 in the validation period. The results clearly showed that the magnitudes of the downscaled 30-year daily maximum/minimum temperature values are statistically similar to those of the observations.Annual cycles of mean monthly maximum/minimum temperatures averaged for observations and downscaled from the NCEP/NCAR dataset over the basin during the period of 1971–2000 are shown in Figure4(a). Downscaled maximum temperatures from the NCEP/NCAR basically reproduced the observations with an absolute error less than 0.3°C, except for January, April, and August. In January and August, downscaled NCEP/NCAR mean monthly maximum temperatures were 2.53 and 28.88°C, respectively, which slightly exceeded the observations of 1.95 and 28.54, respectively. In April, the NCEP/NCAR underestimated the mean monthly maximum temperature with its value of 17.20, while the observed value is 17.70°C. Regarding the mean monthly minimum temperature, the downscaled NCEP/NCAR overestimated it in January, February, March, April, August, and September but underestimated it in the other months. The maximum over- and underestimation occurred in February and December with the absolute error values of 0.86 and 0.59, respectively.Figure 4
Annual cycle of mean monthly (a) maximum (solid lines) and minimum (dashed lines) temperatures and (b) precipitation averaged over the Wabash River Watershed for observation and downscaled from the NCEP/NCAR and from five GCMs during the period 1971–2000.
(a)
(b)Correlation coefficients (CCs) between the observed rainfall records and downscaled daily precipitation from the NCEP/NCAR for the calibration period and validation period are shown in Figure5. The values of the CC ranged from 0.63 to 0.79 (with an average of 0.71) in the calibration period except for one node in the western portion with a value of 0.58. The performance of downscaled precipitation in the validation period was slightly worse than that in the calibration period based on the CC values. The CC values ranged from 0.6 to 0.75 (averaged 0.68) during this period with the exception of two nodes located in the eastern and northern portions of the basin with values of 0.47 and 0.50, respectively. Overall, the results showed that the downscaled daily precipitation from the NCEP/NCAR was highly correlated and moderately correlated with observations, respectively.Figure 5
Pearson correlation coefficient between the observations and downscaled maximum temperature (a, b), minimum temperature (c, d), and daily precipitation (e, f) for the calibration period and validation period.
(a)
(b)
(c)
(d)
(e)
(f)The monthly mean precipitation cycles of downscaled precipitation from the NCEP/NCAR and observations over the Wabash River Watershed during the period of 1971–2000 are compared in Figure4(b). Downscaled precipitation captured the bimodal distribution of the annual cycle fairly well, though with different levels of bias and errors in each month. The mean monthly precipitation was overestimated in most of the calendar months, with the exception of October, during which the mean monthly precipitation was slightly underestimated. The discrepancies of monthly mean precipitation between downscaled precipitation and observations were significant in winter months (November, January, and February) and summer months (June, July, August, and September). As a result, the average annual precipitation was 1133.87 mm of NCEP/NCAR downscaled rainfall, which exceeded the observation (1057.31 mm).Five nodes in the precipitation grids of the observation that were evenly distributed over the study basin were selected for frequency distribution analysis and comparison (Figure5). In terms of frequency, the downscaled precipitation from the NCEP/NCAR to a great extent underestimated no/minimal rain events, ranging from 56.82% (at Node 4) to 59.07% (Node 5) compared with the observations that ranged from 70.33% (at Node 1) to 74.72% (Node 5). In addition, the events of high heavy rain and violent rain were also underestimated in this region, and the discrepancies were considerably discernible for violent rain events. The frequency of violent events ranged from 0.05% (at Nodes 3 and 5) to 0.24% (Node 1) for precipitation downscaled from the NCEP/NCAR but ranged from 0.26% (Node 5) to 0.47% (Node 1) for observations. It is notable that the discrepancy of violent events showed an increasing tendency along a south-north gradient. The violent rain events were 48.9% underestimated at Node 1 but 82.8% underestimated at Node 3. On the contrary, the occurrence days of light rain, low moderate rain, high moderate rain, and low heavy rain were significantly overestimated (Figure 6). The results were consistent with the study [46]. Moreover, the frequency distribution of downscaled precipitation records was very similar to that of the NCEP/CFSR over the Wabash River Basin during 1979 to 2000 (Figure 6), and numerous studies have documented the poor performance of the NCEP/CFSR as the driving force of the hydrologic model [47]. Given the high nonlinearity of watershed systems, even small changes in precipitation can cause tremendous variations in watershed response. Good performance of precipitation records in reproducing the observed annual average rainfall does not guarantee good performance in reproducing the streamflow. Thus, we anticipate here that significant differences in streamflow simulations will be observed when the calibrated hydrologic model uses downscaled GCMs’ precipitation as inputs.Figure 6
Frequency distribution of observed daily precipitation values and downscaled precipitation from the NCEP/NCAR and NCEP/CFSR at selected nodes (Figure5(e)).


### 4.3. Hydrological Evaluation of Downscaled Precipitation and Temperature from GCMs in the Historical Period

The calibrated and validated ANN model was then subjected to the GCMs in the historical period to evaluate the performance of each GCM in representing the observations and driving the SWAT hydrologic model. The annual cycle of mean monthly downscaled maximum/minimum temperatures from the five GCMs in the historical period is plotted in Figure4(a). It is obvious that the mean monthly maximum temperature reproduced by the GCMs is consistent with observations except IPSL-CM5A-LR, which to a large extent overestimated the maximum temperature in January, February, May, November, and December with the bias ranging from 1.68 to 3.33. With respect to the downscaled maximum temperature from the other four GCMs, all the absolute errors of mean monthly maximum temperature lied in the range of 0.01 to 1.18°C. Good matching was also witnessed for minimum temperature, with the absolute errors ranging from 0.001 to 1.74°C (Figure 4).Similar to the downscaled precipitation from the NCEP/NCAR, precipitation from the five GCMs basically captured the bimodal distribution of the annual cycle fairly well, albeit with different amplitudes (Figure4(b)). CanESM2 had a tendency to greatly underestimate the mean monthly precipitation in summer season months (May, June, July, and August), leading to an overall underestimation of the annual average rainfall (977.29 mm), which was contrary to NorESM1-M that to some extent overestimated the precipitation in these months and consequently overestimated the annual rainfall with the value of 1091.59 mm. GFDL-ESM2G exhibited improved performance in the second half of the year than that in the first half, except for September. In the first half year, GFDL-ESM2G overestimated the precipitation in January and June but underestimated that in the other months. The total error that resulted from both the positive bias and negative bias nearly cancelled each other, which led to a slight underestimation of the annual average rainfall (1026.56 mm). MIROC-ESM-CHEM produced an annual average rainfall of 1087.85 mm, which slightly exceeded the observations. The maximum error occurred in January, with a positive bias of 33.78 mm. IPSL-CM5A-LR performed best in reproducing the annual cycle of mean monthly precipitation as well as the annual average rainfall. All the bias was less than 10 mm for mean monthly precipitation except January, September, and December with values of 11.42, −11.96, and 19.05, respectively. The annual average rainfall was predicted to be 1051.42 mm, which is comparable to the observation of 1057.31 mm.The downscaled precipitation and maximum/minimum temperature from the five GCMs were then used to drive the calibrated SWAT model, respectively. The SWAT simulated annual streamflow produced by the five GCM projections indicated a wide range of annual runoff and flowrates (Figure7). As expected, all the GCMs except CanESM2 in the 8th percentile to a large extent overestimated the annual streamflow (Figure 7(a)), which was mainly blamed on the overestimation of the daily streamflow less than 2000 m3/s (Figure 7(b)). As shown in Table 7, the annual streamflow under the 90th percentile was substantially overestimated, while the ensemble mean of simulations with climate variables downscaled from five GCMs as inputs fit well for streamflows observed from 2000 to 4000 m3/s (approximately 90th percentile to 99th percentile). Beyond that, the streamflow rate was slightly underestimated. Overall, the downscaled precipitation and maximum/minimum temperature from the five GCMs as driving forces of the established hydrological model resulted in relatively high overestimation of the observed streamflow throughout the simulation period. The large temporal variation of precipitation (Figure 6) is present and thus leads to streamflow overestimation by impacting the process of runoff generation when using the hydrological model.Figure 7
Q-Q plots of (a) annual runoff and (b) daily streamflow for Station 3.
(a)
(b)Table 7
Annual streamflow statistics simulated using the five GCMs for the historical period and observation.
Annual streamflow (107 m3)
1
2
3
4
Observation
2731.61
1812.68
644.33
274.60
CanESM2
3416.56
2482.26
675.49
258.82
GFDL-ESM2G
3755.29
2855.95
632.29
267.05
IPSL-CM5A-LR
3766.90
2830.17
682.66
254.07
MIROC-ESM-CHEM
4040.99
3108.98
681.87
250.14
NorESM1-M
4262.30
3265.12
709.80
287.38
Ensemble
3848.41
2908.50
676.42
263.49
1: average annual runoff; 2: annual streamflow under the 90th percentile; 3: 90th–99th percentile; 4: 99th exceedance percentile.


### 4.4. Impact of Climate Change on Streamflow

The downscaled future precipitation and temperature for the five GCMs under the three RCP scenarios were inputted to the calibrated SWAT, separately, to investigate corresponding streamflow responses. Figure8 presents the ensemble mean of average annual runoff of Station 3 for the historical, midcentury, and late-century timeframes under three RCP scenarios, from which some interesting implications can be drawn. All projected future changes in annual runoff under RCP scenarios showed a decreasing trend to different degrees. In the midcentury period, the multimodel ensemble means were projected to decrease by 13.00%, 17.59%, and 6.91% under RCP2.6, RCP4.5, and RCP8.5 scenarios, respectively, compared with the historical period. Considerably, more significant decreases were witnessed in the late-century period, with the decrease rates of 25.29%, 27.61%, and 15.04% for the three scenarios (Table 8).Figure 8
Annual cycle of the monthly mean streamflow for the ensemble mean of SWAT simulations in Station 3 with downscaled climate variables from five GCMs as inputs over the period of 1976–2000.Table 8
Statistics of annual runoff under three scenarios for the two future periods and their changes compared to that under the historical scenario.
Scenarios
Annual streamflow (107 m3) and changes
RCP2.6
RCP4.5
RCP8.5
Midcentury
3354.64 (−13.00%)
3177.81 (−17.59%)
3589.39 (−6.91%)
Late century
2880.90(−25.29%)
2791.31 (−27.61%)
3276.05 (−15.04%)The annual cycle of the simulated multimodel ensemble mean of monthly averaged streamflows is also shown in Figure8. Under RCP2.6 and RCP4.5 scenarios, mean monthly streamflows in all calendar months decreased to varying degrees. The decline rates ranged from 0.40% to 26.03% for RCP2.6 scenarios in the midcentury and from 5.04% to 39.00% in the late century. For RCP4.5 scenarios, however, the ranges were 6.68%–34.12% and 7.51%–42.74%, respectively. The biggest decreases all occurred in November to March. Generally, the RCP2.6 and RCP4.5 scenarios projected very similar annual runoff and monthly runoff change, which were very different from those of RCP8.5. In the midcentury, the mean monthly streamflows under the RCP8.5 scenario were slightly decreased in most months, except for April and June, when the values were increased to be comparable with that of the historical scenario. In the late century, the mean monthly streamflows also showed an increasing trend in April and May with a decreasing trend in other months. Similar to RCP2.6 and RCP4.5, the biggest decline under the RCP8.5 scenario appeared in November to March.The projected annual runoff change was neither positively nor negatively related to the CO2 emission scenarios (i.e., RCP scenarios), which may be attributed to the changes in both of the rainfall pattern projected towards the midcentury and late century and an increase in the evaporation demand in terms of the climate warming. To distinguish the impacts of the projected change in precipitation and temperature on the simulated streamflow, the annual cycle of multimodel ensemble means of precipitation and maximum/minimum temperature under three scenarios for the two future periods is shown in Figure 9. The multimodel ensemble means of annual precipitation under RCP2.6 and RCP4.5 scenarios for the two future periods were 1023.37 (RCP26-midcentury), 1079.91 (RCP26-late century), 999.27 (RCP45-midcentury), and 1077.08 (RCP45-late century), respectively, which were very close to those of the historical period (1046.94 mm) and the observation (1057.31 mm). No significant change in precipitation was projected under RCP8.5 for the midcentury period with the annual average rainfall of 1076.87 mm; however, the precipitation under RCP8.5 for the late-century period was projected to have a 14.23% increase compared with the rainfall in the historical period. With respect to the maximum temperature, RCP2.6 and RCP4.5 projected an increase of 1.91°C (RCP26-midcentury), 1.99°C (RCP26-late century), 2.12°C (RCP45-midcentury), and 3.15°C (RCP45-late century), while RCP8.5 projected an increase of 2.32°C for the midcentury period and a significant increase of 5.31°C for the late-century period. Similar results were found for minimum temperature with increased values of 1.64°C (RCP26-midcentury), 1.80°C (RCP26-late century), 1.86°C (RCP45-midcentury), 2.88°C (RCP45-late century), 2.10°C (RCP85-midcentury), and 4.92°C (RCP85-late century). All of the above projection results are consistent with those of the IPCC report over this basin [2]. Thus, we can assert that climate warming plays an important role for RCP2.6 and RCP4.5, especially for the late-century period. For the RCP8.5 cases, during the midcentury period, the projected runoff slightly decreased though the changes in both the rainfall and temperature, very similar to the RCP26-midcentury case. We predicted that the rainfall pattern would be very different from the present pattern. During the late-century period, RCP8.5 projected a decrease of monthly streamflow with the exception of April, May, and August. Both precipitation and temperature are primary factors affecting the decrease in the streamflow.Figure 9
Annual cycle of mean monthly (a) precipitation and (b) maximum (solid lines) and minimum (dashed lines) temperatures averaged over the Wabash River Watershed for observation and multimodel ensembles for the historical period, midcentury, and late century under RCP scenarios.
(a)
(b)Regarding the mean monthly precipitation, an interesting finding was that all the cases, except for RCP4.5 and RCP8.5 during the late-century period, witnessed a decrease of rainfall in the summer season (June, July, and August), during which crops grow. Moreover, all the cases with the exception of RCP2.6 and RCP4.5 for the midcentury period indicated that the rainfall tends to increase in the winter and spring seasons.Overall, all scenarios projected a decrease of the streamflow with different amplitudes. These results were basically consistent with the study by Chien et al. [15]. Most of RCP2.6 and RCP4.5 scenarios projected an increase in precipitation in November to February, while the simulated streamflows will decrease in these months. These results suggest that temperature is a very important factor in the hydrological regime in this area. In summer months, the decrease in the streamflow was not serious with respect to the decrease of precipitation and increase of temperature. Regarding the RCP8.5 scenario, especially for the late-century period, large uncertainties derived from GCMs and the ANN downscaling method were witnessed. Therefore, further quantifying investigation of the uncertainties is necessary for reliable assessment of climate change impacts.


### 4.5. Uncertainty Analysis

Studies on climate change impact are always accompanied with uncertainties derived from GCMs, downscaling methods, and hydrological modeling. In this study, five GCMs were selected and used to drive the hydrological model. In view of the wide range of projected climate variables, especially precipitation, the bias was still large although bias correction was performed. With respect to the annual average rainfall, the absolute error was estimated ranging from 0.01 to 3.24 mm. The downscaling process is another important source of uncertainties. For this study, the frequency distribution of downscaled rainfall from GCMs was very similar to that of the NCEP/CFSR and deviated significantly from the observation. The mismatch of frequency distribution of precipitation led to a dramatic overestimation of the projected streamflow. The uncertainty of SWAT simulations combined with that of observed climate variables provides a threshold that can be compared with uncertainties arising from using different GCM climate projection data. With the climate variables downscaled from GCMs as inputs, the bias percentage of the simulated streamflow can reach up to 56% compared with that driven by using observed variables as inputs. For the sake of eliminating the impact of errors derived from downscaling, the five-GCM-projected future streamflow was compared with that projected during the historical period. Uncertainties derived from the hydrological model were comparatively small, given that the observed streamflow was successfully reproduced although the peak flow was slightly underestimated. What calls for special attention is the fact that the simulated future streamflow is dependent on the assumption that land use, soil, and topographic slope attributes of the basin were left unchanged throughout the whole simulation period.


## 4.1. Calibration and Validation of Hydrological Model Using Observed Climate Datasets

The parameter uncertainty definitely leads to uncertainties in model simulation. The best-fit parameter values and the sensitivities of model parameters are listed in Table3. The main sources of streamflow uncertainty were the SCS runoff curve number for moisture condition II (CN2), soil evaporation compensation factor (ESCO), groundwater “revap” coefficient (GW_REVAP), deep aquifer percolation fraction (RCHRG_DP), baseflow alpha factor for bank storage (ALPHA_BNK), groundwater delay time (GW_DELAY), moist bulk density (SOL_BD), effective hydraulic conductivity in the main channel alluvium (CH_K2), and threshold depth of water in the shallow aquifer required for the return flow to occur (GWQMN). The P values of these factors are less than 0.05. These results were basically consistent with the study [15].As shown in Figure3, the model reproduced observed river streamflow records fairly well albeit with slight underestimation of the peak flow values. The evaluation metrics of the simulated daily streamflow are listed in Table 6. The NS coefficient values ranged from 0.74 to 0.79 and 0.74 to 0.80 in the calibration period and validation period, respectively, indicating that good performance or very good performance was achieved. With respect to PBIAS values, the streamflows were slightly underestimated at Stations 1, 2, and 3, with PBIAS values ranging from 6.7 to 8.8 in calibration and 0.94 to 10.24 in validation. In contrast, system errors at Station 4 suggested that streamflows were highly overestimated in both calibration and validation periods. Similar results were found in P-factor values. Only Station 4 was below the acceptable level of 0.70, which indicated a relatively large uncertainty in this station, blaming the parameter-adjusting scheme provided by SWAT-CUP. The scheme calibrates the parameters using a single global modification term that scales the initial estimates by values, multiplication or additive terms and thus is not sufficient to cope with the heterogeneity of the precipitation input.Figure 3
Comparisons of daily streamflow simulations (black lines) with observations (gray areas) during calibration (1976–1991) and validation (1992–2000) periods.Table 6
Calibration (1976–1991) and validation (1992–2000) performance metrics for daily runoff at the four monitoring points.
ID
Calibration (1976–1991)
Validation (1992–2000)
R
2
NS coefficient
PBIAS
P-factor
R-factor
R
2
NS coefficient
PBIAS
Station 1
0.75
0.74
8.8
0.89
0.92
0.77
0.74
10.24
Station 2
0.79
0.77
7.6
0.87
0.82
0.82
0.80
0.94
Station 3
0.81
0.79
6.7
0.88
0.94
0.82
0.80
3.77
Station 4
0.80
0.78
−18.3
0.51
0.88
0.84
0.83
−12.50


## 4.2. Calibration and Validation of ANN Downscaling Model

The Pearson correlation coefficient between the observed climate and downscaled variables, including daily maximum temperature, daily minimum temperature, and daily precipitation, from the NCEP/NCAR for the calibration period and validation period was used to evaluate the performance of downscaled climate variables. For daily maximum temperature and daily minimum temperature, the correlation coefficients were all greater than 0.96 in the calibration period and greater than 0.94 in the validation period. The results clearly showed that the magnitudes of the downscaled 30-year daily maximum/minimum temperature values are statistically similar to those of the observations.Annual cycles of mean monthly maximum/minimum temperatures averaged for observations and downscaled from the NCEP/NCAR dataset over the basin during the period of 1971–2000 are shown in Figure4(a). Downscaled maximum temperatures from the NCEP/NCAR basically reproduced the observations with an absolute error less than 0.3°C, except for January, April, and August. In January and August, downscaled NCEP/NCAR mean monthly maximum temperatures were 2.53 and 28.88°C, respectively, which slightly exceeded the observations of 1.95 and 28.54, respectively. In April, the NCEP/NCAR underestimated the mean monthly maximum temperature with its value of 17.20, while the observed value is 17.70°C. Regarding the mean monthly minimum temperature, the downscaled NCEP/NCAR overestimated it in January, February, March, April, August, and September but underestimated it in the other months. The maximum over- and underestimation occurred in February and December with the absolute error values of 0.86 and 0.59, respectively.Figure 4
Annual cycle of mean monthly (a) maximum (solid lines) and minimum (dashed lines) temperatures and (b) precipitation averaged over the Wabash River Watershed for observation and downscaled from the NCEP/NCAR and from five GCMs during the period 1971–2000.
(a)
(b)Correlation coefficients (CCs) between the observed rainfall records and downscaled daily precipitation from the NCEP/NCAR for the calibration period and validation period are shown in Figure5. The values of the CC ranged from 0.63 to 0.79 (with an average of 0.71) in the calibration period except for one node in the western portion with a value of 0.58. The performance of downscaled precipitation in the validation period was slightly worse than that in the calibration period based on the CC values. The CC values ranged from 0.6 to 0.75 (averaged 0.68) during this period with the exception of two nodes located in the eastern and northern portions of the basin with values of 0.47 and 0.50, respectively. Overall, the results showed that the downscaled daily precipitation from the NCEP/NCAR was highly correlated and moderately correlated with observations, respectively.Figure 5
Pearson correlation coefficient between the observations and downscaled maximum temperature (a, b), minimum temperature (c, d), and daily precipitation (e, f) for the calibration period and validation period.
(a)
(b)
(c)
(d)
(e)
(f)The monthly mean precipitation cycles of downscaled precipitation from the NCEP/NCAR and observations over the Wabash River Watershed during the period of 1971–2000 are compared in Figure4(b). Downscaled precipitation captured the bimodal distribution of the annual cycle fairly well, though with different levels of bias and errors in each month. The mean monthly precipitation was overestimated in most of the calendar months, with the exception of October, during which the mean monthly precipitation was slightly underestimated. The discrepancies of monthly mean precipitation between downscaled precipitation and observations were significant in winter months (November, January, and February) and summer months (June, July, August, and September). As a result, the average annual precipitation was 1133.87 mm of NCEP/NCAR downscaled rainfall, which exceeded the observation (1057.31 mm).Five nodes in the precipitation grids of the observation that were evenly distributed over the study basin were selected for frequency distribution analysis and comparison (Figure5). In terms of frequency, the downscaled precipitation from the NCEP/NCAR to a great extent underestimated no/minimal rain events, ranging from 56.82% (at Node 4) to 59.07% (Node 5) compared with the observations that ranged from 70.33% (at Node 1) to 74.72% (Node 5). In addition, the events of high heavy rain and violent rain were also underestimated in this region, and the discrepancies were considerably discernible for violent rain events. The frequency of violent events ranged from 0.05% (at Nodes 3 and 5) to 0.24% (Node 1) for precipitation downscaled from the NCEP/NCAR but ranged from 0.26% (Node 5) to 0.47% (Node 1) for observations. It is notable that the discrepancy of violent events showed an increasing tendency along a south-north gradient. The violent rain events were 48.9% underestimated at Node 1 but 82.8% underestimated at Node 3. On the contrary, the occurrence days of light rain, low moderate rain, high moderate rain, and low heavy rain were significantly overestimated (Figure 6). The results were consistent with the study [46]. Moreover, the frequency distribution of downscaled precipitation records was very similar to that of the NCEP/CFSR over the Wabash River Basin during 1979 to 2000 (Figure 6), and numerous studies have documented the poor performance of the NCEP/CFSR as the driving force of the hydrologic model [47]. Given the high nonlinearity of watershed systems, even small changes in precipitation can cause tremendous variations in watershed response. Good performance of precipitation records in reproducing the observed annual average rainfall does not guarantee good performance in reproducing the streamflow. Thus, we anticipate here that significant differences in streamflow simulations will be observed when the calibrated hydrologic model uses downscaled GCMs’ precipitation as inputs.Figure 6
Frequency distribution of observed daily precipitation values and downscaled precipitation from the NCEP/NCAR and NCEP/CFSR at selected nodes (Figure5(e)).


## 4.3. Hydrological Evaluation of Downscaled Precipitation and Temperature from GCMs in the Historical Period

The calibrated and validated ANN model was then subjected to the GCMs in the historical period to evaluate the performance of each GCM in representing the observations and driving the SWAT hydrologic model. The annual cycle of mean monthly downscaled maximum/minimum temperatures from the five GCMs in the historical period is plotted in Figure4(a). It is obvious that the mean monthly maximum temperature reproduced by the GCMs is consistent with observations except IPSL-CM5A-LR, which to a large extent overestimated the maximum temperature in January, February, May, November, and December with the bias ranging from 1.68 to 3.33. With respect to the downscaled maximum temperature from the other four GCMs, all the absolute errors of mean monthly maximum temperature lied in the range of 0.01 to 1.18°C. Good matching was also witnessed for minimum temperature, with the absolute errors ranging from 0.001 to 1.74°C (Figure 4).Similar to the downscaled precipitation from the NCEP/NCAR, precipitation from the five GCMs basically captured the bimodal distribution of the annual cycle fairly well, albeit with different amplitudes (Figure4(b)). CanESM2 had a tendency to greatly underestimate the mean monthly precipitation in summer season months (May, June, July, and August), leading to an overall underestimation of the annual average rainfall (977.29 mm), which was contrary to NorESM1-M that to some extent overestimated the precipitation in these months and consequently overestimated the annual rainfall with the value of 1091.59 mm. GFDL-ESM2G exhibited improved performance in the second half of the year than that in the first half, except for September. In the first half year, GFDL-ESM2G overestimated the precipitation in January and June but underestimated that in the other months. The total error that resulted from both the positive bias and negative bias nearly cancelled each other, which led to a slight underestimation of the annual average rainfall (1026.56 mm). MIROC-ESM-CHEM produced an annual average rainfall of 1087.85 mm, which slightly exceeded the observations. The maximum error occurred in January, with a positive bias of 33.78 mm. IPSL-CM5A-LR performed best in reproducing the annual cycle of mean monthly precipitation as well as the annual average rainfall. All the bias was less than 10 mm for mean monthly precipitation except January, September, and December with values of 11.42, −11.96, and 19.05, respectively. The annual average rainfall was predicted to be 1051.42 mm, which is comparable to the observation of 1057.31 mm.The downscaled precipitation and maximum/minimum temperature from the five GCMs were then used to drive the calibrated SWAT model, respectively. The SWAT simulated annual streamflow produced by the five GCM projections indicated a wide range of annual runoff and flowrates (Figure7). As expected, all the GCMs except CanESM2 in the 8th percentile to a large extent overestimated the annual streamflow (Figure 7(a)), which was mainly blamed on the overestimation of the daily streamflow less than 2000 m3/s (Figure 7(b)). As shown in Table 7, the annual streamflow under the 90th percentile was substantially overestimated, while the ensemble mean of simulations with climate variables downscaled from five GCMs as inputs fit well for streamflows observed from 2000 to 4000 m3/s (approximately 90th percentile to 99th percentile). Beyond that, the streamflow rate was slightly underestimated. Overall, the downscaled precipitation and maximum/minimum temperature from the five GCMs as driving forces of the established hydrological model resulted in relatively high overestimation of the observed streamflow throughout the simulation period. The large temporal variation of precipitation (Figure 6) is present and thus leads to streamflow overestimation by impacting the process of runoff generation when using the hydrological model.Figure 7
Q-Q plots of (a) annual runoff and (b) daily streamflow for Station 3.
(a)
(b)Table 7
Annual streamflow statistics simulated using the five GCMs for the historical period and observation.
Annual streamflow (107 m3)
1
2
3
4
Observation
2731.61
1812.68
644.33
274.60
CanESM2
3416.56
2482.26
675.49
258.82
GFDL-ESM2G
3755.29
2855.95
632.29
267.05
IPSL-CM5A-LR
3766.90
2830.17
682.66
254.07
MIROC-ESM-CHEM
4040.99
3108.98
681.87
250.14
NorESM1-M
4262.30
3265.12
709.80
287.38
Ensemble
3848.41
2908.50
676.42
263.49
1: average annual runoff; 2: annual streamflow under the 90th percentile; 3: 90th–99th percentile; 4: 99th exceedance percentile.


## 4.4. Impact of Climate Change on Streamflow

The downscaled future precipitation and temperature for the five GCMs under the three RCP scenarios were inputted to the calibrated SWAT, separately, to investigate corresponding streamflow responses. Figure8 presents the ensemble mean of average annual runoff of Station 3 for the historical, midcentury, and late-century timeframes under three RCP scenarios, from which some interesting implications can be drawn. All projected future changes in annual runoff under RCP scenarios showed a decreasing trend to different degrees. In the midcentury period, the multimodel ensemble means were projected to decrease by 13.00%, 17.59%, and 6.91% under RCP2.6, RCP4.5, and RCP8.5 scenarios, respectively, compared with the historical period. Considerably, more significant decreases were witnessed in the late-century period, with the decrease rates of 25.29%, 27.61%, and 15.04% for the three scenarios (Table 8).Figure 8
Annual cycle of the monthly mean streamflow for the ensemble mean of SWAT simulations in Station 3 with downscaled climate variables from five GCMs as inputs over the period of 1976–2000.Table 8
Statistics of annual runoff under three scenarios for the two future periods and their changes compared to that under the historical scenario.
Scenarios
Annual streamflow (107 m3) and changes
RCP2.6
RCP4.5
RCP8.5
Midcentury
3354.64 (−13.00%)
3177.81 (−17.59%)
3589.39 (−6.91%)
Late century
2880.90(−25.29%)
2791.31 (−27.61%)
3276.05 (−15.04%)The annual cycle of the simulated multimodel ensemble mean of monthly averaged streamflows is also shown in Figure8. Under RCP2.6 and RCP4.5 scenarios, mean monthly streamflows in all calendar months decreased to varying degrees. The decline rates ranged from 0.40% to 26.03% for RCP2.6 scenarios in the midcentury and from 5.04% to 39.00% in the late century. For RCP4.5 scenarios, however, the ranges were 6.68%–34.12% and 7.51%–42.74%, respectively. The biggest decreases all occurred in November to March. Generally, the RCP2.6 and RCP4.5 scenarios projected very similar annual runoff and monthly runoff change, which were very different from those of RCP8.5. In the midcentury, the mean monthly streamflows under the RCP8.5 scenario were slightly decreased in most months, except for April and June, when the values were increased to be comparable with that of the historical scenario. In the late century, the mean monthly streamflows also showed an increasing trend in April and May with a decreasing trend in other months. Similar to RCP2.6 and RCP4.5, the biggest decline under the RCP8.5 scenario appeared in November to March.The projected annual runoff change was neither positively nor negatively related to the CO2 emission scenarios (i.e., RCP scenarios), which may be attributed to the changes in both of the rainfall pattern projected towards the midcentury and late century and an increase in the evaporation demand in terms of the climate warming. To distinguish the impacts of the projected change in precipitation and temperature on the simulated streamflow, the annual cycle of multimodel ensemble means of precipitation and maximum/minimum temperature under three scenarios for the two future periods is shown in Figure 9. The multimodel ensemble means of annual precipitation under RCP2.6 and RCP4.5 scenarios for the two future periods were 1023.37 (RCP26-midcentury), 1079.91 (RCP26-late century), 999.27 (RCP45-midcentury), and 1077.08 (RCP45-late century), respectively, which were very close to those of the historical period (1046.94 mm) and the observation (1057.31 mm). No significant change in precipitation was projected under RCP8.5 for the midcentury period with the annual average rainfall of 1076.87 mm; however, the precipitation under RCP8.5 for the late-century period was projected to have a 14.23% increase compared with the rainfall in the historical period. With respect to the maximum temperature, RCP2.6 and RCP4.5 projected an increase of 1.91°C (RCP26-midcentury), 1.99°C (RCP26-late century), 2.12°C (RCP45-midcentury), and 3.15°C (RCP45-late century), while RCP8.5 projected an increase of 2.32°C for the midcentury period and a significant increase of 5.31°C for the late-century period. Similar results were found for minimum temperature with increased values of 1.64°C (RCP26-midcentury), 1.80°C (RCP26-late century), 1.86°C (RCP45-midcentury), 2.88°C (RCP45-late century), 2.10°C (RCP85-midcentury), and 4.92°C (RCP85-late century). All of the above projection results are consistent with those of the IPCC report over this basin [2]. Thus, we can assert that climate warming plays an important role for RCP2.6 and RCP4.5, especially for the late-century period. For the RCP8.5 cases, during the midcentury period, the projected runoff slightly decreased though the changes in both the rainfall and temperature, very similar to the RCP26-midcentury case. We predicted that the rainfall pattern would be very different from the present pattern. During the late-century period, RCP8.5 projected a decrease of monthly streamflow with the exception of April, May, and August. Both precipitation and temperature are primary factors affecting the decrease in the streamflow.Figure 9
Annual cycle of mean monthly (a) precipitation and (b) maximum (solid lines) and minimum (dashed lines) temperatures averaged over the Wabash River Watershed for observation and multimodel ensembles for the historical period, midcentury, and late century under RCP scenarios.
(a)
(b)Regarding the mean monthly precipitation, an interesting finding was that all the cases, except for RCP4.5 and RCP8.5 during the late-century period, witnessed a decrease of rainfall in the summer season (June, July, and August), during which crops grow. Moreover, all the cases with the exception of RCP2.6 and RCP4.5 for the midcentury period indicated that the rainfall tends to increase in the winter and spring seasons.Overall, all scenarios projected a decrease of the streamflow with different amplitudes. These results were basically consistent with the study by Chien et al. [15]. Most of RCP2.6 and RCP4.5 scenarios projected an increase in precipitation in November to February, while the simulated streamflows will decrease in these months. These results suggest that temperature is a very important factor in the hydrological regime in this area. In summer months, the decrease in the streamflow was not serious with respect to the decrease of precipitation and increase of temperature. Regarding the RCP8.5 scenario, especially for the late-century period, large uncertainties derived from GCMs and the ANN downscaling method were witnessed. Therefore, further quantifying investigation of the uncertainties is necessary for reliable assessment of climate change impacts.


## 4.5. Uncertainty Analysis

Studies on climate change impact are always accompanied with uncertainties derived from GCMs, downscaling methods, and hydrological modeling. In this study, five GCMs were selected and used to drive the hydrological model. In view of the wide range of projected climate variables, especially precipitation, the bias was still large although bias correction was performed. With respect to the annual average rainfall, the absolute error was estimated ranging from 0.01 to 3.24 mm. The downscaling process is another important source of uncertainties. For this study, the frequency distribution of downscaled rainfall from GCMs was very similar to that of the NCEP/CFSR and deviated significantly from the observation. The mismatch of frequency distribution of precipitation led to a dramatic overestimation of the projected streamflow. The uncertainty of SWAT simulations combined with that of observed climate variables provides a threshold that can be compared with uncertainties arising from using different GCM climate projection data. With the climate variables downscaled from GCMs as inputs, the bias percentage of the simulated streamflow can reach up to 56% compared with that driven by using observed variables as inputs. For the sake of eliminating the impact of errors derived from downscaling, the five-GCM-projected future streamflow was compared with that projected during the historical period. Uncertainties derived from the hydrological model were comparatively small, given that the observed streamflow was successfully reproduced although the peak flow was slightly underestimated. What calls for special attention is the fact that the simulated future streamflow is dependent on the assumption that land use, soil, and topographic slope attributes of the basin were left unchanged throughout the whole simulation period.


## 5. Conclusions

In this study, we applied climate change simulations generated by five GCMs that contributed to CMIP5 for simulating the midcentury and late-century behavior of streamflow conditions over the Wabash River Basin. Large-scale climate variables simulated by the GCMs are first bias-corrected by adopting a quantile mapping method. Then, an ANN method was used for statistical downscaling to bridge the scale mismatch between the large-scale climate variables and the regional hydrological model. The potential impacts of climate change on runoff were quantified by driving the SWAT hydrological model with climate information during the historical and future periods under RCP2.6, RCP4.5, and RCP8.5 scenarios.Sensitivity analysis of the SWAT model indicated that CN2, ESCO, GW_REVAP, RCHRG_DP, ALPHA_BNK, GW_DELAY, SOL_BD, CH_K2, and GWQMN are the most sensitive parameters during the calibration process. Based on NS coefficient values, the simulations of the calibrated SWAT model fit the observed streamflow well in both the calibration and validation periods. TheP factor values suggest that the simulated predictions were considered acceptable with the exception of that for Station 4. The observed streamflow was significantly overestimated in this station. Overall, the calibrated SWAT model basically captured the hydrological process characteristics but underestimated the peak flow values to some extent.The downscaled maximum temperature and minimum temperature from five GCMs using the calibrated ANN model reproduced the observation very well based on the Pearson correlation coefficient. With respect to precipitation, the downscaled annual average precipitation from five GCMs compared fairly well with the observation albeit with different amplitudes. In addition, the annual cycle of mean monthly precipitation downscaled from GCMs also matched the observation well. However, downscaled precipitation from the NCEP/NCAR underestimated the events of no/minimal rain, high heavy rain, and violent rain, while it overestimated the occurrence days of light rain, low moderate rain, high moderate rain, and low heavy rain.The downscaled future climate variables were used to drive the calibrated SWAT model, and the results showed the potential for a significant decrease of the streamflow. The projected streamflow could decrease by 13.00%, 17.59%, and 6.91% in the midcentury period and 25.29%, 27.61%, and 15.04% in the late-century period under RCP2.6, RCP4.5, and RCP8.5 scenarios, respectively. Changes in temperature play a dominant role in shaping the streamflow simulations under RCP2.6 and RCP4.5 scenarios, while both of the changes in temperature and precipitation are primary factors that affect the decrease in the streamflow. Given the unexpected streamflow simulation under the RCP8.5 scenario for the midcentury period, it is predicted that the rainfall pattern would be radically shifted [48].



---

*Source: 9698423-2020-02-10.xml*