Scenario development:
Two projection period, both under RCP4.5 (moderate) and RCP6.0 (severe) were presented, mid-century (2040-2069) and late-century (2070-2099). The results of hydrological simulation were shown in monthly, seasonal, and annual time scales. The seasons were defined as DJF (winter: December, January, February), JJA (summer: June, July, August), MAM (spring: March, April, May), and SON (fall: September, October, November).
SWAT incorporates CO2 to account for its impact on plant water requirements and on level of the potential evapotranspiration (PET) (Neitsch, Arnold et al. 2011). It takes the CO2concentration amount as a single input value for each subbasin. The CO2 concentration values for the historical and future projections are shown in Table 5. The values are derived from Meinshausen, Smith et al. (2011).
Projected population for the conterminous US indicates significant increase in demand for food, energy, and urban development (Sohl, Sayler et al. 2014). From 2001 to 2011, the Southeast region (AL, AR, FL, GA, KY, LA, MS, NC, SC, TN, and VA) has lost more than 100 and 1400 sq2 agricultural land and forest, respectively, and gained 600 sq2 developed land cover (Sleeter, Loveland et al. 2018). For UCS and Pea and Yellow River Subbasin, farming land decreased 27.21% and urban area increased 42.55% from 1992 to 2011 (Hinson, Rogers et al. 2015). Regarding land use condition in the future, there have been few studies (national US and global scale) based on different scenarios including Special Report on Emission Scenarios (SRES) (96) , RCP, and Shared Socioeconomic Pathways (SSP) (95) of Intergovernmental Panel for Climate Change (IPCC) (Nakicenovic, Alcamo et al. 2000; Wear 2011; Sohl, Sayler et al. 2014; Sohl, Wimberly et al. 2016; Riahi, Van Vuuren et al. 2017; Sleeter, Loveland et al. 2018). Sohl et al. (2014) using different land use forecasting model, has predicted 22.9% to 61% increase in urban land cover for conterminous US by the year 2050. They projected noticeable loss of natural covers which was due to expansion of anthropogenic land uses. The fourth National Climate Assessment reported 50% and 80% increase in urban land use allocation by 2100 under SSP2 and SSP5 respectively with 2010 land use condition as the baseline (Sleeter, Loveland et al. 2018). To account for these changes, we obtained and analyzed projected land covers for each decade till 2100 from Fourth National Climate Assessment dataset through Global Change Explorer (GCX) ((GCX) 2020). These maps are based on SSP scenarios with 19 land cover classes (Bierwagen, Theobald et al. 2010; (EnvironmentalProtectionAgency) 2017). One possible caveat though, is that there is not much agreement between different forecasting models and they appear to be at the beginning stage of development (Sohl, Wimberly et al. 2016). Future projection results for all models in this study were presented in monthly, seasonal, and annual average and compared to the historical result to analyze the future hydrological condition within the UCS. Then we did the same comparison for SWAT simulation, mainly discharge. Corresponding land use projections to projection periods (mid-century and end-century) were used in SWAT modeling to simulate the discharge and evapotranspiration (ET). Finally, we used box plots to show the changes of the climate and hydrological variables.
SWAT: For this study we used Soil & Water Assessment Tool (SWAT). SWAT is assemblages of mathematical equations representing different parts of hydrological cycle including movement, fate, and transport of water, sediments and nutrients in and on soil, through groundwater, and in river streams and reservoirs (Arnold, Srinivasan et al. 1998; ASABE Jun. 2017). The development of the model started in early 1990s and it has been evolving by USDA (United States Department of Agriculture) agricultural Research Service (ARS) (Gassman, Reyes et al. 2007; Arnold, Moriasi et al. 2012). It is a process based and semi-distributed continuous-time river basin scale model (Arnold, Kiniry et al. 2011; Arnold, Moriasi et al. 2012). It has been written in Fortran language including more than 310 subroutines representing different parts of the hydrological and bio-geochemical processes (Arnold and Fohrer 2005; Arnold, Kiniry et al. 2011). It was originally developed to evaluate water resources management and Nonpoint Source (NPS) pollution in large river basin (Arnold, Moriasi et al. 2012). It has proven to be effective for its purposes and computationally efficient and can be used for long term continuous simulation including climate change impact studies (Gassman, Reyes et al. 2007; Arnold, Moriasi et al. 2012). It operates on a daily time step and outputs daily, monthly, and yearly results (Gassman, Reyes et al. 2007; Arnold, Moriasi et al. 2012). SWAT splits a watershed into subwatersheds that are further split into hydrologic response units (HRUs) (Gassman, Reyes et al. 2007; Arnold, Moriasi et al. 2012). HRUs are nonspatial units and unique combination of homogeneous land use, soil, slope and management characteristics (Gassman, Reyes et al. 2007; Arnold, Moriasi et al. 2012). This gives SWAT the capability to model surface runoff, infiltration, soil water movement, ET, in-stream transformations, sediment movement, canopy interception, plant uptake, and nutrients circulation including biogeochemical processes at HRU level (Neitsch, Arnold et al. 2011). Main components of a SWAT model for a given watershed are weather, hydrology, erosion/sedimentation, plant growth, nutrients, pesticides, agricultural management, stream routing and pond/reservoir routing (Arnold and Fohrer 2005). Simulation in SWAT has two parts, land phase and routing phase; in land phase, the amount of water, sediment, nutrient, and pesticide loadings are regulated into the main channel in each subbasin, and in the routing phase, in- stream processes including water movement, sediment transport and the nutrients loading are simulated (Neitsch, Arnold et al. 2011; Arnold, Moriasi et al. 2012). In SWAT, Water balance is the base of all the processes and the hydrological cycle is climate driven, thus, SWAT requires precipitation, minimum and maximum temperature, solar radiation, relative humidity, and wind speed in daily time scale (Arnold, Moriasi et al. 2012).
SWAT uses equation (1) to simulate water balance (Arnold, Srinivasan et al. 1998; Neitsch, Arnold et al. 2011):
\begin{equation} SW=\text{SW}_{0}+\ \sum_{i=1}^{t}{(R_{\text{day}}-Q_{\text{surf}}-}E_{a}-w_{\text{seep}}-Q_{\text{gw}})\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1)\nonumber \\ \end{equation}
where SW and \(\text{SW}_{0}\) are soil water content for beginning and end of the model, respectively. t (day) is time.\(R_{\text{day}}\)is rainfall; \(Q_{\text{surf}}\)is surface runoff;\(E_{a}\)is evapotranspiration; \(w_{\text{seep}}\)is percolation to vadose zone, and \(Q_{\text{gw}}\ \)is return flow amount, and all variables are in mm (Arnold, Srinivasan et al. 1998; Neitsch, Arnold et al. 2011).
Water yield as part of subbasins blue water is the amount of water after leaving HRUs and entering the main channel is calculated with equation (2) (Neitsch, Arnold et al. 2011; Arnold, Kiniry et al. 2013; Veettil and Mishra 2016):
\begin{equation} WYLD=Q_{\text{surf}}+Q_{\text{lat}}+Q_{\text{gw}}-tloss\ -Pond\ abstractions\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2)\nonumber \\ \end{equation}
where WYLD is the amount of water yield, \(Q_{\text{surf}}\)is surface runoff, \(Q_{\text{gw}}\ \)is return flow amount, \(Q_{\text{lat}}\) is the amount of lateral flow, tloss is transmission losses, and the abstracted water from the pond; all variables are in mm (Arnold, Kiniry et al. 2013; Chanapathi, Thatikonda et al. 2018; Pandey, Khare et al. 2019).
In SWAT surface runoff can be estimated in two ways: SCS (Soil Conservation Service) runoff curve number method (USDA-SCS, 1972) and the Green & Ampt infiltration method (1911) (Neitsch, Arnold et al. 2011)(45). We used the former. Equation (3) is SCS runoff curve number method (Arnold, Srinivasan et al. 1998; Neitsch, Arnold et al. 2011):
\begin{equation} \text{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }Q_{\text{surf}}=\ \frac{\left(R_{\text{day}}-0.2\ S\right)^{2}}{\left(R_{\text{day}}+0.8\ S\right)}\text{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }R_{\text{day}}>0.2S\ \ \ \ \ \ \ \ \ \ \ \ \left(3\right)\text{\ \ \ }\nonumber \\ \end{equation}\begin{equation} Q_{\text{surf}}=0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ R_{\text{day}}\leq 0.2S\nonumber \\ \end{equation}
where \(Q_{\text{surf}}\)is surface runoff; \(R_{\text{day}}\)is rainfall; and S is retention parameter. 0.2S is estimated as the initial abstraction including surface storage (Neitsch, Arnold et al. 2011). Retention parameter varies through the watershed and time owing to changes in soil, land use and management (Neitsch, Arnold et al. 2011). S is estimated as follows:
\begin{equation} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ S=254\ (\frac{100}{\text{CN}}-10)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left(4\right)\text{\ \ \ }\nonumber \\ \end{equation}
where CN is the curve number which is adjusted for different soil moisture level and slope (Neitsch, Arnold et al. 2011).
For surface runoff equation 3 was used. For flow routing the variable storage coefficient method were used (Williams 1969; Neitsch, Arnold et al. 2011). Since our modeling required simulation of CO2 climate change effects, the Penman-Monteith method was used for calculation of potential evapotranspiration(Monteith 1965; Allen 1986; Allen, Jensen et al. 1989). Actual evapotranspiration (AET) the was calculated by procedure established by Richtie (1972) (Ritchie 1972). The UCS was delineated into 54 subbasin and 1821 HRUs. ‘SWAT2012 rev64’ version was used to perform the modeling.
Calibration :
calibration can be manually done or through a combination of manual and auto calibration procedures (Moriasi, Wilson et al. 2012). Our approach to calibrate and validate was the later through a split-sample strategy (Moriasi, Wilson et al. 2012). To evaluate the performance of the model we have first carried out a sensitivity analysis (SA) manually and then using SWAT Calibration and Uncertainty Procedures (SWAT-CUP) to filter out insensitive parameters to reduce the computational workload of the calibration (Gupta, Sorooshian et al. 1999; Saltelli, Tarantola et al. 2004; Abbaspour, Vejdani et al. 2007; Ghoraba 2015). Sensitivity analysis is to estimate how much model outputs change with respect to each model parameter (input) change (Saltelli, Tarantola et al. 2004; Arnold, Moriasi et al. 2012). First a set the parameters were selected according to UCS hydrologic characteristics and the literature (Abbaspour, Yang et al. 2007; Joh, Lee et al. 2011; Sudheer, Lakshmi et al. 2011; Arnold, Moriasi et al. 2012; Abbaspour, Rouholahnejad et al. 2015; Osei, Amekudzi et al. 2019; Qiu, Shen et al. 2019). Then using one-factor-at-a time sensitivity analysis initial parametrization was carried out and parameters were optimized and their initial ranges were predicted (Morris 1991; Green and Van Griensven 2008; Abbaspour, Rouholahnejad et al. 2015). After a set of manual calibration using first set of the parameters, we used SWAT-CUP to modify the selected parameters and perform sensitivity and uncertainty analysis (Arnold, Moriasi et al. 2012). We used the Sequential Uncertainty Fitting version algorithm (SUFI2) within SWAT-CUP (Abbaspour, Johnson et al. 2004; Abbaspour, Yang et al. 2007). The SUFI2 is based on the invers modeling and is to estimate parameters using observed data (Abbaspour, Johnson et al. 2004; Abbaspour, Yang et al. 2007). In other words, it uses initial large parameter uncertainty and through steps, decrease the uncertainty until the uncertainty range falls within a range/band called 95% Prediction Uncertainty (95PPU) (Abbaspour, Johnson et al. 2004; Abbaspour 2015). SUFI2 uses a global search approach to carry out optimization and uncertainty analysis and it can handle many parameters (Abbaspour, Johnson et al. 2004; Abbaspour 2015). For accuracy quantification of the model, our objective function includes Nash–Sutcliffe Efficiency (NSE), Coefficient of Determination (R2), Percent Bias (PBAIS), RMSE-observations standard deviation ratio (RSR) (Nash and Sutcliffe 1970; Gupta, Sorooshian et al. 1999; Legates and McCabe Jr 1999; Moriasi, Arnold et al. 2007). The metrics for satisfactory thresholds were selected based on the literature (Santhi, Arnold et al. 2001; Moriasi, Arnold et al. 2007; Abbaspour, Rouholahnejad et al. 2015). Table 6 shows the objective function and the thresholds of the metrics and final results for calibration and validation period. Following are the formulas used for these metrics.
\begin{equation} NSE=1-\ \left[\frac{\sum_{i=1}^{n}\left(Y_{i}^{\text{obs}}\ -\ Y_{i}^{\text{sim}}\right)^{2}}{\sum_{i=1}^{n}\left(Y_{i}^{\text{obs}}\ -\ Y^{\text{mean}}\right)^{2}}\right]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ,\ \ \ \ \ \ \ \ \ \ \ \ -\infty\ <NSE\ \leq\ 1\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left(5\right)\nonumber \\ \end{equation}\begin{equation} R^{2}\ =\ \left[\frac{\sum_{i=1}^{n}\left(Y_{i}^{\text{obs}}-\ Y_{i}^{\text{mean}}\right)\left(Y_{i}^{\text{sim}}\ -\ Y^{\text{simmean}}\ \right)}{\left[\sum_{i=1}^{n}\left(Y_{i}^{\text{obs}}-\ Y_{i}^{\text{mean}}\right)^{2}\right]^{0.5}\left[\sum_{i=1}^{n}\left(Y_{i}^{\text{sim}}\ -\ Y^{\text{simmean}}\ \right)^{2}\right]^{0.5}}\right]^{2}\ \ ,\ \ \ \ 0\leq R^{2}\ \leq\ 1\ \ \ \ \ \ \ \ \ \ \ \left(6\right)\text{\ \ }\nonumber \\ \end{equation}\begin{equation} PBAIS\ =\ \left[\frac{\sum_{i=1}^{n}{\left(Y_{i}^{\text{obs}}\ -\ Y_{i}^{\text{sim}}\right)\ \times\ \left(100\right)}}{\sum_{i=1}^{n}\left(Y_{i}^{\text{obs}}\right)}\right]\ \ \ \ \ \ \ \ \ \ ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ -\infty\ <\ PBAIS\ \ <\ +\infty\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left(7\right)\text{\ \ \ \ \ \ }\nonumber \\ \end{equation}\begin{equation} RSR\ =\ \frac{\left[\sqrt{\sum_{i=1}^{n}\left(Y_{i}^{\text{obs}}-\ Y_{i}^{\text{sim}}\right)^{2}}\right]}{\left[\sqrt{\sum_{i=1}^{n}\left(Y_{i}^{\text{obs}}-\ Y^{\text{mean}}\right)^{2}}\right]}\ \ \ \ \ \ \ \ \ \ \ \ \ \ ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\ \leq\ \ RSR\ \ <\ +\infty\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left(8\right)\text{\ \ \ \ }\nonumber \\ \end{equation}
where \(Y_{i}^{\text{obs}}\) is the i th measured stream flow;\(Y_{i}^{\text{sim}}\) is i th simulated stream flow;\(Y^{\text{mean}}\) is the mean of observed stream flow data;\(Y^{\text{simmean}}\) is the mean of simulated data, n the total number of observations.
NSE value between 0.0 and 1 is considered acceptable with 1 being the optimal value indicating the plot of observed versus simulated fits perfectly (Nash and Sutcliffe 1970; Moriasi, Arnold et al. 2007). For R2, the higher the value the lesser the error variance and values greater than 0.5 are considered to be acceptable (Santhi, Arnold et al. 2001; Moriasi, Arnold et al. 2007). Negative PBIAS means overestimation and positive PBIAS means underestimation, with 0 being the optimal value. RSR ranges from 0.0 to a positive large number, with 0.0 being the optimal condition meaning the RMSE is zero. The lower RSR, the lower residual variation, indicating better model performance (Moriasi, Arnold et al. 2007). NSE, R2, and RSR are unitless and PBAIS has the unit of the constituent being evaluated which for our case is cms (m3/s). The global sensitivity analysis (where all parameters are allowed to change through analysis) were carried out through SWAT_CUP to prioritize the most responsive parameters and remove parameters with smaller sensitivity from further sampling (Abbaspour, Yang et al. 2007). Through the SA, a multiple regression analysis was used to determine the parameter sensitivity statistics (t-stat and p-value) (Abbaspour 2015). t-stat and p-value were used to describe the relative significance and significance of sensitivity, respectively (Chanapathi, Thatikonda et al. 2018). Sensitive parameters correspond to larger absolute t-stat values among the parameters and to smaller p-values (close to 0; commonly accepted threshold is 0.05) (Abbaspour 2015). The SA results are shown in Table 7. Based on the sensitivity analysis the following parameter were identified as responsive: SOL_AWC (Available water capacity of the soil layer), RCHRG_DP (Deep aquifer percolation fraction), CH_K2 (Effective hydraulic conductivity in main channel alluvium), SLSUBBSN (Average Slope Length), ESCO (Soil evaporation compensation coefficient), GWQMN (Threshold water level in the shallow aquifer for the base flow), and ALPHA_BF (Baseflow alpha factor). Table 8 shows final parameters and their fitted values. It is worth to mention that the SA result is the prediction of the average changes of the objective function being produced by changes in a given parameter as other parameters are also changing (Khalid, Ali et al. 2016). p-factor and r-factor were also considered for SUFI2 performance evaluation as well as measuring the goodness of calibration (Abbaspour, Johnson et al. 2004; Abbaspour, Yang et al. 2007). p-factor is the percentage of observations covered by the 95PPU and r-factor is the ratio of the 95PPU average thickness and the standard deviation of the observations (Abbaspour 2015; Abbaspour, Rouholahnejad et al. 2015). The extent from 2.5% to 97.5% of the cumulative distribution of the simulated variable resulting from the Latin hypercube sampling is 95PPU (Abbaspour, Johnson et al. 2004; Abbaspour, Yang et al. 2007; Abbaspour, Rouholahnejad et al. 2015). p-factor ranges from 0 to 1; p-factor greater than 0.7 means acceptable goodness of fit; r-factor varies between 0 and infinity; here r-factor less than 1.5 was considered satisfactory (Abbaspour, Johnson et al. 2004; Abbaspour, Yang et al. 2007; Abbaspour 2015; Abbaspour, Rouholahnejad et al. 2015). Table 6 shows the goodness of fit metrics. The selected parameters were used to calibrate SWAT at two USGS site (Newton (USGS2361000) and Bellwood (USGS2361500)) for subbasins 29 and 52, respectively. Observed data from the later site is used for SA and uncertainty analysis. The number of simulations was 450 with 4 iteration. SWAT was performed for period from 1998 to 2013. From 1998 to 2000 was considered as warm up; from 2001 to 2010 was the calibration period and from 2010 to 2013 was the validation period. Given the final values for the model performance metrics (Table 6) and the accepted thresholds, it was determined SWAT stream flow estimation for the UCS was efficient. Figure 2 shows the calibration result for both calibration (2001-2010) and validation period (2011-2013).