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).