Development of a Distributed Hydrologic Model for a Region with Fragipan Soils to Study Impacts of Climate on Soil Moisture : A Case Study on the Obion River Watershed in West Tennessee

Previous land surface modeling efforts to predict and understand water budgets in the U.S. Southeast for soil water management have struggled to characterize parts of the region due to an extensive presence of fragipan soils for which current calibration approaches are not adept at handling. This study presents a physically based approach for calibrating fragipan-dominated regions based on the “effective” soil moisture capacity concept, which accounts for the dynamic perched saturation zone effects created by the low hydraulic capacities of the fragipan layers. The approach is applied to the Variable Infiltration Capacity model to develop a hydrologic model of the Obion River Watershed (ORW), TN, which has extensive fragipan coverage. Model calibration was performed using observed streamflow data, as well as evapotranspiration and soil moisture data, to ensure correct partitioning of surface and subsurface fluxes. Estimated Nash-Sutcliffe coefficients for the various sub-drainage areas within ORW were all greater than 0.65, indicating good model performance. The model results suggest that ORW has a high responsivity and high resilience. Despite forecasted temperature increases, the simulation results suggest that water budget trends in the ORW are unlikely to change significantly in the near future up to 2050 due to sufficient precipitation amounts.


Introduction
The soil water balance (Figure 1) is a prime control of the structure and function of intensively managed agroecosystems, driving the dynamics of both soil microbial communities, and crops [1].Stress resulting from soil moisture deficits nearly halts nutrient cycling [2], while anoxic conditions that develop if the soil is too wet can alter microbial decomposition processes [3].Both conditions can negatively impact the physical and chemical make-up of the soil leading to crop yield gaps.Thus, proper management of the soil water balance is essential to maintain both the productivity and sustainability of any agroecosystem [4].
Even in historically "water-rich" regions, like the U.S. Southeast where there are abundant fresh water supplies and an extensive water management infrastructure, there is a risk for soil water deficits considering increasing demands for clean water e.g., [5] and escalating occurrences of climate extremes, e.g., [6].Recent studies project decreases in precipitation amounts during the summer months in the Southeast, as well as a potential for more frequent, intense storm events; however, there is high uncertainty associated with these forecasts [6].Understanding the water budget and predicting potential changes in the near future are crucial for optimizing soil water management scenarios and developing more effective and equitable allocation procedures for crop production and domestic water use.
Geosciences 2018, 8, x FOR PEER REVIEW 2 of 25 proper management of the soil water balance is essential to maintain both the productivity and sustainability of any agroecosystem [4].Even in historically "water-rich" regions, like the U.S. Southeast where there are abundant fresh water supplies and an extensive water management infrastructure, there is a risk for soil water deficits considering increasing demands for clean water e.g., [5] and escalating occurrences of climate extremes, e.g., [6].Recent studies project decreases in precipitation amounts during the summer months in the Southeast, as well as a potential for more frequent, intense storm events; however, there is high uncertainty associated with these forecasts [6].Understanding the water budget and predicting potential changes in the near future are crucial for optimizing soil water management scenarios and developing more effective and equitable allocation procedures for crop production and domestic water use.Physically based hydrologic models are useful for simulating complex physical landatmosphere interactions and may be applied to predict water budgets for current climate conditions, as well as for projected conditions [9].An example of such a model is the Variable Infiltration Capacity (VIC) model, which has been used successfully to simulate water budgets across the U.S. and other regions around the globe [10][11][12][13].VIC is a physically based semi-distributed, macroscale hydrologic model used to quantify water and energy balances at a daily or sub-daily time step [14,15].It has been widely used for different time and basin scales [10][11][12][13] to monitor and forecast water budget changes.VIC is responsive to both seasonal and event dynamics, and can consider spatial heterogeneities in climate and land use, making it a particularly useful model for predictions related to climate and land use changes, e.g., [16,17].
Past hydrological modeling studies in the conterminous U.S. have shown bias in predicting streamflows in the Southeast.This bias has been attributed to calibration issues [18,19].Since other regions with similar climates have not shown a similar bias, landscape characteristics in the Southeast are believed to be the cause of the bias.Though general physical ranges for model calibration parameters have been provided in the literature, there are certain parameters and conditions that can be challenging to characterize.For example, the b-parameter of the VIC model is an empirical parameter used to describe the distribution of the infiltration capacity of the landscape, where the infiltration capacity is defined here as the total volumetric capacity of a vertical soil column to hold water [10].The b-parameter dictates the hydrologic behavior of the study area in terms of rainfallrunoff-infiltration partitioning.While a few studies have demonstrated the determination of the bparameter directly from observed data of soil depth and porosity, e.g., [10,20], these studies have made the implicit assumption of a "static" distribution of soil moisture capacity, which is not always true.This is illustrated by Huang et al. [20] using the soil moisture capacity on a typical hillslope.Despite the soil thickness, and consequently total soil moisture capacity, generally increasing downslope, the "effective" moisture capacity is lowest at the footslope because these soils tend to be Physically based hydrologic models are useful for simulating complex physical land-atmosphere interactions and may be applied to predict water budgets for current climate conditions, as well as for projected conditions [9].An example of such a model is the Variable Infiltration Capacity (VIC) model, which has been used successfully to simulate water budgets across the U.S. and other regions around the globe [10][11][12][13].VIC is a physically based semi-distributed, macroscale hydrologic model used to quantify water and energy balances at a daily or sub-daily time step [14,15].It has been widely used for different time and basin scales [10][11][12][13] to monitor and forecast water budget changes.VIC is responsive to both seasonal and event dynamics, and can consider spatial heterogeneities in climate and land use, making it a particularly useful model for predictions related to climate and land use changes, e.g., [16,17].
Past hydrological modeling studies in the conterminous U.S. have shown bias in predicting streamflows in the Southeast.This bias has been attributed to calibration issues [18,19].Since other regions with similar climates have not shown a similar bias, landscape characteristics in the Southeast are believed to be the cause of the bias.Though general physical ranges for model calibration parameters have been provided in the literature, there are certain parameters and conditions that can be challenging to characterize.For example, the b-parameter of the VIC model is an empirical parameter used to describe the distribution of the infiltration capacity of the landscape, where the infiltration capacity is defined here as the total volumetric capacity of a vertical soil column to hold water [10].The b-parameter dictates the hydrologic behavior of the study area in terms of rainfall-runoff-infiltration partitioning.While a few studies have demonstrated the determination of the b-parameter directly from observed data of soil depth and porosity, e.g., [10,20], these studies have made the implicit assumption of a "static" distribution of soil moisture capacity, which is not always true.This is illustrated by Huang et al. [20] using the soil moisture capacity on a typical hillslope.Despite the soil thickness, and consequently total soil moisture capacity, generally increasing downslope, the "effective" moisture capacity is lowest at the footslope because these soils tend to be wetter due to higher groundwater levels.Consequently, footslopes tend to be saturated more frequently than ridges.This saturation behavior cannot be accounted for by considering only the total volume that the soil can hold.One must also consider additional "dynamic" factors or processes that influence moisture-holding capacity.
Another example where current calibration approaches may be inappropriate is for regions with fragipan soils [21].Fragipan soils are dense pans that are seemingly cemented when dry, brittle when moist, and slake when submerged in water containing redox features [22,23].Observations suggest that they cover nearly 17 million ha in the Eastern and Southeastern United States [24].Subsurface layers containing fragipan soils are water-restricting horizons that can lead to perched saturation zones, thereby affecting the effective moisture capacity and controlling both runoff and baseflow.The effective moisture capacity in this case is regulated by the low saturated hydraulic conductivity and depth to the fragipan layer, which together dictate how much the overlying soil can hold before runoff is observed.Therefore, estimating the b-parameter in this case requires an approach based on the effective moisture capacity that considers the hydraulic conductivity in addition to the depth to the restrictive layer.Such an approach has not been examined to date, despite the extensive coverage of fragipan soils in Southeastern U.S.This could potentially explain past difficulties/biases in simulating streamflows in the U.S. Southeast with land surface models.
Recent developments for calibrating land surface models like VIC have included the application of evolutionary algorithms for determining a set of parameters that best match model simulations with observed data, e.g., [19,25].Despite significant strides in optimization routines, the "search space" or range of values examined by these algorithms for each parameter is still user provided.An intrinsic problem in hydrologic systems modeling is that a given end state of the system can be reached by several potential means, corresponding to different sets of calibrated parameters [26].Thus, unless the ranges provided are physically based and restricted to capture for example the effective moisture capacity needed to represent fragipan soils, uncertainties in model predictions related to the calibration parameters will persist even if observed data are perfectly matched by model simulations [27][28][29].Consequently, the model's ability to forecast system response under climate and land use changes will be questionable since there is no guarantee that the model's effectiveness in predicting water budget trends is due to replication of the underlying physical processes, and not the model structure and parameter inputs.
Further uncertainty in the forecasting capability of land surface models may arise from the type of data that is used for model calibration.Currently, most water budget studies primarily utilize observed streamflow data at select locations along the stream network (typically the watershed outlet) for calibration purposes.Whereas, these data capture the integrated effects of rainfall-runoff-subsurface processes across the entire watershed, they do not necessarily capture the relative partitioning of rainfall into runoff, soil moisture, evapotranspiration (ET), groundwater recharge, and baseflow across spatially heterogeneous landscapes [19].Given the number of variables, it is plausible that rainfall can be partitioned differently into these components/fluxes and the same integrated streamflow fluxes would be observed at the calibration point, i.e., an equifinality issue [27].Hence, although the model may be able to replicate fluxes at the calibration point, inferences on the water budget within the watershed may be questionable.This is particularly true for regions with fragipan soils, where the layers are distributed over the watershed, both laterally and with depth, and are thus difficult to characterize in terms of their effects on both the surface and subsurface water fluxes.Thus, for forecasting studies, it is important also that other data, such as soil moisture and ET, are calibrated along with the streamflow data.The selection of additional calibration variables should not be arbitrary but rather guided by their ability to account for the major components of the water balance and closure of the budget.In so doing, model uncertainty is reduced and inferences may be made regarding potential changes in future water budgets with higher levels of confidence.This level of model calibration is now feasible given the availability of various spatiotemporal satellite data products for water budget components across the globe.
The main goal of this study is to examine the impacts of present and near-future climate on the water budget in the Obion River Watershed (ORW), Tennessee, a predominantly agricultural watershed with fragipan soils in the Southeastern U.S. The specific objectives are threefold: (1) Develop and provide a methodological blueprint for setting up, calibrating, and validating a VIC model in a fragipan region, (2) account for the effects of landscape heterogeneity on the partitioning of both surface and subsurface fluxes by considering ET and soil moisture as additional calibration parameters to streamflow, and (3) compare present and near-future trends in soil water fluxes in ORW and their implications for crop production.The calibration and validation blueprint presented herein calls for detailed soil information, including the extent of fragipan soils and infiltration data for each sub-drainage areas (SDA) within ORW.This information is used along with Horton's infiltration relationship [30] to provide physical ranges of the b-parameter.While the approach is used for a watershed with fragipan soils, it is not limited to these landscapes and offers a simple physical development process for the application of VIC to other types of watersheds as well.The study results provide insight on the water balance into the near-future to help farmers, watershed managers, and policymakers in decision making and resource management.

Study Area
The Obion River Watershed (HUC 08010202-3) is approximately 6400 km 2 and representative of the lower Mississippi River Basin, which is dominated by row crop agriculture [31] (Figure 2).The ORW has a humid, continental climate and loess-derived soils, making it ideal for agricultural production [32].As a result, intensive row crop agriculture now dominates the watershed once covered by forest and wooded wetlands.Intensive row crop agriculture and pastured grasslands cover 68.2% of the watershed while forests cover only 28.7%.Residential areas cover the remaining 3.1% (AVHRR land cover data).There are generally two major crop rotations in the cultivated areas; corn-soybean with winter wheat planted as a cover crop or as a double crop, and continuous cotton.Approximately 80 to 90% of the management is continuous no-till.watershed with fragipan soils in the Southeastern U.S. The specific objectives are threefold: (1) Develop and provide a methodological blueprint for setting up, calibrating, and validating a VIC model in a fragipan region, (2) account for the effects of landscape heterogeneity on the partitioning of both surface and subsurface fluxes by considering ET and soil moisture as additional calibration parameters to streamflow, and (3) compare present and near-future trends in soil water fluxes in ORW and their implications for crop production.The calibration and validation blueprint presented herein calls for detailed soil information, including the extent of fragipan soils and infiltration data for each sub-drainage areas (SDA) within ORW.This information is used along with Horton's infiltration relationship [30] to provide physical ranges of the b-parameter.While the approach is used for a watershed with fragipan soils, it is not limited to these landscapes and offers a simple physical development process for the application of VIC to other types of watersheds as well.The study results provide insight on the water balance into the near-future to help farmers, watershed managers, and policymakers in decision making and resource management.

Study Area
The Obion River Watershed (HUC 08010202-3) is approximately 6400 km 2 and representative of the lower Mississippi River Basin, which is dominated by row crop agriculture [31] (Figure 2).The ORW has a humid, continental climate and loess-derived soils, making it ideal for agricultural production [32].As a result, intensive row crop agriculture now dominates the watershed once covered by forest and wooded wetlands.Intensive row crop agriculture and pastured grasslands cover 68.2% of the watershed while forests cover only 28.7%.Residential areas cover the remaining 3.1% (AVHRR land cover data).There are generally two major crop rotations in the cultivated areas; corn-soybean with winter wheat planted as a cover crop or as a double crop, and continuous cotton.Approximately 80 to 90% of the management is continuous no-till.The annual average precipitation in the watershed is approximately 1310 mm.The majority of the rain occurs in the spring and winter, although heavy convective thunderstorms occur in the summer.The annual average temperature is 14.5 • C.During summer, ET can become the dominant flux in the hydrologic cycle, exceeding incoming precipitation.
The topography and soil properties of the watershed transition from its eastern boundary to the western boundary, which borders the Mississippi River (Figure 2).Rolling topography in the eastern headwaters contain several North-South bands of sand and clay formations.The local streams in this region have a moderately high gradient over generally sandy substrates.Loess plains in the middle of the watershed have a gently rolling topography and contain sand, clay, silt, and lignite capped by loess that can be 50-60 feet thick, generally deeper in the bluff regions towards the West.Streams in this part of the watershed are low-gradient and murky, with silty sand bottoms.Many of the streams have been deforested for agricultural purposes.Channel plugs have formed where aggradation and driftwood accumulate to form blockages and alter flow patterns.Along the Mississippi alluvial plain on the West, there are predominantly clayey soils with poor drainage that may contain wooded swampland and oxbow lakes.Most of the very large farms in ORW exist along the Mississippi River, with smaller farms located in other areas of the watershed.
Fragipan soils cover approximately 30% of ORW [33].These soils are mostly found in the uplands of the Eastern three-fourths of the watershed (Figure 2).The depth to the fragipan layer in ranges between 0.1 m and 2.0 m below the surface.In most locations, the fragipan layer is below historical tillage depths, typically 30 cm, and is therefore intact [34][35][36].Furthermore, no-till practices are currently practiced throughout the watershed to prevent additional disturbances to the fragipan layer.The thickness of most fragipan soils is at least 15 cm [24].The saturated hydraulic conductivity of the fragipan layer in ORW, based on the U.S. Natural Resources Conservation Service (NRCS) soil survey, is around 7 mm/h.The combination of poor drainage and relatively shallow depth to the fragipan layer creates a very flashy hydrologic system [37], especially during more intensive rainfall events, in areas where the fragipan soils are close to the surface.

Hydrologic Model
The three-layer variable infiltration capacity (VIC) model [14,38] was used for this study.VIC has been extensively calibrated and applied for a variety of applications on water resources management, land-atmosphere interactions, and climate change in many large and small basins over the U.S. and worldwide [10][11][12][13][39][40][41][42][43].Over the past two decades, VIC has been used in several multi-institutional projects such as Project for Intercomparison of Land surface Parameterization Schemes (PILPS) and the North American Land Data Assimilation System (NLDAS) project, and has been shown to compare well with other land surface models and available observed data [18,44,45].
The VIC model framework is described in References [14,46].VIC balances the water budget within individual grid cells (Figure 1), while considering multiple soil layers with variable infiltration and non-linear base flow.The soil characteristics are defined for each cell and held constant over time.One or more vegetation tiles describes the surface of each grid cell, and the vegetation characteristics such as Leaf Area Index, albedo, resistance, roughness root depth, and its relative fraction in each soil layer, are assigned for each tile.The top two soil layers represent the dynamic response of soil to infiltration.The model allows diffusion between the layers depending on the relative wetness.The model estimates ET using the Penman-Monteith equation.Moisture from the middle layer reaches the bottom layer through gravity drainage, which is regulated by a Brooks-Corey relationship.The bottom layer represents seasonal soil moisture responses except when the other soil layers are saturated.The base flow is only from the bottom layer and is based on the Arno model.
In this study, VIC was run at a daily time step in the water balance mode and results were aggregated at a monthly scale.A grid cell size of 1 km 2 was adopted to match the Advanced Very High Resolution Radiometer (AVHRR) land cover data used in the study, yielding a total of 6417 cells.
For each cell and time step, VIC computed water balance components such as runoff, evapotranspiration, soil moisture, and base flow (Figure 1).

Data Sources
ORW has a wealth of in-situ and remotely sensed data available for model development, calibration, and validation.The University of Tennessee Institute of Agriculture's Research and Education Center (REC) at Milan was a primary source of watershed and land-use parameters, as well as measurements of ET and soil moisture.Since 1962, the Milan station has been an active research facility in the heart of the watershed focusing on soil conservation, no-till agriculture, cover crops, and irrigation.
In addition to the REC data, six U.S. Geological Survey (USGS) stream gages within the watershed provided daily streamflow records for at least several years and up to several decades.Additional monthly ET data were obtained from the USGS Simplified Surface Energy Balance (SSEBop) data products [47], downloaded from the Geo Data Portal (https://cida.usgs.gov/gdp/).SSEBop combines reference ET values (from the Penman-Monteith equation) with remotely sensed thermal and vegetation imagery from the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite to generate actual evapotranspiration (AET) data.Soil moisture measurements were also obtained from nearby sites in the NRCS Soil Climate Analysis Network (SCAN).A 30-m Digital Elevation Model of the watershed was obtained from the Tennessee GIS data server (http://www.tngis.org/elevation.htm),and resampled to a 1 km resolution to match the model grid-cell size.
The observed historical climate data used in the study were obtained from DAYMET (http: //daymet.ornl.gov/)[48] for 1980-2011, at a spatial resolution of 1 km × 1 km.The weather parameters used included daily precipitation, minimum and maximum temperatures, humidity, and shortwave radiation.Modeled historical and future climate data were also used in the study.These were obtained from the fifth phase of the Coupled Model Intercomparison Project (CMIP5) [49], which provides long-term (hindcast/forecast) climate simulations for 1850-2205 with atmosphere-ocean global climate models (AOGCMs) coupled with Earth system models of intermediate complexity.The historical period of CMIP5 simulations ends in 2005 and the projection begins in 2006.For this study, data was obtained for 1980-2005 for the historical period and for 2006-2050 for the projection period.
Two Representative Concentration Pathways (RCPs) [50], RCP 4.5 and RCP 8.5, were considered for future climate projections.RCP 4.5 is a moderate climate scenario, where an additional radiative forcing of 4.5 W/m 2 compared to preindustrial conditions will be trapped in the atmosphere by 2100 RCP 8.5 is an extreme climate scenario, where an additional 8.5 W/m 2 will be trapped by 2100 relative to preindustrial conditions.For each RCP, climate predictions from the following four Global Climate Models (GCMs) were used: (1) The National Oceanic and Atmospheric Administration-Geophysical Fluid Dynamics Laboratory model (GFDL-ESM2M), ( 2

Selection of Calibration Parameters
A calibration process is presented herein for ORW based on its physical characteristics.The study builds on the Troy et al. work [19], which showed that calibration for hydrological modeling (at a spatial resolution of 1/8) using seven key parameters within an allowed range could significantly improve the model results for the conterminous U.S.These seven parameters are presented in Table 1 along with the ranges adopted in Reference [19].The method of determination of each of the parameters in this study is also specified in the table.For small-scale hydrological modeling, heterogeneity within sub-drainage areas should be considered in calibration for better representation of the watershed.The fine spatial resolution (1 km 2 grid cell) employed in this study provided the basis to account for the heterogeneity.Properties of each cell such as soil, vegetation, and climate were determined and assigned accordingly to represent the watershed heterogeneity.
According to Reference [14], the b-parameter has the largest effect on runoff, while the D s and W s parameters are critical in influencing the baseflow.Estimation of these three parameters is thus crucial for model performance.Unlike the other four parameters (i.e., D smax , Exp, and depths of soil layers 2 and 3), which are readily obtained from soil maps, b, D s , and W s are usually calibrated since they are generally unknown and can vary considerably within a watershed [52].In this study, b was estimated using soil data and Horton's infiltration relationship [30], while D s and W s were estimated by evaluating plausible physical values for both parameters.The parameters D smax and Exp were also calculated using empirical relationships (see Table 1).The depths of the second and third soil layers were chosen to be consistent with the soil layers of the CONUS-Soil data [53].For the third layer, however, a sensitivity analysis was further performed to determine a minimum required depth of 1.4 m.Cells with layer 3 depths lower than this value had fast drainage and unrealistically low soil moisture contents.

Estimation of Calibration Parameters
As discussed above, the b-parameter is an important model parameter and should be determined carefully to make sure that model results are physically meaningful.Based on its definition, b controls the partitioning of rainfall into runoff and infiltration, and ranges commonly between 0 and 0.4, though larger values have been estimated for some watersheds [19,54].Higher values produce more runoff and less infiltration whereas lower values produce less runoff and more infiltration.An extensive literature search for the parameterization of the b-parameter (see Table A1 in the Appendix A) reveals two readily noticeable facts: (1) The range of values used in the literature is wide; and (2) there is a lack of a consistent methodology used to derive this parameter.Sivapalan and Woods [10] and Huang et al. [20] are two of the few studies that have attempted to determine b directly from soil observations.However, as explained in the introduction, a shortcoming of these studies is the assumption of a static soil moisture capacity distribution, which is not directly applicable to watersheds with distributions of fragipan soils that vary in space, both laterally and with depth.To address this issue, this study adopts a dynamic or effective soil moisture capacity distribution assumption for determining the value of b.
Since fragipan layers regulate infiltration, leading to perched saturation zones and increased runoff [55], the effective soil moisture capacity for fragipan regions is estimated herein with the following relationship derived from integration of Horton's curve: where i e is the effective soil moisture capacity (L), f c is the final constant infiltration capacity, i.e., saturated hydraulic conductivity (LT −1 ), f 0 and k are, respectively, the initial infiltration capacity (LT −1 ) and decay rate (T −1 ) and t s is the time to steady infiltration (T).Equation ( 1) accounts for dynamic or temporal effects that control water storage and saturation in the soil.Field data for fragipan sites within the watershed were used to determine the variables on the right hand side of Equation 1, from which i e was then computed.The computed values of i e across the watershed were used to establish the cumulative soil moisture capacity curve, from which b was determined by fitting the following distribution [54]: where i m is the maximum soil moisture capacity (L), and A is the fraction of the area for which the infiltration capacity is less than i e .Further details on the estimation of i e and b for ORW is provided below.
For the fragipan areas in the watershed, values of f 0 and k averaged 210 mm/h and 2 min −1 , respectively [56].The average f c value for fragipan soils, based on the NRCS soil survey, was also approximately 7 mm/h.The t s determined from field experiments was on average 60 h for an approximate soil depth of 2 m, with fragipan soils occurring around 600 mm below the surface.Using Equation ( 1), an average i e value of 432 mm could then be estimated for the fragipan regions (i.e., an effective soil moisture capacity approximately 51% of the total capacity of the soil column, which has a porosity of 0.42).Repeating the step for different fragipan depths, the cumulative distribution of effective soil moisture capacities were plotted and Equation ( 2) was fitted to the data to obtain an estimate of 0.53 for the b value.This is shown in Figure 3.Following a similar approach, a b value of 0.04 was determined for the non-fragipan regions.Since fragipan layers regulate infiltration, leading to perched saturation zones and increased runoff [55], the effective soil moisture capacity for fragipan regions is estimated herein with the following relationship derived from integration of Horton's curve: where ie is the effective soil moisture capacity (L), fc is the final constant infiltration capacity, i.e., saturated hydraulic conductivity (LT −1 ), f0 and k are, respectively, the initial infiltration capacity (LT −1 ) and decay rate (T −1 ) and ts is the time to steady infiltration (T).Equation ( 1) accounts for dynamic or temporal effects that control water storage and saturation in the soil.Field data for fragipan sites within the watershed were used to determine the variables on the right hand side of Equation 1, from which ie was then computed.The computed values of ie across the watershed were used to establish the cumulative soil moisture capacity curve, from which b was determined by fitting the following distribution [54]: where im is the maximum soil moisture capacity (L), and A is the fraction of the area for which the infiltration capacity is less than ie.Further details on the estimation of ie and b for ORW is provided below.
For the fragipan areas in the watershed, values of f0 and k averaged 210 mm/h and 2 min −1 , respectively [56].The average fc value for fragipan soils, based on the NRCS soil survey, was also approximately 7 mm/h.The ts determined from field experiments was on average 60 h for an approximate soil depth of 2 m, with fragipan soils occurring around 600 mm below the surface.Using Equation ( 1), an average ie value of 432 mm could then be estimated for the fragipan regions (i.e., an effective soil moisture capacity approximately 51% of the total capacity of the soil column, which has a porosity of 0.42).Repeating the step for different fragipan depths, the cumulative distribution of effective soil moisture capacities were plotted and Equation ( 2) was fitted to the data to obtain an estimate of 0.53 for the b value.This is shown in Figure 3.Following a similar approach, a b value of 0.04 was determined for the non-fragipan regions.Since VIC requires a single value of b per grid cell, it was necessary to determine an effective b value for each cell based on the relative influences of fragipan and non-fragipan areas within the cell.This required a sensitivity analysis of the b-parameter to determine first whether or not variations in b had a significant impact on rainfall-runoff-baseflow flux proportioning, and then to determine its relationship with the runoff coefficient (i.e., ratio of runoff to rainfall) to establish weighting factors for computing the effective b value from the fragipan and non-fragipan b values and area proportions.As mentioned above, fragipan soils cover approximately 30% of ORW.However, the actual Since VIC requires a single value of b per grid cell, it was necessary to determine an effective b value for each cell based on the relative influences of fragipan and non-fragipan areas within the cell.This required a sensitivity analysis of the b-parameter to determine first whether or not variations in b had a significant impact on rainfall-runoff-baseflow flux proportioning, and then to determine its relationship with the runoff coefficient (i.e., ratio of runoff to rainfall) to establish weighting factors for computing the effective b value from the fragipan and non-fragipan b values and area proportions.As mentioned above, fragipan soils cover approximately 30% of ORW.However, the actual proportion of fragipan to non-fragipan area within each grid cell will vary from cell to cell depending on the landscape characteristics [57].For the sake of simplicity, an assumption of 30% fragipan coverage was made for each grid cell falling within a sub-drainage areas where fragipan soils were present.Since the VIC flux predictions are aggregated herein to the sub-drainage scale, this assumption is valid because fragipan effects at that scale are aggregated in reality and will reflect the approximate coverage in each sub-drainage area.
The sensitivity analysis was performed for a two-year period (2010-2011).The b-parameter was varied between 0.15 and 0.9 and its effects on the runoff, baseflow, and streamflow were analyzed.The results are presented in Figure 4 along with the simulated precipitation.As seen, higher b values produced more runoff (Figure 4a) and less baseflow (Figure 4b).For streamflow, this resulted in higher hydrograph peaks and lower troughs (Figure 4c).Overall, runoff differed by as much as 150% (a factor of 2.6) in May of the first year, while baseflow differed by as much as 170% (a factor of 2.7), indicating a sensitivity of the water fluxes to the value of the b-parameter.In terms of streamflow, the differences in peaks were only as much as 50% (a factor of 1.5) in May of the first year while differences in troughs were only as much as 100% (a factor of 2) in August of the first year.Given the large difference between the estimated b-parameter for fragipan (~0.53) and non-fragipan regions (~0.04) in ORW, the sensitivity analysis suggests that the separate estimation for fragipan and non-fragipan areas is justified since significantly different responses in terms of runoff, baseflow, and streamflow are expected.
Geosciences 2018, 8, x FOR PEER REVIEW 9 of 25 proportion of fragipan to non-fragipan area within each grid cell will vary from cell to cell depending on the landscape characteristics [57].For the sake of simplicity, an assumption of 30% fragipan coverage was made for each grid cell falling within a sub-drainage areas where fragipan soils were present.Since the VIC flux predictions are aggregated herein to the sub-drainage scale, this assumption is valid because fragipan effects at that scale are aggregated in reality and will reflect the approximate coverage in each sub-drainage area.
The sensitivity analysis was performed for a two-year period (2010-2011).The b-parameter was varied between 0.15 and 0.9 and its effects on the runoff, baseflow, and streamflow were analyzed.The results are presented in Figure 4 along with the simulated precipitation.As seen, higher b values produced more runoff (Figure 4a) and less baseflow (Figure 4b).For streamflow, this resulted in higher hydrograph peaks and lower troughs (Figure 4c).Overall, runoff differed by as much as 150% (a factor of 2.6) in May of the first year, while baseflow differed by as much as 170% (a factor of 2.7), indicating a sensitivity of the water fluxes to the value of the b-parameter.In terms of streamflow, the differences in peaks were only as much as 50% (a factor of 1.5) in May of the first year while differences in troughs were only as much as 100% (a factor of 2) in August of the first year.Given the large difference between the estimated b-parameter for fragipan (~0.53) and non-fragipan regions (~0.04) in ORW, the sensitivity analysis suggests that the separate estimation for fragipan and nonfragipan areas is justified since significantly different responses in terms of runoff, baseflow, and streamflow are expected.The relationship between the b-parameter and the normalized runoff coefficient is plotted in Figure 5 based on the sensitivity analysis results (each dataset is normalized with the runoff coefficient corresponding to b = 0.9).The overall trend for the mean runoff coefficient (solid red line in Figure 5) follows an approximate power law (Hoerl curve).Using this relationship, an effective normalized runoff coefficient of 0.4 was estimated for a cell with 30% fragipan (b = 0.53) and 70% nonfragipan (b = 0.04) coverage.From this effective normalized runoff coefficient value, a corresponding effective b parameter of 0.15 was estimated, which was adopted in the model.The relationship between the b-parameter and the normalized runoff coefficient is plotted in Figure 5 based on the sensitivity analysis results (each dataset is normalized with the runoff coefficient corresponding to b = 0.9).The overall trend for the mean runoff coefficient (solid red line in Figure 5) follows an approximate power law (Hoerl curve).Using this relationship, an effective normalized runoff coefficient of 0.4 was estimated for a cell with 30% fragipan (b = 0.53) and 70% non-fragipan (b = 0.04) coverage.From this effective normalized runoff coefficient value, a corresponding effective b parameter of 0.15 was estimated, which was adopted in the model.The common ranges of Ds and Ws adopted from the literature are, respectively, 0.001-1 and 0.2-1 [9] (see Table 2).The selected values for calibrating Ds, and Ws used in this study (presented in Table 1) were chosen after reviewing the literature for VIC simulations of regions close to west Tennessee or with relatively similar soil and/or hydrological regime to ORW (see Table 2).
Based on the soil properties, model parameter values for Dsmax ranged from 0 to 100 for the 6147 cells, with an average value of around 10. Likewise, parameter values for Exp ranged from 3 to 20, also with an average value of around 10. Finally, the depths of the 2nd and 3rd layers generally fell between 0.4 m and 1.5 m.The common ranges of D s and W s adopted from the literature are, respectively, 0.001-1 and 0.2-1 [9] (see Table 2).The selected values for calibrating D s , and W s used in this study (presented in Table 1) were chosen after reviewing the literature for VIC simulations of regions close to west Tennessee or with relatively similar soil and/or hydrological regime to ORW (see Table 2).
Based on the soil properties, model parameter values for D smax ranged from 0 to 100 for the 6147 cells, with an average value of around 10. Likewise, parameter values for Exp ranged from 3 to 20, also with an average value of around 10. Finally, the depths of the 2nd and 3rd layers generally fell between 0.4 m and 1.5 m.

Model Calibration
The model was calibrated by comparing simulated monthly streamflow fluxes for each of the six SDAs (Figure 2) against corresponding observed monthly USGS streamflow data for the period of 1980-2011.Relative bias and the Nash-Sutcliffe efficiency coefficient [60] were used as the measures of performance for each calibration case.The best set of calibration parameters and model performance values for each SDA, along with the USGS gauging stations used, are presented in Table 3.A comparison between simulated streamflow hydrographs and corresponding USGS streamflow hydrographs for two of the SDAs is also shown in Figure 6.As seen, the calibrated model correctly predicts the overall trends in streamflow and captures well the low flows.However, it tends to underestimate high peak flows.The Nash-Sutcliffe Efficiency coefficient above 0.65 for all SDAs suggests that the model results are good [61][62][63].0.9 0.9 0.73 0.12 07026040 6 0.9 0.9 0.73 0.16 07026300

Model Calibration
The model was calibrated by comparing simulated monthly streamflow fluxes for each of the six SDAs (Figure 2) against corresponding observed monthly USGS streamflow data for the period of 1980-2011.Relative bias and the Nash-Sutcliffe efficiency coefficient [60] were used as the measures of performance for each calibration case.The best set of calibration parameters and model performance values for each SDA, along with the USGS gauging stations used, are presented in Table 3.A comparison between simulated streamflow hydrographs and corresponding USGS streamflow hydrographs for two of the SDAs is also shown in Figure 6.As seen, the calibrated model correctly predicts the overall trends in streamflow and captures well the low flows.However, it tends to underestimate high peak flows.The Nash-Sutcliffe Efficiency coefficient above 0.65 for all SDAs suggests that the model results are good [61][62][63].
The effective soil moisture capacity approach used herein for determining b gives confidence that the model is able to partition correctly precipitation to infiltration and surface runoff.To validate the model's ability to simulate and partition subsurface fluxes, the calibrated results were also compared with observed data of AET and soil moisture.Figure 7 compares the MODIS-derived AET with the VIC-simulated AET for the entire watershed from 2000 to 2012.As seen, the MODIS values are generally within one standard deviation of the VIC values and the Nash-Sutcliffe efficiency coefficient is 0.77, which suggests good agreement.Furthermore, the simulated AET values are also similar to the average range of 25-125 mm per month reported for the Arkansas-Red River basin (which neighbors ORW on the west) in Reference [58].The effective soil moisture capacity approach used herein for determining b gives confidence that the model is able to partition correctly precipitation to infiltration and surface runoff.To validate the model's ability to simulate and partition subsurface fluxes, the calibrated results were also compared with observed data of AET and soil moisture.

Climate Impact on the Obion River Watershed Water Budget
Figure 9 presents annual precipitation and temperature time series in ORW for the historical and projection periods examined herein.In each plot, the dotted blue lines represent the maximum and minimum model predictions (among all four GCMs), and the solid black curve shows the observed historical data.The solid cyan line is the ensemble average of the four GCMs.
While there are no obvious trends in precipitation for both the RCP 4.5 and RCP 8.5 scenarios, RCP 4.5 shows a slight decline, whilst RCP 8.5 shows a slight increase.This is accordance with reported trends for the region for the 20th century.According to Hanson and Weltzin [64], the region had about a 1% increase of precipitation per decade over a 100-year period.It was also reported that the number of moderate and extreme dry months (according to the Palmer Index) in the second half of the same 100-year period was less than the first half [64].
On the contrary, historic and projected temperatures clearly exhibit a rising trend with an approximate increase of 3 °C over the 70-year period from 1980 to 2050, for both RCP 4.5 and RCP 8.5 scenarios.Consistent with these are the annual time series of potential evapotranspiration (PET), which also show an increasing trend for both RCP 4.5 and 8.5 till 2050 (~1.4 mm/year), with no

Climate Impact on the Obion River Watershed Water Budget
Figure 9 presents annual precipitation and temperature time series in ORW for the historical and projection periods examined herein.In each plot, the dotted blue lines represent the maximum and minimum model predictions (among all four GCMs), and the solid black curve shows the observed historical data.The solid cyan line is the ensemble average of the four GCMs.
While there are no obvious trends in precipitation for both the RCP 4.5 and RCP 8.5 scenarios, RCP 4.5 shows a slight decline, whilst RCP 8.5 shows a slight increase.This is accordance with reported trends for the region for the 20th century.According to Hanson and Weltzin [64], the region had about a 1% increase of precipitation per decade over a 100-year period.It was also reported that the number of moderate and extreme dry months (according to the Palmer Index) in the second half of the same 100-year period was less than the first half [64].
On the contrary, historic and projected temperatures clearly exhibit a rising trend with an approximate increase of 3 °C over the 70-year period from 1980 to 2050, for both RCP 4.5 and RCP 8.5 scenarios.Consistent with these are the annual time series of potential evapotranspiration (PET), which also show an increasing trend for both RCP 4.5 and 8.5 till 2050 (~1.4 mm/year), with no

Climate Impact on the Obion River Watershed Water Budget
Figure 9 presents annual precipitation and temperature time series in ORW for the historical and projection periods examined herein.In each plot, the dotted blue lines represent the maximum and minimum model predictions (among all four GCMs), and the solid black curve shows the observed historical data.The solid cyan line is the ensemble average of the four GCMs.
While there are no obvious trends in precipitation for both the RCP 4.5 and RCP 8.5 scenarios, RCP 4.5 shows a slight decline, whilst RCP 8.5 shows a slight increase.This is accordance with reported trends for the region for the 20th century.According to Hanson and Weltzin [64], the region had about a 1% increase of precipitation per decade over a 100-year period.It was also reported that the number of moderate and extreme dry months (according to the Palmer Index) in the second half of the same 100-year period was less than the first half [64].
On the contrary, historic and projected temperatures clearly exhibit a rising trend with an approximate increase of 3 • C over the 70-year period from 1980 to 2050, for both RCP 4.5 and RCP 8.5 scenarios.Consistent with these are the annual time series of potential evapotranspiration (PET), which also show an increasing trend for both RCP 4.5 and 8.5 till 2050 (~1.4 mm/year), with no significant differences between the scenarios (see Figure 10).This is to be expected since PET is a function of the mean temperature, among others.
function of the mean temperature, among others.
Unlike PET, the annual time series of AET show different trends (see Figure 11).While the projected AET for RCP 8.5 is increasing, the projected AET for RCP 4.5 is declining.This phenomenon can be explained by the precipitation trends shown in Figure 8.Consistent with AET trends, the projected precipitation declines for RCP 4.5 and rises for RCP 8.5.Studies indicate that changes in wind and air temperature may not always be the root cause for AET changes, and that changes in precipitation characteristics can dominate AET trend changes [65].The decline in precipitation for RCP 4.5 results in a soil moisture limitation that primarily drives declines in the ET [66].Thus, the available water restricts variation of AET for this region in West Tennessee.significant differences between the scenarios (see Figure 10).This is to be expected since PET is a function of the mean temperature, among others.
Unlike PET, the annual time series of AET show different trends (see Figure 11).While the projected AET for RCP 8.5 is increasing, the projected AET for RCP 4.5 is declining.This phenomenon can be explained by the precipitation trends shown in Figure 8.Consistent with AET trends, the projected precipitation declines for RCP 4.5 and rises for RCP 8.5.Studies indicate that changes in wind and air temperature may not always be the root cause for AET changes, and that changes in precipitation characteristics can dominate AET trend changes [65].The decline in precipitation for RCP 4.5 results in a soil moisture limitation that primarily drives declines in the ET [66].Thus, the available water restricts variation of AET for this region in West Tennessee.Unlike PET, the annual time series of AET show different trends (see Figure 11).While the projected AET for RCP 8.5 is increasing, the projected AET for RCP 4.5 is declining.This phenomenon can be explained by the precipitation trends shown in Figure 8.Consistent with AET trends, the projected precipitation declines for RCP 4.5 and rises for RCP 8.5.Studies indicate that changes in wind and air temperature may not always be the root cause for AET changes, and that changes in precipitation characteristics can dominate AET trend changes [65].The decline in precipitation for RCP 4.5 results in a soil moisture limitation that primarily drives declines in the ET [66].Thus, the available water restricts variation of AET for this region in West Tennessee.Figure 12 shows the soil moisture time series for the historical and projected periods.The projected time series for both RCPs are within the same range as the historical period.In addition, the slope of the trend for both RCPs is negligible, suggesting no significant changes due to climate in the near future till 2050.It is worth noting that variation of soil moisture through ET also depends on land use/land cover [67].In this study, the land use/land cover was fixed to current values for the projection period.Any land use/land cover changes in future may directly affect ET, and consequently the soil moisture.Figure 12 shows the soil moisture time series for the historical and projected periods.The projected time series for both RCPs are within the same range as the historical period.In addition, the slope of the trend for both RCPs is negligible, suggesting no significant changes due to climate in the near future till 2050.It is worth noting that variation of soil moisture through ET also depends on land use/land cover [67].In this study, the land use/land cover was fixed to current values for the projection period.Any land use/land cover changes in future may directly affect ET, and consequently the soil moisture.
The above trends for precipitation, temperature, PET, and AET, along with trends in runoff and baseflow, are summarized in Figure 13.Box plots of the parameters for the historical period are compared to the associated box plots projected for RCPs 4.5 and 8.5.The comparisons are for the 20-year periods of 1986-2005 in the past and 2031-2050 in the future.For comparison purposes, all parameters are presented as deviations from the average value during the historical period.Precipitation and temperature are known drivers of change in the water balance.For ORW, while there are slight changes in precipitation, there are significant changes in the temperature for the projected period.However, as noted previously, the water balance in this region is more restricted by water availability.Therefore, the projected changes in the other water budget parameters, including AET, soil moisture, baseflow, and runoff are slight, except some changes in the most extreme data points manifested by extension of the whiskers in Figure 13.
projected time series for both RCPs are within the same range as the historical period.In addition, the slope of the trend for both RCPs is negligible, suggesting no significant changes due to climate in the near future till 2050.It is worth noting that variation of soil moisture through ET also depends on land use/land cover [67].In this study, the land use/land cover was fixed to current values for the projection period.Any land use/land cover changes in future may directly affect ET, and consequently the soil moisture.The above trends for precipitation, temperature, PET, and AET, along with trends in runoff and baseflow, are summarized in Figure 13.Box plots of the parameters for the historical period are compared to the associated box plots projected for RCPs 4.5 and 8.5.The comparisons are for the 20year periods of 1986-2005 in the past and 2031-2050 in the future.For comparison purposes, all parameters are presented as deviations from the average value during the historical period.Precipitation and temperature are known drivers of change in the water balance.For ORW, while there are slight changes in precipitation, there are significant changes in the temperature for the projected period.However, as noted previously, the water balance in this region is more restricted by water availability.Therefore, the projected changes in the other water budget parameters, including AET, soil moisture, baseflow, and runoff are slight, except some changes in the most extreme data points manifested by extension of the whiskers in Figure 13.Statistical tests were used to assess quantitatively changes in the mean and standard deviation of each parameter from the historical period to the projected period.The same 20-year periods considered in Figure 13 were used for these tests.T-tests were used to measure differences in the mean while F-tests were used to compare standard deviations.The results of the tests are presented in Table 4.Only two parameters, i.e., temperature and PET, show significant changes in their mean values into the projected period for both RCP 4.5 and RCP 8.5.Further, the only parameter that shows a statistically significant change in its standard deviation is the temperature for RCP 8.5.The findings suggest that there will be a significant increase in extreme temperature conditions under the RCP 8.5 scenario.Overall, however, the statistical results suggest that changes in projected water budget will be not significant in near future for the ORW.
The characteristics of ORW were also examined with the Budyko curve [68], which is a common representation of the land-water balance that describes the mean-annual partitioning of precipitation into streamflow and evaporation.The curve summarizes the relation between the ratios of AET/P (evaporative index) and PET/P (dryness index), where P is precipitation.The distribution of annual values on the Budyko curve (Figure 14) can be used to determine the following characteristics of a Statistical tests were used to assess quantitatively changes in the mean and standard deviation of each parameter from the historical period to the projected period.The same 20-year periods considered in Figure 13 were used for these tests.T-tests were used to measure differences in the mean while F-tests were used to compare standard deviations.The results of the tests are presented in Table 4.Only two parameters, i.e., temperature and PET, show significant changes in their mean values into the projected period for both RCP 4.5 and RCP 8.5.Further, the only parameter that shows a statistically significant change in its standard deviation is the temperature for RCP 8.5.The findings suggest that there will be a significant increase in extreme temperature conditions under the RCP 8.5 scenario.Overall, however, the statistical results suggest that changes in projected water budget will be not significant in near future for the ORW.The characteristics of ORW were also examined with the Budyko curve [68], which is a common representation of the land-water balance that describes the mean-annual partitioning of precipitation into streamflow and evaporation.The curve summarizes the relation between the ratios of AET/P (evaporative index) and PET/P (dryness index), where P is precipitation.The distribution of annual values on the Budyko curve (Figure 14) can be used to determine the following characteristics of a watershed [69]: 1.
Resistance (or responsivity), which measures the degree to which streamflow is synchronized with precipitation, and; 2.
Resilience (or elasticity), which measures the degree to which a watershed can return to normal functioning following perturbations.
Geosciences 2018, 8, x FOR PEER REVIEW 17 of 25 1. Resistance (or responsivity), which measures the degree to which streamflow is synchronized with precipitation, and; 2. Resilience (or elasticity), which measures the degree to which a watershed can return to normal functioning following perturbations.
Resistance was estimated from the deviation in the AET/P of annual values (i.e., Δy-axis), while resilience was estimated as the ratio of deviation in PET/P to the deviation in AET/P (i.e., Δx-axis/Δyaxis).The data suggest that ORW has a relatively low resistance (i.e., high responsivity = 0.9), meaning that streamflow is largely synchronous with precipitation.This result is consistent with the effects of the fragipan soils in the region, which systematically transfers precipitation into streamflow discharge with less storage.Further, the ORW has a high resilience (i.e., elasticity = 2.3 (> 1)), so it has the capacity to sustain its functionality after extreme climatic perturbation such as drought or extreme precipitation [69].Resistance was estimated from the deviation in the AET/P of annual values (i.e., ∆y-axis), while resilience was estimated as the ratio of deviation in PET/P to the deviation in AET/P (i.e., ∆x-axis/∆y-axis).The data suggest that ORW has a relatively low resistance (i.e., high responsivity = 0.9), meaning that streamflow is largely synchronous with precipitation.This result is consistent with the effects of the fragipan soils in the region, which systematically transfers precipitation into streamflow discharge with less storage.Further, the ORW has a high resilience (i.e., elasticity = 2.3 (> 1)), so it has the capacity to sustain its functionality after extreme climatic perturbation such as drought or extreme precipitation [69].

Discussion and Conclusions
Biases in streamflow predictions in the U.S. Southeast from past hydrological modeling studies have been attributed to calibration issues [18,19].Since other regions with similar climate have not shown a similar bias, landscape characteristics are believed to be the cause of the bias.Specifically, the occurrence and extensive coverage of fragipan soils in the region are thought to control hydrologic response by affecting the partitioning of rainfall into runoff and infiltrated water, especially during more intensive rainfall events, in areas where the fragipan soils are close to the surface.
This study presented a methodology for calibrating regions with significant fragipan coverage based on the effective soil moisture capacity associated with these features.Estimation of the b-parameter of the VIC model, which controls rainfall-runoff-infiltration partitioning, was based on fragipan layer characteristics such as the saturated hydraulic capacity, the time to saturation, and the depth to the fragipan layer.This approach avoided the need to evaluate the wide range of b-parameter values recommended in literature (see Table A1) to match the observed flows.Further, since the approach is based on the physical characteristics of the soil, simulated trends in runoff and baseflow, as well as the behavior of the watershed in terms water budget response can be attributed to the processes that take place.
The proposed approach was successfully applied to the Obion River Watershed located in the U.S. Southeast, specifically the Northwestern part of Tennessee.Estimated b-parameter values of 0.53 and 0.04 for fragipan and non-fragipan soils, respectively, in the watershed produced good model predictions with Nash-Sutcliffe coefficient values greater than 0.65 and relative bias values ranging between 0 and 0.32 for sub-drainage areas within the watershed.A sensitivity analysis on the b-parameter validated the need to examine fragipan and non-fragipan soils separately in the watershed and confirmed that fragipan soils were resulting in flashier streamflow response with higher peaks and lower troughs.Further, the calibrated model was also shown to correctly predict trends in evapotranspiration and soil moisture, confirming that it was able to correctly partition subsurface fluxes in the various components.
The calibrated model was used to study the impacts of climate on the water balance in ORW.Water budgets were simulated for past and future projected climates up to 2050 (based on four GCMs for RCPs 4.5 and 8.5).While no significant changes in the water balance were forecasted for the near future (based on T-tests and F-tests), the results suggested a significant change in extreme temperature conditions for the RCP 8.5 scenario.Trends in the fluxes also suggested that the water budget in ORW is controlled to a large degree by water availability.This explains the lack of sensitivity of the water budget to forecasted temperature increases.The effect of a temperature increase will be to increase PET.However, since water availability is limiting, the AET will be regulated by the water input (i.e., precipitation) and not temperature [66], thereby rendering the watershed less sensitive to temperature changes.Being a cropland, this scenario occurs primarily during the growing season in ORW when most of the ET occurs.Nonetheless, the effects are reflected in the annual water budget trends since ET is the primary system water output on an annual basis (see Table 4).
Annual changes in the evaporative (AET/P) and dryness (PET/P) indices of ORW on the Budyko curve were used to characterize its responsivity (i.e., the degree to which streamflow is synchronized with precipitation) and resilience (i.e., the degree to which a watershed can return to normal functioning following perturbations).The ORW was found to have a high responsivity of 0.9, suggesting that stream flow was largely synchronized with precipitation.This finding is consistent with the effects of the fragipan soils in the region, which tends to systematically transfer precipitation into streamflow discharge with relatively less storage.The ORW was also found to have a high resilience of 2.3, also indicating that it has the capacity to sustain its functionality after extreme climatic perturbation such as drought or extreme precipitation.
The projection of no significant changes in the water budget of ORW in the near future has implications to the region's agriculture.While the irrigation for agricultural purposes in Tennessee has been increasing over past decades and is expected to increase over the next decades [70], it is well below the average in the country.The average acres of farms under irrigation was 17.8% for Tennessee and 25.4% for the United States with average acre-feet irrigation per acre of 0.4 for Tennessee and 1.6 for the United States [71].Since the mean soil moisture and evapotranspiration are not expected to change profoundly until 2050, agriculture in Western Tennessee can be expected to continue to perform similarly to recent trends, holding economic factors constant.Furthermore, water sources for irrigation are not expected to be impacted in the near future from climate changes according to the projected water balance.

Figure 1 .
Figure 1.Depiction of soil water balance based on the continuous streamflow model [7,8].

Figure 1 .
Figure 1.Depiction of soil water balance based on the continuous streamflow model [7,8].
) the Met Office Hadley Centre model (HadGEM2-CC), (3) the Institut Pierre-Simon Laplace model (IPSL-CM5A-MR), and (4) the Japan Agency for Marine-Earth Science and Technology, Atmosphere and Ocean Research Institute at the University of Tokyo, and National Institute for Environmental Studies model (MIROC-ESM-CHEM).

Figure 3 .
Figure 3. Calibrated variable infiltration curve for fragipan sites in ORW showing data points computed with the effective soil moisture capacity approach.

Figure 3 .
Figure 3. Calibrated variable infiltration curve for fragipan sites in ORW showing data points computed with the effective soil moisture capacity approach.

Figure 4 .
Figure 4. Sensitivity analysis of b and its effects on runoff (a), baseflow (b), and streamflow (c).The plots shows the monthly fluxes for the two years of 2010-2011.

Figure 4 .
Figure 4. Sensitivity analysis of b and its effects on runoff (a), baseflow (b), and streamflow (c).The plots shows the monthly fluxes for the two years of 2010-2011.

Geosciences 2018, 8 , 25 Figure 5 .
Figure 5. Variation of normalized runoff coefficient with the b parameter.The black circles are data points from the sensitivity analysis and the solid red line is a Hoerl model trendline fitted to the data.

Figure 5 .
Figure 5. Variation of normalized runoff coefficient with the b parameter.The black circles are data points from the sensitivity analysis and the solid red line is a Hoerl model trendline fitted to the data.

Figure 7
compares the MODIS-derived AET with the VIC-simulated AET for the entire watershed from 2000 to 2012.As seen, the MODIS values are generally within one standard deviation of the VIC values and the Nash-Sutcliffe efficiency coefficient is 0.77, which suggests good agreement.Furthermore, the simulated AET values are also similar to the average range of 25-125 mm per month reported for the Arkansas-Red River basin (which neighbors ORW on the west) in Reference [58].

Figure 7 .
Figure 7. Validation of monthly AET for the whole watershed.The model average and variability represents, respectively, spatially averaged monthly flux and the spatial standard deviation.

Figure 8
Figure 8 shows the comparison between the simulated and observed soil moisture data from 2004 to 2010.The observed data in this case are field observations obtained from the NRCS-SCAN network.The comparison confirms that the model was able to replicate the magnitude and trends in soil moisture in the watershed, with a fair Nash-Sutcliffe Efficiency Coefficient of 0.4.

Figure 8 .
Figure 8. Validation of monthly soil moisture for the whole watershed.The model average and variability represents, respectively, spatially averaged monthly flux, and the spatial standard deviation.

Figure 7 .
Figure 7. Validation of monthly AET for the whole watershed.The model average and variability represents, respectively, spatially averaged monthly flux and the spatial standard deviation.

Figure 8 25 Figure 7 .
Figure 8 shows the comparison between the simulated and observed soil moisture data from 2004 to 2010.The observed data in this case are field observations obtained from the NRCS-SCAN network.The comparison confirms that the model was able to replicate the magnitude and trends in soil moisture in the watershed, with a fair Nash-Sutcliffe Efficiency Coefficient of 0.4.

Figure 8
Figure 8 shows the comparison between the simulated and observed soil moisture data from 2004 to 2010.The observed data in this case are field observations obtained from the NRCS-SCAN network.The comparison confirms that the model was able to replicate the magnitude and trends in soil moisture in the watershed, with a fair Nash-Sutcliffe Efficiency Coefficient of 0.4.

Figure 8 .
Figure 8. Validation of monthly soil moisture for the whole watershed.The model average and variability represents, respectively, spatially averaged monthly flux, and the spatial standard deviation.

Figure 8 .
Figure 8. Validation of monthly soil moisture for the whole watershed.The model average and variability represents, respectively, spatially averaged monthly flux, and the spatial standard deviation.

Figure 9 .
Figure 9. Annual precipitation (top) and temperature (bottom) time series based on observations for 1980-2011 and modeled climate from the fifth phase of the Coupled Model Intercomparison Project (CMIP5) [49] for historical (1980-2005) and projected (2006-2050) period according to two Representative Concentration Pathways (RCPs) [50] of RCP 4.5 (left) and RCP 8.5 (right).The average curve is the average of four GCM models, and the variability curves are the minimum and maximum of four models.The red line is the slope of the time series during the projected period.The percent change from 2006 to 2050 is shown in red above each line.

Figure 9 .
Figure 9. Annual precipitation (top) and temperature (bottom) time series based on observations for 1980-2011 and modeled climate from the fifth phase of the Coupled Model Intercomparison Project (CMIP5) [49] for historical (1980-2005) and projected (2006-2050) period according to two Representative Concentration Pathways (RCPs) [50] of RCP 4.5 (left) and RCP 8.5 (right).The average curve is the average of four GCM models, and the variability curves are the minimum and maximum of four models.The red line is the slope of the time series during the projected period.The percent change from 2006 to 2050 is shown in red above each line.

Figure 9 .
Figure 9. Annual precipitation (top) and temperature (bottom) time series based on observations for 1980-2011 and modeled climate from the fifth phase of the Coupled Model Intercomparison Project (CMIP5) [49] for historical (1980-2005) and projected (2006-2050) period according to two Representative Concentration Pathways (RCPs) [50] of RCP 4.5 (left) and RCP 8.5 (right).The average curve is the average of four GCM models, and the variability curves are the minimum and maximum of four models.The red line is the slope of the time series during the projected period.The percent change from 2006 to 2050 is shown in red above each line.

Figure 10 .
Figure 10.Time series of annual PET based on four GCM models for the historical period (1980-2005) and projected period (2006-2050) considering RCP 4.5 (left) and RCP 8.5 (right).The average curve is the average of four GCM models, and the variability curves are the minimum and maximum of four models.The red line is the slope of the time series during the projected period.The percent change from 2006 to 2050 is shown in red above each line.

Geosciences 2018, 8 , 25 Figure 10 .
Figure 10.Time series of annual PET based on four GCM models for the historical period (1980-2005) and projected period (2006-2050) considering RCP 4.5 (left) and RCP 8.5 (right).The average curve is the average of four GCM models, and the variability curves are the minimum and maximum of four models.The red line is the slope of the time series during the projected period.The percent change from 2006 to 2050 is shown in red above each line.

Figure 11 .
Figure 11.Time series of annual AET based on four GCM models for the historical period (1980-2005) and projected period (2006-2050) considering RCP 4.5 (left) and RCP 8.5 (right).The average curve is the average of four GCM models, and the variability curves are the minimum and maximum of four models.The red line is the slope of the time series during the projected period.The percent change from 2006 to 2050 is shown in red above each line.

Figure 12 .
Figure 12.Time series of annual soil moisture based on four GCM models for the historical period (1980-2005) and projected period (2006-2050) considering RCP 4.5 (left) and RCP 8.5 (right).The average curve is the average of four GCM models, and the variability curves are the minimum and maximum of four models.The red line is the slope of the time series during the projected period.The percent change from 2006 to 2050 is shown in red above each line.

Figure 11 .
Figure 11.Time series of annual AET based on four GCM models for the historical period (1980-2005) and projected period (2006-2050) considering RCP 4.5 (left) and RCP 8.5 (right).The average curve is the average of four GCM models, and the variability curves are the minimum and maximum of four models.The red line is the slope of the time series during the projected period.The percent change from 2006 to 2050 is shown in red above each line.

Figure 12 .
Figure 12.Time series of annual soil moisture based on four GCM models for the historical period (1980-2005) and projected period (2006-2050) considering RCP 4.5 (left) and RCP 8.5 (right).The average curve is the average of four GCM models, and the variability curves are the minimum and maximum of four models.The red line is the slope of the time series during the projected period.The percent change from 2006 to 2050 is shown in red above each line.

Figure 12 .
Figure 12.Time series of annual soil moisture based on four GCM models for the historical period (1980-2005) and projected period (2006-2050) considering RCP 4.5 (left) and RCP 8.5 (right).The average curve is the average of four GCM models, and the variability curves are the minimum and maximum of four models.The red line is the slope of the time series during the projected period.The percent change from 2006 to 2050 is shown in red above each line.

Figure 13 .
Figure 13.Box plots of average annual deviations from the historical mean for two 20-year periods of 1996-2005 (historical: black) and 2031-2050 (projected RCP 4.5: green and projected RCP 8.5: Red).Boxes are the median, 25% and 75% percentiles, and whiskers extend to the most extreme data points.The historical mean value from 1996-2005 is used as the reference.

Figure 13 .
Figure 13.Box plots of average annual deviations from the historical mean for two 20-year periods of 1996-2005 (historical: black) and 2031-2050 (projected RCP 4.5: green and projected RCP 8.5: Red).Boxes are the median, 25% and 75% percentiles, and whiskers extend to the most extreme data points.The historical mean value from 1996-2005 is used as the reference.

Figure 14 .
Figure 14.Annual distribution in the ratios of AET/P and PET/P, and the Budyko curve for the watershed for the past period based on observed climate (1980-2011) as well as the historical period (1980-2005) and projected period (2006-2050) based on GCM climate annual values.

Figure 14 .
Figure 14.Annual distribution in the ratios of AET/P and PET/P, and the Budyko curve for the watershed for the past period based on observed climate (1980-2011) as well as the historical period (1980-2005) and projected period (2006-2050) based on GCM climate annual values.

Table 1 .
[51]ription and determination method of model parameters1.isthevariable infiltration curve parameter, D smax is the maximum base flow velocity, D s is the fraction of D smax where nonlinear base flow begins, W s is the fraction of maximum soil moisture content above which nonlinear base flow occurs, the Exp parameter characterizes the variation of saturated hydraulic conductivity with soil moisture, and the other two parameters are the depths of 2nd and 3rd soil layers; 2 Average slope of cell; 3 W s > D s ; 4 b 1 is slope of soil retention curve as used in Campbell's equation[51].

Table 2 .
Physical ranges for calibration parameters adopted in past studies.

Table 3 .
Streamflow calibration: Nash-Sutcliffe Efficiency (NSE) Coefficient and Relative Bias for the best combination of calibration parameters each Sub Drainage Area (SDA).

Table 3 .
Streamflow calibration: Nash-Sutcliffe Efficiency (NSE) Coefficient and Relative Bias for the best combination of calibration parameters each Sub Drainage Area (SDA).

Table 4 .
Statistical comparison of average annual fluxes for two 20-year periods of 1996-2005 (historical) and 2031-2050 (projected; RCP 4.5 and RCP 8.5).Results of 0 and 1 show, respectively, statistically insignificant and significant changes (significance level 0.01).The statistically significant changes are highlighted in bold.