Future shift of the relative roles of precipitation and temperature in controlling annual runoff in the conterminous United States

This study examines the relative roles of climatic variables in altering annual runoff in the conterminous United States (CONUS) in the 21st century, using a monthly ecohydrological model (the Water Supply Stress Index model, WaSSI) driven with historical records and future scenarios constructed from 20 Coupled Model Intercomparison Project Phase 5 (CMIP5) climate models. The results suggest that precipitation has been the primary control of runoff variation during the latest decades, but the role of temperature will outweigh that of precipitation in most regions if future climate change follows the projections of climate models instead of the historical tendencies. Besides these two key factors, increasing air humidity is projected to partially offset the additional evaporative demand caused by warming and consequently enhance runoff. Overall, the projections from 20 climate models suggest a high degree of consistency on the increasing trends in temperature, precipitation, and humidity, which will be the major climatic driving factors accounting for 43–50, 20–24, and 16–23 % of the runoff change, respectively. Spatially, while temperature rise is recognized as the largest contributor that suppresses runoff in most areas, precipitation is expected to be the dominant factor driving runoff to increase across the Pacific coast and the southwest. The combined effects of increasing humidity and precipitation may also surpass the detrimental effects of warming and result in a hydrologically wetter future in the east. However, severe runoff depletion is more likely to occur in the central CONUS as temperature effect prevails.

Surface and subsurface (shallow aquifers) runoff is a critical source of fresh water for humans (Vörösmarty et al., 2000). The impacts of temperature and precipitation changes on the magnitude and variability of runoff (Ficklin et al., 2009;Arnell and Gosling, 2013;Nash and Gleick, 1991;Vano et al., 2012) have drawn particular attention due to its importance for water supplies. Future changes in precipitation, evaporation, and plant water use are direct driving forces of runoff generation. Climate change alters both precipitation and the partitioning of precipitation into evapotranspiration (E T ) and runoff since a warmer climate generally provides more energy for water fluxes between the land and the atmosphere. Although an increase in precipitation may cause increase in both E T and runoff, the enhanced evaporative demand can result in decreases in runoff efficiency (ratio of runoff to precipitation) (McCabe and Wolock, 2016). Both observation and simulation studies in the US suggest that higher E T induced by rising temperature is unlikely to be counterbalanced by the increase in precipitation, and will lead to less runoff at large scales (Duan et al., 2016b;Jackson et al., 2005;Duan et al., 2017). Further, warming may also cause precipitation decrease in some regions and exacerbate the effects of temperature on runoff change.
Several studies have examined the relative contributions of historical changes in precipitation and temperature to runoff variation at watershed (Karl and Riebsame, 1989), regional (Ryberg et al., 2014;Gupta et al., 2015), and continental (McCabe and Wolock, 2011) levels across the conterminous US (CONUS). These studies all agree that precipitation, instead of temperature, explains most of the longterm change and variability in runoff during the past century. McCabe and Wolock (2011) suggested that the effects of temperature on runoff may become more substantial under a warming climate. However, no study in the literature has rigorously investigated the potential changes in the roles of precipitation and temperature under future climate scenarios. According to the Parameter-elevation Relationships on Independent Slopes Model (PRISM) dataset (http://prism.oregonstate.edu/) (Daly et al., 2008), the rate of decadal change in temperature over the CONUS fluctuated between −0.03 and +0.28 • C from the 1960s to the 2000s. The rate of warming is likely to accelerate under intermediate or high-emission scenarios and increase the pressure of water scarcity in many regions in this century (IPCC, 2014;Schewe et al., 2014). In addition, future change in climate is projected to vary spatiotemporally in both direction and magnitude in the CONUS (Mearns et al., 2012). Thus, sensitivity of water budget to climate change may be discrepant across time and space. Although the possible underestimation of the influence of temperature in altering regional water resources has been discussed recently (Sospedra-Alfonso et al., 2015;Woodhouse et al., 2016), a comprehensive evaluation under different climate backgrounds and land-cover compositions is still lacking.
We aim to address two questions in this study: (1) to what extent, if any, will the relative roles of precipitation and temperature in controlling runoff shift, if future climate changes follow the projections of climate models instead of the tendencies documented in the recent decades, and (2) how will runoff change in the future and what are the potential roles of other climatic driving forces besides precipitation and temperature? In the remainder of the paper, we first describe the methodology of runoff simulation and sensitivity assessment, and the hydro-climatic datasets used, followed by the results. Then, the advantages, limitations, and implications of this study are discussed and the conclusions are drawn.

Study area
The CONUS covers the 48 adjoining states and the District of Columbia. In the hydrologic unit system developed by the US Geological Survey (USGS) (http://water.usgs.gov/ GIS/huc.html), the nation is divided into six levels of hydrologic units and each unit is identified by a unique hydrologic unit code (HUC) consisting of 2-12 digits. The first level of classification divides the CONUS into 18 two-digit HUC areas that are also commonly referred to as water resource regions (WRRs) (Fig. 1). These regions can be further divided into 2099 eight-digit HUC areas, or HUC-8 watersheds. This study investigates climate and runoff variations at the resolution of HUC-8 watershed, as well as the aggregations in each WRR and the entire CONUS. The full lists and boundaries of hydrologic units at different levels can be found in the Watershed Boundary Dataset (https://datagateway.nrcs.usda.gov/).

Runoff modeling
The runoff responses to climate change and variability were modeled with the Water Supply Stress Index model (WaSSI). WaSSI is a monthly ecohydrological model that was developed to capture land-cover specific large-scale water balance in the CONUS based on empirical and physically based parameters (Caldwell et al., 2012;Sun et al., 2011b). It was integrated from a snow model, an E T model, and a soil moisture accounting model. A conceptual snow sub-model (Mc-Cabe and Markstrom, 2007) is used to partition the total precipitation into rainfall and snowfall, and to estimate snowpack melt/accumulation and snow water equivalent with concern of the mean elevation, latitude, and air temperature in the watershed. E T is calculated with an ecosystem E T model developed from the empirical relationships between E T and precipitation, potential evapotranspiration (PET), and leaf area index (LAI) (Sun et al., 2011a, b). These E T functions were established for 10 different land-cover classes independently to account for the different water demand within different vegetation, ranging from cropland, deciduous forest, evergreen forest, mixed forest, grassland, shrubland, wetland, open water, urban area, to barren land. Then, this E T estimation is further constrained by soil water availability, which is simulated using the algorithms of Sacramento Soil Moisture Accounting model (SAC-SMA) (Burnash, 1995), as well as the processes of infiltration and runoff generation at monthly basis. SAC-SMA is a classic rainfall-runoff conceptual model that has been successfully used by the US National Weather Service (NWS) to issue river forecasts across the country for decades. Necessary inputs for WaSSI include monthly precipitation, air temperature, PET, LAI, and landcover composition. In this study, the spatial distribution of LAI and the 10 land-cover classes were assumed to be static over time. Monthly climate data were first scaled to watersheds by the area-weighted averages. All the water balance components were calculated independently for each landcover class within each watershed, and then were aggregated into monthly means. The model parameters were acquired from several previous studies, including (1) the parameters of snow sub-model (four parameters for each WRR) were estimated for each WRR by comparing regional monthly mean snow water equivalent to remotely sensed values from the Snow Data Assimilation System (McCabe and Markstrom, 2007;Caldwell et al., 2012); (2) the parameters of E T submodel (three empirical parameters for each land-cover type) were estimated by empirical relationships derived from eddy covariance or sapflow measurements at multiple sites (Sun et al., 2011a, b); and (3) SAC-SMA parameters (11 parameters for each watershed) used to drive the soil water balance sub-model were developed from soil physical characteristics documented by the State Soil Geographic Database (http: //soildatamart.nrcs.usda.gov) (Anderson et al., 2006;Koren et al., 2003).
The WaSSI model has been validated against observations at USGS gauged sites at the levels of both 8-digit (Caldwell et al., 2012) and 12-digit HUC watersheds (S. . We here verify the model performance at CONUS and WRR scales to complement to previous validations. The simulated annual runoff, driven by monthly precipitation and temperature from the PRISM dataset, was compared against the USGS measurements over the entire CONUS ( Fig. 2a and c) and in the 18 WRRs ( Fig. 2b and d) for the time period of 1961-2010. Despite a slight overestimation of the minima, WaSSI shows reliable accuracy in capturing annual runoff at both CONUS and WRR scales, with R-square statistic reaching 0.91 and 0.95, and root mean squared error (RMSE) limited to 29 and 55 mm yr −1 , respectively.

Quantifying the independent effects of climatic variables
Large-scale water balance can be described as runoff (R) equals precipitation (P ) minus E T and changes in soil moisture (S M ) and the hydrologically connected snowpack (S P ): While P is the primary water input, changing temperature (T ) and other climatic factors interact with each other and affect R by altering the melt/accumulation of snowpack and controlling E T with the constraints of vegetation and soil moisture.
Here we develop a simple approach of sensitivity test to examine the relative roles of climatic variables in R variation, as where R denotes the change in R, which equals the combined effect of variations in all the climatic variables (Ci, i = 1, 2, . . ., N ). R can be decomposed into the independent effects of each variable (E Ci ) and the effect of interactions between them (E Int ). From a pre-change period (t 1 ) to a post-change period (t 2 ), R is quantified by R change (%) driven by changes in all the variables, as the difference between R(C1 t 2 , . . ., Ci t 2 , . . ., CN t 2 ) and R(C1 t 1 , . . ., Ci t 1 , . . ., CN t 1 ), while E Ci is estimated by R change driven by changes in the variable Ci only, as the difference between R(C1 t 1 , . . ., Ci t 2 , . . ., CN t 1 ) and R(C1 t 1 , . . ., Ci t 1 , . . ., CN t 1 ). E Int is calculated as the R minus N i=1 E Ci , representing the changes in R that cannot be accounted for by the independents effects. Given that the changing climatic variables may cause either positive or negative effects on R, their contributions (%) are quantified by the relative weights, as  Livneh et al., 2013, available at http://maca.northwestknowledge.net/) were used to test the potential future changes in R. This dataset includes the CMIP5 experiments of "historical", Representative Concentration Pathways (RCP) 4.5, and RCP8.5, which correspond to the climate forcings (i.e., greenhouse gases emissions, aerosols, land use feedbacks, etc.) observed in the history and projected in a future with the radiative forcing reaching 4.5 and 8.5 W m −2 in 2100 (equivalent to 650 ppm and 1370 ppm CO 2 ), respectively (Moss et al., 2010;IPCC, 2014). The used climatic variables include monthly P , maximum and minimum T , solar radiation (R s ), wind speed (W s ), and specific humidity (Sh) spanning from 1950 to 2099 (Fig. 3).
To evaluate the R responses to various changes in future climates, we conducted four 30-year simulation experiments: (i) RCP4.5/2030s (S1 scenario), near-future 2020-2049 under RCP4.5; (ii) RCP4.5/2080s (S2), far-future 2070-2099 under RCP4.5; (iii) RCP8.5/2030s (S3), near-future 2020-2049 under RCP8.5; (iv) RCP8.5/2080s (S4), far-future 2070-2099 under RCP8.5. These four future scenarios cover two post-change time periods (2030s and 2080s) and are compared to the historical condition in 1970-1999 (1980s) that represents the baseline level. Traditional sensitivity test methods usually assume a fixed amount of change (Karl and Riebsame, 1989) or allow one (or more) of the variables to remain constant over time (McCabe and Wolock, 2011). In this study, the 30-year-long continuous climate series were used to examine the long-term patterns while implicitly incorporating the inter-and intra-annual variations. This large set of climate projections was collected to enable a robust quantification of the major uncertainties from GCM structure and emission scenario.

Estimation of potential evapotranspiration
Hamon's PET equation has been used for PET estimation in previous WaSSI simulations because it only requires mean temperature as input and has shown reliable correlation with actual E T in historical periods (Lu et al., 2005;Vörösmarty et al., 1998). Essentially, temperature-based methods perform well because T is correlated with radiation and humidity at monthly timescale (Sheffield et al., 2012). Such correlations are the physical bases of the empirical E T functions, through which variability in P , T , and LAI was able to explain the Table 1. List of the 20 climate models and the changes in mean annual precipitation and temperature over the conterminous United States (CONUS) from the baseline scenario (B) to future scenarios S1 (RCP4.5/2030s), S2 (RCP4.5/2080s), S3 (RCP8.5/2030s), and S4 (RCP8.5/2080s). main controls of evaporation and transpiration fluxes without including the radiative and aerodynamic variables. However, recent studies revealed that the bias in temperature-based methods could be amplified in future scenarios of global warming, leading to overestimation of PET and ultimately E T and the severity of land surface drying (Milly and Dunne, 2011;Sheffield et al., 2012). Penman-Monteith (PM) reference E T (Allen et al., 1998), as a commonly used alternative PET model, incorporates the effects of surface temperature, humidity, wind, and radiation, and is considered the most reliable PET approach where sufficient meteorological data exist (Kingston et al., 2009;Feng and Fu, 2013). In this case, using the Hamon equation would lead to 130 mm yr −1 larger PET increase from the baseline to RCP8.5/2080s than that using PM equation (Fig. 4). We assume that the PM PET projections are more reasonable because the effects of future changes in R s , W s , and Sh are included as well as T . In the remainder of this paper, we focus on analyzing the R changes and the independent effects of five climatic variables based on PM PET, i.e., P , T (including changes in maximum T , minimum T , and mean T that was estimated as the average of maximum and minimum), R s , W s , and Sh. Effects of P and T evaluated from simulations of Hamon PET will also be investigated to address the consistency and discrepancy caused by using different PET methods.

Projected changes in runoff
Changes in mean annual R under future climate change scenarios vary among HUC-8 watersheds (Fig. 5) and WRRs (Fig. 6) across the CONUS. Runoff depletion is projected to cover most part of the central CONUS across WRR7-WRR12, with largest decreases over 50 % found in the south of WRR10 (Missouri) under RCP8.5. Increases are mainly projected in the southwest, the north of Missouri, and regions along the Atlantic coast and Pacific coast. Extreme increases over 100 % are projected in several arid watersheds in WRR15 (Lower Colorado) and WRR16 (Great Basin). However, this may be caused by the inability of GCMs in reproducing the low P values in these extremely dry areas. Although the general spatial patterns appear to be similar in the four scenarios, there is an evident expansion of the areas showing either extreme increasing or decreasing trend from the 2030s to the 2080s under both RCP4.5 (Fig. 5a, b) and RCP8.5 (Fig. 5c, d) scenarios.
The large variability of regional changes in R (Fig. 6) indicates considerable uncertainties from GCM structure. In most cases, the uncertainty range is limited to −30-+30 %, showing both positive and negative changing signals. The distributions of the median lines and interquartile ranges (IQRs) suggest a hydrologically drier future in WRR7-12 and WRR14 (Upper Colorado), where consistent decreasing signal is found in all the scenarios. Increasing trend can be found in WRR1 (New England), WRR2 (Mid-Atlantic), WRR17 (Pacific Northwest), and WRR18 (California). Generally, the uncertainty ranges tend to increase from 2030s to 2080s under both RCPs, and reach a particularly high level under RCP8.5/2080s. There is a noticeable consistency in the pattern that the GCMs agree more on the simulations in 2030s while the uncertainty aggregates over time toward the 2080s, which implies the limitation of the state-of-the-art GCMs in predicting the farther future.

Independent effects of climate variables
The changes in R discussed above are under the combined impact of changing P , T , R s , W s , and Sh. The independent effects of these factors over the entire CONUS are illustrated in Fig. 7a, b. The P and T are clearly the two most influential factors, which are projected to cause divergent changes in R due to the increase in P (+15-31 mm yr −1 ) and T (+1.8-5.3 • ). The median values show that annual R under the independent P effect is expected to increase by 13 mm yr −1 (4 %) in 2030s and 24 mm yr −1 (8 %) in 2080s under RCP4.5, and by 21 (7 %) and 30 (10 %) mm yr −1 at the same time under RCP8.5. In contrast, the independent effects of T reach −32 (−11 %), −50 (−17 %), −34 (−12 %), and −80 (−28 %) mm yr −1 in the scenarios S1-S4. The negative effect of rising T is expected to exceed the positive effect of increasing P and lead to overall decrease in R. However, Sh, the third largest contributor, will enhance R by 3-12 % and largely offset the T effects. A significant increasing trend in Sh is projected under both RCP4.5 and RCP8.5 (Fig. 3e), which will suppress vapor pressure deficit and thus partially counterbalance the increasing evaporative demand caused by warming. Meanwhile, the effects of R s (slightly negative), W s (slightly positive), and interactions among the factors (Int) are relatively minimal (< 3 %), suggesting that the variations in T , P , and Sh can explain the major changes in R.
It is worth noticing that much larger uncertainty ranges can be found in the P effects. Compared to the highly consistent increases in T and Sh, the 20 GCMs constantly disagree on the changing direction of P . Under RCP8.5/2080s, the multimodel result of P effect ranges from −11 to 24 %, and the IQR also reaches the highest level (13 %). This indicates that uncertainty in P projection is still the largest contributor to the uncertainty in R simulations, especially in the far future.
We also compared these results with those evaluated based on Hamon PET (Fig. 7c), and found some similar features. The differences in independent effects of P and T between the two sets of results are mostly smaller than 5 %, and both results show that the T effect would be twice as large as the P effect at the CONUS scale. This suggests that the bias in PET model structure is not likely to turn over the relative importance of P and T effects as long as E T model is properly calibrated. However, the projected decreases in R (i.e., the "Total" effects) are obviously more severe when using Hamon PET because the positive effect of increasing humidity is not considered.  Table 2 summarizes the relative contributions of P and T to R change for the historical and future periods in 18 WRRs and the entire CONUS. Historical changes in P , T , and their effects on R were tested using PRISM climate data spanning from January 1960 to December 2010. Given the significant spatial and temporal variability in R trend across the CONUS (Mauget, 2003;Wolock, 2002, 2011;Gupta et al., 2015), a consistent breakpoint is statistically unavailable. We hereby took 1985 as the breakpoint year for all the watersheds and evaluated the multi-decadal mean changes from 1961-1985 (pre-change period) to 1986-2010 (post-change period). Although the selection of different breakpoints may cause certain deviations, the analysis can provide a comparable benchmark for exploring the shifts in future scenarios at a multi-decadal scale. Unsurprisingly, the results of these latest decades show the prevailing role of P in nearly all the regions, with WRR14 being the only exception. In the future periods (from baseline to S1-S4), however, results derived from both PM and Hamon PET suggest that the role of T rise will surpass P and become the largest driver in most of the regions (15-16 out of 18 WRRs) in the future. In contrast, a larger mean contribution of P can be occasionally found in the Atlantic coast (WRR1, 2), Pacific coast (WRR18), and the southwest (WRR12, 15). Considering that the inconsistency among GCMs may make the recognition of larger contributor dubious, we used Wilcoxon signed-rank test (Gibbons and Chakraborti, 2011) to assess the statistical significance of the difference between each pair of P and T contributions (i.e., 20 samples from the 20 GCMs). The test results reveal high agreement among GCMs on the prominent role of T across most regions (underlined in Table 2). At CONUS level, the mean contributions of P and T are projected to lie within 20-24 and 43-50 % using PM PET, and 33-40 and 55-62 % using Hamon PET, suggesting a similar shift in the relative importance of these two key driving factors. However, future changes in Sh, R s , and W s account for another 16-23, 2-7, and 1-4 % of R change respectively, and indirectly affect the attributions to P and T . For example, the R increase in WRR1 would be completely attributed to P increase if Sh was not considered, and thus lead to an overestimation of P contribution.

Spatial distribution of the major driving factors
To further investigate the spatial pattern of future climatic controls on annual R, we mapped the coverage of dominant driving factors (Fig. 8) and examined its consistency with the changing trend in R at watershed scale (Table 3). Judging by multi-model ensemble means, P and T are the largest driving factor in 10-22 and 68-89 % of the CONUS area, respectively. High consistency on their dominant roles (80 % or more of the 20 GCMs agree on the sign) can be found in Figure 8. Relative importance of P and T in affecting runoff change across the HUC-8 watersheds in the future scenarios of S1 (RCP4.5/2030s) (a), S2 (RCP4.5/2080s) (b), S3 (RCP8.5/2030s) (c), and S4 (RCP8.5/2080s) (d). The watersheds under larger influence of P and T are marked with blue and red colors, respectively. The dark colors denote the areas where 80 % or more of the 20 GCMs agree on the sign, while the light colors denote the results of ensemble average. 4-7 and 21-41 % of the CONUS, respectively. As P and T are projected to keep increasing, the coverages of P -and Tdominant areas are also expected to expand from the 2030s to the 2080s. A directional change suggests that rising T will become more influential in the east (WRR1-6), while P will prevail in more watersheds across the west (WRR13-18). Although the aggregated effect of Sh is quite close to that of P at large scales, it is only expected to play a dominant role in several watersheds (1 % in area) across the borders between WRR10 and WRR11 under RCP8.5/2080s.
The P -dominant areas that are mainly distributed in the southwest (WRR13,15) and Pacific coast (WRR17,18) show clear signals of increasing R, driven by the widespread increase in P . One the other hand, only two-thirds of the Tdominant areas coincide with the areas of decreasing R, covering a large part of the central CONUS (WRR7, 9, 10, 11) and a number of watersheds scattered in the northwest (WRR14, 16, 17). Although T is also identified as the most influential factor in the eastern regions WRR1-5, the combined effect of the other four factors, primarily P and Sh, is projected to exceed the T effect and lead to an increase in R.

Spatial patterns of future runoff change
This study characterizes and generalizes large-scale relationships among changing P , T , and R despite the large geographic differences. The coherence in the spatial dynamics of R trend and the corresponding climatic drivers shows a Table 2. Comparison of multi-model averaged contributions (%) of precipitation (P ) and temperature (T ) to changes in runoff in the 18 water resource regions (WRRs) and entire CONUS in the historical period  and future scenarios S1 (RCP4.5/2030s), S2 (RCP4.5/2080s), S3 (RCP8.5/2030s), and S4 (RCP8.5/2080s). A larger contributor identified by multi-model ensemble means is in bold, and a significantly larger contributor identified by Wilcoxon signed-rank test (at 5 % significance) is in italic.

WRR
Historical Projections based on PM PET Projections based on Hamon PET S1 S2 S3 S4 S1 S2 S3 S4   P  T  P  T  P  T  P  T  P  T  P  T  P  T  P  T Table 3. Cross-comparison of the areal proportions (%) with different dominant driving factors and changing directions of runoff (R) in the future scenarios S1 (RCP4.5/2030s), S2 (RCP4.5/2080s), S3 (RCP8.5/2030s), and S4 (RCP8.5/2080s). The areas where a variable is the largest driving factor identified by multi-model ensemble means is marked in parentheses. The areas where a variable is identified as the "dominant" factor are in bold. A climate variable is identified as the "dominant" one only when 80 % or more of the 20 GCMs agree that it is the largest driving factor of runoff change.
Scenario S1 S2 S3 S4 Precipitation R * 4 (10) 7 (17) 6 (15) 6 (21) rough pattern: T change dominates R decrease while P and Sh changes dominate R increase. However, it should be interpreted with limitations on timescale and underlying surface features. This pattern does not hold true in all the watersheds due to the nonlinear complexity of R response to climate change at various timescales, as well as the influence of other watershed characteristics (e.g., topography, land use, soil property). For example, slight decreases in annual P but increases in annual R are projected in south Texas due to the changes in intra-annual climate variability. The role of T may also become more positive in regions where water avail-ability is dominated by snow melting (Barnett et al., 2005;Lutz et al., 2014). In addition, local R can be affected by other factors, such as land-cover evolution and the direct effects of atmospheric composition on transpiration (Gedney et al., 2006;Zhang et al., 2001Zhang et al., , 2015.

The role of land cover and land use
Land cover, LAI, and soil are important controls on catchment water balance and R sensitivity to climate change (Zhang et al., 2001;Bosch and Hewlett, 1982;Cheng et al., 2014). This study specifically focused on evaluating the separate and combined effects of changing climates on R within a static land cover/land use. We did not consider the potential evolution of land cover and its interactions with water balance. We made no explicit tabulation of the impact of land cover/land use on the R responses to climate change, but we did incorporate it as a key factor by estimating E T with a set of functions of climate, LAI, and soil moisture capacity and deficit. Across the land-cover classes, the uncertainty ranges of independent contributions of P (13-30 %) and T (39-51 %) are relatively small compared to the ranges across WRRs (18-47 and 29-52 %). This may be because the discrepancy across different land covers is largely offset by the different climate backgrounds across the country. Evaluation of future land-cover change and its impact on R is beyond the scope of this study. However, our results imply that the potential impact of land-cover change might not be large enough to alter the relative significance of P and T in controlling future continental water availability.

Implications for water and land management
Our results have important implications for water and land management across the CONUS. Water resource planning may need to prepare different management strategies for areas facing contrasting future hydrological conditions. Additional water storage such as reservoirs may be needed in regions expecting more R, while inter-basin water transfer, improving water use efficiency, and other water conservation measures such as rain harvesting, and waste water recycling should be implemented for areas expecting water shortages. The vast croplands across central US are likely to be threatened by rising T and diminishing water availability for irrigation and food production. Adaptations in cropping systems and irrigation strategy are needed to secure food supply and increase resiliency to drought and changing climate (Challinor et al., 2014;Teixeira et al., 2013). The drier and hotter conditions may also result in increasing water stress, higher risks of tree insects and disease outbreaks, and catastrophic wildfires in forests (Dale et al., 2001) (e.g., national forests in WRR14, 16, 17) and grasslands (e.g., in WRR10-11). Innovative land management practices such as forest thinning and fuel management, irrigation, and planting drought-tolerant species are vital to minimize the potential risk and vulnerability to climate change and reduce the threats to ecosystems and society (G. Grant et al., 2013;Vose et al., 2016).

Uncertainties and caveats
Considerable uncertainty lies in the projection of future climate changes from the 20 GCMs. The uncertainty ranges under both RCP4.5 and RCP8.5 show significant expansions over time from the 2030s to the 2080s. In particular, the large uncertainty in predicting future P may substantially compromise the reliability in evaluating either R change or the roles of P and T (Karl and Riebsame, 1989;Piao et al., 2010). Although the results allow us to draw some conclusions on the general patterns, uncertainties are large and vary differently across space and time. There are certain limitations in this evaluation that should be noted when interpreting the results. First, we did not incorporate other sources of uncertainty, such as the methodology of downscaling (Duan and Mei, 2014a;Chen et al., 2011), and structure and parameters of hydrologic model (Jung et al., 2012). Although the selections of GCM and emission scenario are more likely to be the largest sources of uncertainty in hydro-climatic modeling (Kay et al., 2009;Wilby and Harris, 2006;Duan and Mei, 2014b), the other sources may also affect the results to different extents. The roles of uncertainties from different sources can be particularly equivocal when investigating seasonal/monthly variability and extreme events (Bosshard et al., 2013;Giuntoli et al., 2015;Bae et al., 2011;Kay et al., 2009). Second, we focused on the independent effects of potential climate changes, while assuming that the interrelationship among the meteorological variables and waterbalance components remains the same as in historical periods. In future studies, improved climate datasets and better representation of the physical mechanisms of climatic factors (e.g., radiation: Bohn et al., 2013;wind speed: McVicar et al., 2012) are needed to reduce uncertainties.

Conclusions
This study evaluates the relative roles of precipitation and air temperature, as well as solar radiation, wind speed, and air humidity, in altering annual runoff across the CONUS based on a large ensemble of simulations using data from both historical measurements and CMIP5 GCM projections. Despite the large uncertainty and spatial variability involved in the results, two robust conclusions can be drawn at the CONUS and regional scales on a multi-decadal basis. First, the role of temperature will outweigh that of precipitation in a continued warming future in the 21st century, in spite of the fact that precipitation has been the primary control of runoff variation during the latest decades. The projections from 20 climate models suggest a high degree of consistency on the increasing trends in both precipitation and temperature, but the negative effect of temperature is expected to exceed the positive effect of precipitation on runoff change in most regions. Over the entire CONUS, temperature is projected to be the largest contributor (43-50 %), followed by precipitation (20-24 %), humidity (16-23%), solar radiation (2-7 %), and wind speed (1-4 %). Spatially, precipitation is likely to be the dominant driving factor for runoff increase across the Pacific coast and the southwest, while temperature will be more influential in the central CONUS and parts of the northwest and cause runoff decreases. Second, increasing humidity is expected to partially offset the additional evaporative demand caused by warming, and consequently enhance runoff widely across the country. Although the rising temperature is projected to be the largest control of runoff change in the eastern CONUS, the combined effects of increasing humidity and precipitation will surpass the detrimental effects of warming and result in a hydrologically wetter future. This study also raises concern on the choice of PET method. It has been well acknowledged in hydrometeorological communities that temperature-based PET methods tend to be oversensitive to temperature change.
Our results further demonstrate that the main risk of using temperature-based PET is overlooking the effects of other changing climatic variables (mainly humidity in this case), which have not been as widely measured as temperature and are relatively understudied, rather than overestimating the effects of temperature.
Data availability. All data used in the analysis and conclusions in this paper are available in the references.
Author contributions. KD and GS designed the study. KD performed the analysis. EC, HA, SS, DC, and XZ helped collect the data. All authors contributed to the interpretation of the results. KD wrote the initial draft, and GS, PC, SM, and YZ refined the paper.
Competing interests. The authors declare that they have no conflict of interest.