Projecting water yield and ecosystem productivity across the United States by linking an ecohydrological model to WRF dynamically downscaled climate data

Quantifying the potential impacts of climate change on water yield and ecosystem productivity is essential to developing sound watershed restoration plans, and ecosystem adaptation and mitigation strategies. This study links an ecohydrological model (Water Supply and Stress Index, WaSSI) with WRF (Weather Research and Forecasting Model) using dynamically downscaled climate data of the HadCM3 model under the IPCC SRES A2 emission scenario. We evaluated the future (2031–2060) changes in evapotranspiration (ET), water yield (Q) and gross primary productivity (GPP) from the baseline period of 1979–2007 across the 82 773 watersheds (12-digit Hydrologic Unit Code level) in the coterminous US (CONUS). Across the CONUS, the future multi-year means show increases in annual precipitation (P ) of 45 mm yr (6 %), 1.8 C increase in temperature (T ), 37 mm yr (7 %) increase in ET, 9 mm yr (3 %) increase in Q, and 106 gC m yr (9 %) increase in GPP. We found a large spatial variability in response to climate change across the CONUS 12-digit HUC watersheds, but in general, the majority would see consistent increases all variables evaluated. Over half of the watersheds, mostly found in the northeast and the southern part of the southwest, would see an increase in annual Q (> 100 mm yr or 20 %). In addition, we also evaluated the future annual and monthly changes of hydrology and ecosystem productivity for the 18 Water Resource Regions (WRRs) or two-digit HUCs. The study provides an integrated method and example for comprehensive assessment of the potential impacts of climate change on watershed water balances and ecosystem productivity at high spatial and temporal resolutions. Results may be useful for policy-makers and land managers to formulate appropriate watershed-specific strategies for sustaining water and carbon sources in the face of climate change.


Introduction
Due to human activities, such as emissions of greenhouse gas, aerosol, and land use/cover change (LUCC), the Earth's climate system has been significantly altered over the past 100 years.The Intergovernmental Panel on Climate Change (IPCC, 2014) concludes that global surface temperature has increased 0.85 • C during 1880-2012, and increased 0.78 • C  during 2003-2012 when compared to 1850-1900.Additionally, extreme precipitation and droughts have increased (Tebaldi et al., 2006;Trenberth, 2011;Bony et al., 2013;Hegerl et al., 2014).The global climate is projected to continue to change over this century and beyond (IPCC, 2014).In comparison to the period of 1986-2005, the period 2018-2100 is projected to see 0.3 to 4.8 • C increase in global surface temperature (IPCC, 2014).Future changes in precipitation show a small increase in the global average, but a substantial shift in where and how precipitation falls (Noake et al., 2012;Scheff and Frierson, 2012;Liu et al., 2013a).
In the US, average temperature has dramatically increased since the record keeping began in 1895.The most recent decade was believed to be the warmest on record (see the website: http://www.nasa.gov/home/hqnews/2010/jan/HQ_10-017_Warmest_temps.html).Mean precipitation over the US has increased overall since 1900; some areas have increased with a higher rate than the national average, and some areas have decreased (Groisman et al., 2004;Meehl et al., 2005;Anderson et al., 2015).Over the past century, climate change in the US has caused severe water stress, floods and droughts as well as forest morality (Xu et al., 2013), leading to serious economic losses in some regions.Quantifying the impacts on future climate change on water and ecosystem productivity has become a major research area in hydrology and ecosystem sciences (Lettenmaier et al., 1994;Lins and Slack, 1999;Groisman et al., 2001;McCabe and Wolock, 2011;Sagarika et al., 2014).
Because climate change patterns are not uniform across space or time (Sankarasubramanian et al., 2001;Sankarasubramanian, 2003;Wang and Hejazi, 2011;Xu et al., 2013;Brikowski, 2014) climate change impacts on water cycle and ecosystem productivity vary from region to region, and variability will be even bigger across small watersheds.To support future water resource planning, watershed management and to develop sound adaptation strategies over the continen-tal US (CONUS), tools are needed to integrate various climate scenarios from a variety of atmospheric ocean general circulation models (AOGCMs) and community Earth system models (CESMs), and hydrological and vegetation dynamic models (Brown et al., 2013;Blanc et al., 2014;Yu et al., 2014).
Two major research gaps exist in past climate change studies that aim at quantifying the interactions among climate, hydrology and ecosystem productivity.First, few studies have provided projections of future climate change impacts on water and carbon balances at watershed scale using a consistent approach.Various land surface models (LSMs) simulate and predict water fluxes for a large region, but the scale is often too coarse with a spatial resolution ranging from 0.25 to 2.5 • .The water budget within each grid cell in LSMs may not be balanced because it is not a closed watershed system.Key hydrological processes (e.g., lateral surface and sub-surface flows among grid boxes) have been rarely considered, potentially resulting in uncertainties in water balance projections (Overgaard et al., 2006;Li et al., 2011).Second, to save computational resource and enhance the computational efficiency, a statistical (or empirical) downscaling method has been mostly used to generate climate forcing to land surface models or watershed ecosystem models.However, this type of method does not consider the effects of atmospheric dynamical processes (Xue et al., 2014) and could introduce uncertainties into the crucial land surface variables.
Therefore, the general goal of this study is to explore how dynamically downscaled climate data can be used to drive a common ecosystem model for climate change assessment at a fine spatial scale (i.e., 12-digit HUC watersheds, whose detailed information can be found in the following text).The specific objectives of this study are to (1) evaluate future climate changes in precipitation and temperature during 1979-2007 and 2031-2060 for one emission scenario over the CONUS using dynamically downscaled climate projections from the WRF (Weather Research and Forecasting) model; (2) project future changes of water yield (Q), ET, and GPP for the study area by linking the WRF dynamically downscaled climate change scenarios and the WaSSI model.The goal is to generate information that can be useful for policy makers to plan for potential shifts in water resources and ecosystem productivity at the watershed to national level.

Study area
The research area includes the conterminous continental US covering 82 773 12-digit HUC watersheds within the 18 Water Resources Regions (WRRs; Fig. 1a).The size of these HUC12 watersheds ranges from 0.16 to 9238 km 2 , with the median and the mean values of 88.2 and 95.0 km 2 , respectively.Moreover, the area of the overwhelming majority of the watersheds (> 80 000) is between 50 and 170 km 2 .The WRRs vary in size with the maximum of 1.3 × 10 6 km 2 (WRR10) and the minimum of 1.1 × 10 5 km 2 (WRR6).In addition, climatology and land surface characteristics (e.g., land cover; Fig. 1b) vary dramatically among these WRRs.From the east to the west CONUS, multiyear mean  annual precipitation as estimated by the Parameter-elevation Regressions on Independent Slopes Model (PRISM) shows longitudinal decreases ranging from 1300 mm yr −1 to 341 mm yr −1 .For the multi-year mean temperature , the spatial distribution displays the latitudinal characteristic decreasing from the south to the north CONUS, with a range from a high of 18 • C to a low of 4 • C. The WRRs in the east had the larger percentages (around 10 %) of urban use with WRR2 (13 %) and WRR4 (11 %) ranked as the top two.The wetlands are mainly located in the WRRs in the eastern US, while the western regions had the higher percentages of shrubland (> 30 %).The WRRs in the east generally had higher forest (including mixed, evergreen and deciduous forests) percentages (> 33 %) than the southwest (< 30 %).The deciduous and the evergreen forests were mainly found in the east and the west, respectively. Mos of the crop lands were located in the east and central CONUS (Fig. 1b).

Dynamically downscaled climate by WRF
The IPCC Special Report on Emissions Scenarios (SRES) scenarios were designed to project future global environment with a special reference to the production of greenhouse gases and aerosol precursor emissions (Nakicenvoic and Swart, 2000).The SRES scenarios mainly include four narrative storylines (i.e., A1, A2, B1, and B2), which describe the relationships between the forces affecting greenhouse gas and aerosol emissions and their evolution in the 21st century for large regions and the globe.Each storyline represents a specific and typical demographic, economic, technological, social, and environment progresses with di-vergence in increasingly irreversible ways.The A2 storyline represents the high end of the SRES emission scenarios (but not the highest) and has been widely used by the scientific communities (Seneviratne et al., 2006;Wi et al., 2012).Therefore, the SRES A2 emission scenario was selected in this study.From an impact and adaptation point of view, if one can adapt to a larger climate change, then the smaller climate changes of the lower-end scenarios can also be adapted to.Moreover, the historic emissions (1990 to present) correspond to a relatively high emission trajectory (http://www.narccap.ucar.edu/about/emissions.html).
The global circulation models (GCMs) have significant issues in representing local climates, mountains in particular, because of their coarse spatial resolution (Leung and Qian, 2003).To downscale the GCMs climate data to a higher spatial resolution for regional and local applications, two types of downscaling method are available: dynamical and statistical (or empirical) downscaling (Huang et al., 2011).Due to better representation of finer-scale physical processes in climate variables (Gao et al., 2011;Xue et al., 2014), dynamical downscaling was used here for generating the current and the future climate.
The HadCM3 (Hadley Centre Coupled Model, Version 3) is a coupled atmosphere-ocean general circulation model (AOGCM) developed by the Hadley Centre in the United Kingdom (Gordon et al., 2000;Pope et al., 2000;Collins et al., 2001), which has been used extensively for climate prediction, detection, and attribution, and other climate sensitivity studies, e.g., the Third, the Fourth, and the Fifth IPCC Assessments reports.For the atmospheric component, the model dynamics and physics are solved on a 3.75 • (longitude) × 2.5 • (latitude) grid with 19 hybrid vertical levels, and a horizontal resolution of 1.25 • (longitude) × 1.25 • (latitude) with 20 vertical levels in the oceans.The reader is referred to Pope et al. (2000) for details of the HadCM3 dynamical and physical processes.Generally speaking, although flux adjustments are not utilized by the HadCM3, it still ranks highly compared to other models in the respect of current climate simulation (Reichler and Kim, 2008).In addition, among the many GCMs, the HadCM3 model was believed to have the most realistic description of the ENSO mechanisms in the current climate, and reasonably capture ENSO-associated precipitation anomalies over North America (van Oldenborgh et al., 2005;Joseph and Nigam, 2006;Dominguez et al., 2009).Based on the importance of precipitation in hydrology and ecosystem productivity assessment, we chose the HadCM3 model to provide forcing fields for running the Advanced Research version (ARW) of the Weather Research and Forecasting (WRF) regional climate model (Skamarock et al., 2005).
Data generated from the WRF model are described below.The WRF model was run for the years 1969 to 2079 at a 35 km resolution.HadCM3 inputs with 6 h time resolution were used, and the dynamically downscaled output by the WRF model was also stored at 6 h time intervals.For the model domain, the CONUS and northern Mexico were included (Wi et al., 2012).The model's physical parameterizations mainly included: WRF Single-Moment three-class microphysics (Hong et al., 2004), Kain-Fritsch cumulus parameterization (Kain and Fritsch, 1993), Goddard Shortwave radiation (Chou and Suarez, 1994), Rapid Radiative Transfer Model (RRTM), Longwave (Mlawer et al., 1997), Eta surface layer (Janjic, 2002), Mellor-Yamada-Janjic (MYJ) planetary boundary layer (Janjic, 2002), and the Noah land surface model Version 1.0 (Chen and Dudhia, 2001).To ensure the maintenance of synoptic-scale circulation features, like ridges and troughs, in the RCM (regional climate model), we performed spectral nudging on the zonal and meridional winds, the temperatures and the geo-potential height fields for all pressure levels below 0.36 of the surface pressure (for a surface pressure of 1000 mb this would mean all pressures below 360 mb) effectively nudging only at very high elevations above the surface.

Climate data bias corrections
The dynamically downscaled precipitation and temperature simulations by WRF were sufficient for a hydrological study  by Wi et al. (2012) in the Colorado River Basin).Our comparison study showed that although downscaled climate simulations agreed well with the observations (PRMS data) in a climatological sense, some large regional biases were found.Therefore, bias correction was performed using a monthly Bias Correction Spatial Disaggregation (BCSD; Wood et al., 2002Wood et al., , 2004) ) approach.The method has been applied for hydrologic forecasting in the eastern US (Wood et al., 2002).Basically, the bias correction includes the following procedures: (1) scale up the PRISM monthly precipitation and temperature with 4 km × 4 km resolution to match the simulated WRF data (35 km × 35 km) for the time period of 1978-2007; (2) construct cumulative distribution functions (CDFs) for climate variables in each grid cell and month for both historic WRF and upscaled PRISM data sets; (3) the paired CDFs are combined to form a "quantile map", where at each rank probability or percentile, the bias between the WRF and the PRISM (at that location, for that variable, and during that month) is calculated; (4) the computed bias in each month, grid cell and variable are applied to the WRF future outputs .The detailed procedures can be found in Brekke et al. (2013), see; http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/techmemo/downscaled_climate.pdf.Both the corrected WRF monthly precipitation and temperature in historic and future periods were scaled to the 12-digit HUC watershed scale because the WaSSI model operated on the 12-digit HUC watershed level.

The WaSSI model
The WaSSI model is an integrated, water-centric processbased ecohydrological model designed for modeling water and carbon balance and water supply stress at a broad scale (Sun et al., 2011a;Caldwell et al., 2012;Sun et al., 2015a, b).It operates on a monthly time step at the 8-digit HUC or 12-digit HUC watershed scale for the CONUS.The WaSSI model simulates the full monthly water (ET, Q, and soil moisture storage) and carbon balances (GPP, ecosystem respiration, and net ecosystem productivity) for each land cover class at the given watershed scale.This model has been tested in a variety of geographical regions, and widely used for quantitatively assessing combined or individual effects of climate change, land use/cover change (LUCC), and population dynamics on water supply stress and ecosystem productivity (i.e., carbon dynamic) over the CONUS (Sun et al., 2008(Sun et al., , 2011a;;Lockaby et al., 2011;Caldwell et al., 2012;Averyt et al., 2013;Tavernia et al., 2013;Marion et al., 2014;Sun et al., 2015a, b).The model has also been applied internationally in Mexico, China (Liu et al., 2013b), andAfrica (McNulty et al., 2016).
The key algorithms of the WaSSI model were derived from accumulated knowledge of ecosystem carbon and water cycles gained through the global eddy covariance flux monitoring networks and watershed-based ecohydrological studies across the US.The ecosystem ET sub-module, the core of the WaSSI model, is described as a function of potential ET (PET), LAI, precipitation, and soil water availability by land cover type (Sun et al., 2011a).The snow model embedded within WaSSI (McCabe and Wolock, 1999;McCabe and Markstrom, 2007) estimates snow melt rates, mean monthly snow water equivalent (SWE), mean watershed elevation, and monthly air temperature.Infiltration, surface runoff, soil moisture, and baseflow processes for each watershed are simulated by the Sacramento Soil Moisture Accounting Model (SAC-SMA; Burnash, 1995).The ecosystem productivity module computes carbon dynamics (GPP and respiration) using linear relationships between ET and GPP derived from global eddy covari-ance flux measurements (Sun et al., 2011a, b).The User Guide of WaSSI Ecosystem Services Model, Version 2.1 (www.forestthreats.org/research/tools/WaSSI)provides detailed description of model algorithms and data requirements (Caldwell et al., 2012).
The WaSSI has been evaluated at multiple scales using gaging station data for streamflow and remote sensing products for evapotranspiration across the US (Sun et al., 2011a;Caldwell et al., 2012;Sun et al., 2015a).At the 12-digit HUC scale, the model was validated using monthly and annual water yield data collected at 72 selected USGS watersheds, and ET and GPP data for 170 National Forests over the CONUS (Sun et al., 2015a).Overall, the validation results suggested that this model could capture characteristics of water and carbon balances at the selected spatial levels under various climatic conditions (Sun et al., 2015a, b).

Impact analysis
We first examined modeled changes in monthly ET and GPP at the 12-digit HUC watershed scale using the WRF dynamically downscaled, bias-corrected historic and future climate data, respectively.Then, we computed future annual changes at three spatial levels: the entire CONUS as a whole, the 12digit HUC watershed, and the individual WRR.The multiyear means of annual precipitation, temperature, ET, Q, and GPP averaged across the whole CONUS, WRR, or each 12digit HUC watershed for the 1979-2007 time period were compared to those for the 2031-2060 period.The absolute or percentage (except for temperature) changes for each variable were calculated.Herein, the absolute differences were expressed as the future means minus those in the historical period, while the percentage differences were calculated using the absolute difference divided by baseline mean in the 1979-2007.In addition, the future monthly changes of these ecosystem flux variables were also assessed for the whole CONUS and each WRR.

Baseline characteristics of hydro-climatology and ecosystem productivity (1979-2007)
For the baseline period, multi-year means of annual precipitation (Fig. 2a), ET (Fig. 3a), Q (Fig. 3e), and GPP (Fig. 4a) all generally showed longitudinal decreases from east to west across the CONUS.The Pacific Northwest region has the highest precipitation (> 1800 mm yr −1 ), followed by the larger values for precipitation in the southeast (> 1200 mm yr −1 in Fig. 2a).For ET, the maximum (> 750 mm yr −1 in Fig. 3a) mainly appeared in the southeast.
The largest Q higher than 600 mm yr −1 (Fig. 3e), mainly existed in the Pacific Northwest region, the Rocky Mountains, and the Appalachian Mountains, especially for some 12digit HUC watersheds in the Pacific Northwest region being greater than 1000 mm yr −1 .For GPP (Fig. 4a), the 12-digit HUC watersheds with higher values (> 1000 gC m −2 yr −1 ) were mainly located in the areas of the southeast and the Pacific Northwest.By contrast, the average annual temperature climatology of the CONUS presented a clear latitudinal increase ranging from −0.8 • C in the north to 22 • C in the south.Because of topographical effects, temperature in the Rocky Mountains was lower than 4 • C relative to the surrounding regions.
The baseline intra-annual precipitation presented a complicated pattern (Fig. 5).Except in February, precipitation in all the months was more than 65 mm yr −1 , and peaked in May with 78 mm yr −1 .Overall, temperature (Fig. 5b), ET (Fig. 5c), and GPP (Fig. 5e) all increased gradually starting from January, peaked (24.8 • C, 80 mm yr −1 and 205 gC m −2 yr −1 , respectively) in July, and then decreased sharply.Fluctuations of Q clearly differed from other variables (Fig. 5d) following a pattern similar to a sine function.Q increased in January, peaked in April (36 mm yr −1 ), decreased to the lowest (15 mm yr −1 ) in August, and then rose.
We also explored multi-year mean monthly precipitation, temperature, ET, Q, and GPP for each WRR (not shown).Generally, the intra-annual distribution was different (e.g., phases and magnitudes) among the 18 WRRs, due to the complex differences in topography and climate among them.For WRR16-18, most precipitation fell in January-April and October-December, while precipitation in other WRRs mainly concentrated in May-September.In all the WRRs, the intra-annual temperature followed a unimodal curve, with peaks in July or August and the lowest values in January or December.For ET and GPP, the higher values were mainly found from May to November, except for WRR18.Com- paring the monthly distributions among the 18 WRRs, they could be divided into three categories: unimodal, sine, and trough curves.

Future climate change
Future precipitation and temperature followed a similar pattern as the baseline (Fig. 2).Precipitation showed a longitudinal decrease from the east to the west, but temperature presented a clear latitudinal decrease.However, for each 12-digit HUC watershed, these two climate variables would increase or decrease by different magnitudes in the future (Fig. 2c and d for precipitation, and Fig. 2g).During 2031-2060, annual precipitation would increase in 82 % of the CONUS 12-digit HUC watersheds, while decreasing in the rest of the watersheds that were mainly located in the southeast and the west coastal regions.The northeast and the northwest coastal regions would generally have a greater increase (> 150 mm yr −1 ) or decrease (> 200 mm yr −1 ), respectively, in P (Fig. 2c).The greater percent increases in  precipitation (> 18 %) were found in some watersheds in the southwest and the northeast regions (Fig. 2d).Future temperature would increase consistently across watersheds, ranging  be 844 mm yr −1 and 13.1 • C, respectively (Table 1).The mean annual P for the entire CONUS would increase by 45 mm yr −1 (6 %) and T increase by 1.8 • C , respectively (Table 2).Except for the WRR17 with a slight decrease in P (13 mm yr −1 or 1 %), the other 17 WRRs all exhibited increases.A large absolute increment of precipitation (> 60 mm yr −1 ) could be found in WRR2, 4, 5, and 7, while the WRR8 and 14 have lower increases (< 15 mm yr −1 ).For the percentage increment, the higher increases in precipitation (≥ 10 %) existed in the WRR2, 5, 15, and 16; however, WRR1 and 8 showed lower increases (≤ 1 %).For future temperature, all the 18 WRRs show increase relative to the past period, especially WRR9, 10, 14, and 16 (≥ 2 • C).
Both future P and T had similar intra-annual fluctuations to those of the baseline period (top panels in Fig. 5a and b).However, the magnitudes of differences in both P and T differed in different seasons (bottom of Fig. 5a and b).In most months, precipitation would increase ranging from 3 to 11 mm yr −1 , especially in January, May, and September (> 7 mm yr −1 ).For February, March, October, and November, P would have slight reduction with a range from −5 to −1 mm yr −1 .The temperatures for each month would significantly increase by at least 1.5 • C, particularly for January and June-October (> 2.0 • C) (Fig. 5b).
The comparisons of seasonal climatic change patterns among the 18 WRRs suggested that timings agreed well among WRRs (not shown).However, the magnitudes of changes varied greatly.The future monthly precipitation shows an increase in January and May-October in more than 10 WRRs.The differences are most pronounced for January, July, and September (Fig. 6a).In other months, however, future monthly precipitation would reduce to some extent in most of the WRRs.The future monthly temperature for all the WRRs is shown to increase in a range from 0.5 to 3.0 • C. In January and June-October, temperatures in most WRRs show increase at a relatively high rate (> 1.5 • C) compared to other months for most WRRs.

Annual change
The spatial patterns in ET and Q for the baseline are similar to those in the future (Fig. 3).However, the changes of annual ET (Fig. 3c and d) and Q (Fig. 3g and h) for each 12-digit HUC watershed vary in magnitude spatially.Overwhelmingly, the majority (98 %) of the CONUS 12-digit HUC watersheds are suggested to increase in annual ET, and the watersheds with annual ET reduction are mainly concentrated in the northwest coastal region.Regarding the absolute difference of ET (Fig. 3c and d), annual ET shows a relatively higher increase (> 32 mm yr −1 ) in the northeast CONUS, especially in the southeast coastal region and the south part of the northeast CONUS (> 48 mm yr −1 ).Different from the absolute changes, relative changes ( %) in most of the western regions (excluding the west coast) and the northeast had high values (> 6 %) with the highest increments (> 12 %) found in the south of the southwest CONUS.
Across the CONUS, annual Q in 52 and 48 % of the CONUS 12-digit HUC watersheds is projected to increase and decrease by 2031-2060, respectively (Fig. 3g and 3h).In general, the northeast and the south part of the south CONUS shows increase in annual Q, while other regions show decrease (Fig, 3g and h).The positive (> 100 mm yr −1 ) and the negative (> 100 mm yr −1 ) changes in Q are mainly in the northeast, and the west coastal and the southeast regions, respectively.Q in the south part of the southwest CONUS is projected to significantly increase (> 20 %), while the central part of the west CONUS generally decreases by more than 20 %.

Seasonal change
The variations of future CONUS-wide multi-year mean monthly ET and Q are presented in Fig. 5c and d.Although these two variables had similar intra-annual fluctuations to those of the baseline period, their monthly magnitudes changed to some degree.Overall, the future monthly ET would increase with the largest increments (> 2 mm month −1 ) in January.April-October had higher values than the other 5 months.For monthly Q, most of the months (9 months) would increase, especially in January and September (increase > 3 mm month −1 ).
We also compared the future intra-annual fluctuations of ET and Q to those of the baseline period, and found that each WRR agreed well in their flow timings for the baseline and the future periods (not shown here).Fig. 6c and d present the number of the WRR within a given difference interval for ET or Q by month respectively.Generally, future monthly ET would increase by different rates for each month at each WRR (Fig. 6c).Moreover, ET from May to September (roughly the growing season) would have greater increments (> 2.4 mm yr −1 ) in most of the 18 WRRs.Q in most of WRRs would increase in January, February, July, September, and December, but would decrease in April and November.

Annual change
The overall spatial distribution of GPP did not change in the future (Fig. 4b) when compared to the baseline (Fig. 4a).For each 12-digit HUC watershed, GPP would change with great spatial variations (Fig. 4c and d).In the future, the overwhelming majority (98 %) of the CONUS 12-digit HUC watersheds would increase in annual GPP.The watersheds with annual GPP reduction were mainly located in the northwest coastal region.A relatively high increase (> 120 gC m −2 yr −1 ) were found in the northeast, especially in the south part of the region (> 180 gC m −2 yr −1 ; Fig. 4c).In contrast to the absolute difference, most of the west CONUS (excluding the coastal regions) show great increase (> 12 %) in relative change (%) of annual GPP.The highest changes (> 20 %) were mainly located in the south of the southwest region.

Seasonal change
Figure 5e shows the future multi-year mean monthly GPP averaged over the whole CONUS.Despite the similar intraannual fluctuations of multi-year mean monthly GPP during the baseline and the future periods, the future magnitude in each month is projected to change to some degree (the bottom of Fig. 5e).Overall, the future monthly ET would have larger increments (> 9 gC m −2 yr −1 ) in January and May-October than other months.The future intra-annual fluctuation patterns of GPP for each WRR are similar to the baseline periods (not shown here).As indicated by the number of the WRR within a given GPP difference interval (Fig. 6e), the future monthly GPP generally would increase by different rates for each WRR.Moreover, GPP from May to September would have greater increments (> 4 gC m −2 yr −1 ) in most of the 18 WRRs.

Uncertainties
In the present study, we have assumed that the water balance and ecosystems at each 12-digit HUC watershed are unaffected by human activities as represented by a fixed land cover (year 2000), and ecosystem flux changes are attributed to climate change alone.However, one way or another, most catchments in the US have experienced some level of human influence (National Research Council, 2002).Hydrology and ecosystems can be influenced significantly by human activities on various temporal and spatial scales (Foley et al., 2005;Harding et al., 2011).Hydraulic projects such as dam constructions, reservoir management (Hu et al., 2008), groundwater withdrawals for irrigation and domestic use, and land use/cover change all affect watershed balances (Foley et al., 2005;Piao et al., 2007;Wang and Hejazi, 2011;Schilling et al., 2008) and ecosystem productivity (Zhang, Y. et al., 2014).
Similarly, natural disturbances (e.g., wildfire, climate extremes, and pest and pathogen outbreak) also impact water balance and ecosystem productivity -in both the past and the future.For example, the direct effects of wildfire include plant mortality and thus exert adverse impacts on vegetation productivity, consequently leading to a decrease in carbon uptake and stocks (Lenihan et al., 2008;Dore et al., 2010;Lee et al., 2015).Wildfires alter the watershed hydrologic processes through reducing vegetation canopy interception, transpiration, and infiltration rate (Yao, 2003;Neary et al., 2005;Bond-Lamberty et al., 2009;Brookhouse et al., 2013;Nolan et al., 2014Nolan et al., , 2015)).As an important natural disturbance, droughts generally increase vapor pressure gradient between leaves and atmosphere and thus cause stress on plant hydraulic systems (Anderegg et al., 2012;Reichstein et al., 2013).As a result, high tension in the xylem can trigger embolism and partial failure of hydraulic transport in the stem, and even tend to result in vegetation mortality, which can aversely impact on water yield and carbon sink capability (Cook et al., 2007;Allen et al., 2010;Guardiola-Claramonte et al., 2011;Adams et al., 2012).Droughts often lead to pest and pathogen outbreaks (Overpeck et al., 1990;Hason and Weltzin, 2000;Marengo et al., 2008;DeRose and Long, 2012;Jactel et al., 2012), and thus predispose an individual plant species to disease or mortality (Schoeneweiss, 1981;Ayres and Lombarder, 2000).Although our modeling approach has considered water stress on productivity, tree mortality has not been dealt with, so the impacts of droughts on GPP might be underestimated, and water yield may be underestimated as well.
Additionally, elevated CO 2 and climate change can also execrate impacts on hydrological and ecosystem productivity through changing water use efficiency (Miller-Rushing et al., 2009;de Kauwe et al., 2013;Zhang, F. et al., 2014;Liu et al., 2015) and vegetation processes (e.g., stomatal conductance and LAI; Sun et al., 2014).However, the WaSSI model did not consider these effects, potentially resulting in errors in estimating ET, GPP, or water yield (Cox et al., 2000;Gedney et al., 2006;Oki and Kanae, 2006;Betts et al., 2007;Piao et al., 2007).Human activities aside, natural disturbances and their couplings may introduce uncertainties into our results.However, the potential errors are largely dependent on specific trajectories of climate change and land cover change (Qi et al., 2009;Thompson et al., 2011;Alkama et al., 2013).The complex interactions of climate, disturbance, and ecohydrological processes require a more mechanistic integrated modeling approach.

Land management implications
Numerous modeling studies around the world have shown that future climate change could increase or decrease the water availability to certain specific ecosystems and human populations under different climate scenarios (Arnell, 1999;Blanc et al., 2014;Ingjerd et al., 2014;Kundzewicz and Gerten, 2015).Our analyses have shown that, over the whole CONUS, P would increase by 45 mm (6 %) leading to a small increase in Q of 9 mm yr −1 (3 %).So, climate change under the SRES A2 scenario has little influence on water shortage for the entire CONUS.However, there are large regional differences in Q responses to future climate change among the 18 WRRs.The magnitude is large, from a de-crease of −32 mm yr −1 to an increase of 113 mm yr −1 or from −12 to 21 %.Despite the increase in annual P , annual Q in the WRR1, 3, 8, 11, 14, 16 and 18 decreased by various degrees, due to the increased ET.Consequently, the climate scenario studied will likely increase stress on the water supply in these WRRs.In addition, it is worth noting that monthly responses of Q to future climate also vary among watersheds.Water yield in about half of the 18 WRRs (mainly located in the west CONUS) decreases and water yield in the WRR2-8 increases.The increased Q in the wet months tends to intensify the flooding risk, while decreased Q in the major dry months would be likely to aggravate the water shortage conditions.Taking California (mostly in the WRR18) as an example, the monthly Q would decrease by around 5 mm during spring through early summer (the major runoff generation season) due to coupling changes in P and ET.The decrease in flow may cause severe water shortage similar to what happened in 2014-2015in California (Aghakouchak et al., 2014;Mao et al., 2015).Hydrological changes will bring many impacts on water-related economic sectors.For example, droughts would reduce low flows and degrade water quality (high water temperature and nutrient concentrations), thus bringing harmful influence on fisheries (Magoulick and Kobza, 2003;Dolbeth et al., 2008;Gillson et al., 2009), navigation (Theiling et al., 1996), and recreations (Thomas et al., 2013).
The modeling results suggest that GPP over the whole CONUS would increase 106 gC m −2 yr −1 (9 %) in the future.The increase by WRR ranged from 49 to 202 gC m −2 yr −1 or from 5 to 17 % among the 18 WRRs.These findings suggest that carbon stock and vegetation capacity to sequester atmospheric CO 2 for the entire CONUS and each WRR tend to be enhanced under the SRES A2 climate scenario.For the intra-annual GPP changes to climate change, most WRRs show GPP increases, particularly during late spring to summer with higher rates, which implies that the capability of ecosystem to sequestrate carbon in these months will be significantly enhanced in future.By contrast, several WRRs would decrease GPP in several months.For example, during August and September, GPP in WRR17 decreased.The ecosystem sequestration carbon capability would be weakened in these months under the SRES A2 climate scenario.For forests, variations of GPP caused by climate change will be ultimately reflected in timber production, soil carbon storage, and other ecosystem such as dissolved carbon loading in aquatic ecosystems.According to this study, under the SRES A2 climate scenario, the forest biomass and timber production is expected to increase, thus climate change may have implications for timber prices in timberland-dominated regions (Sohngen and Mendelsohn, 1998;Irland et al., 2001;Alig et al., 2004).At the same time, densification of forest lands under a warming climate may provide conditions of increased wildfire potential (Liu et al., 2013c).

Conclusions
We assessed the impacts of future climate change on the hydrological cycle and GPP over the CONUS by linking an ecohydrology model (i.e., WaSSI) with WRF using dynamically downscaled HadCM3 model climate data under the IPCC SRES A2 emission scenario.The current study represents a coupling of bias-corrected, dynamically downscaled climate data with an ecohydrological model to address regional ecosystem issues.The study provides a potential scenario of likely impacts of future climate change on watershed hydrology and productivity across the CONUS, including 82 773 12-digit HUC watersheds.Although only one future climate scenario (the SRES A2 emission scenario) and one GCM (HadCM3 model) was employed here, the methodology applies to other scenarios when more climate change scenarios generated from the WRF are available.
Future climate change will not likely change the spatial patterns of precipitation, temperature, ET, Q, and GPP.However, a large spatial variability in the hydrological and ecosystem productivity responses is expected among the watersheds at both 12-digit and 2-digit HUC scales.The assessment results provide a benchmark of water yield and ecosystem productivity across the whole CONUS, the 18 WRRs and even the 82 773 12-digit HUC watersheds.This type of information will be useful for prioritizing watershed restoration and developing specific measures to mitigate the negative impacts of future climate to sustain the terrestrial ecosystem on different spatial scales (i.e., 12-digit HUC and WRR).

Figure 1 .
Figure 1.Location of the Water Resource Regions (WRRs) over the CONUS (a) with the percentage of each land use/cover type within each WRR.The numeral from 1 to 18 in panel (a) represents the number of WRR.In panel (b), the rectangle size notes the percentage of each land use/cover type within each WRR.Note that the percentages of each land use/cover type were calculated based on the 2006 National Land Cover Database (NLCD) of CONUS.

Figure 2 .
Figure 2. Characteristics of precipitation and temperature during the baseline (1979-2007) and the future periods, and the future changes (future -baseline).

Figure 3 .
Figure 3. Spatial distribution of ET and Q during the baseline and the future periods, and the future changes.

Figure 4 .
Figure 4. Spatial distribution of GPP during the baseline and the future periods, and climate change impacts (future -baseline).

Figure 5 .
Figure 5. Monthly precipitation (a), temperature (b), ET (c), Q (d), and GPP (e) for the whole CONUS during 1979-2007 and 2031-2060 (the top of each panel), and their differences (future -baseline) between the two periods (the bottom of each panel).

Figure 6 .
Figure 6.Number of the WRR within a given interval of change (future minus past) for each month: (a-e) precipitation (P ), temperature (T ), ET, Q, and GPP, respectively.The rectangle size for each month represents the number of the WRR that fall in a given interval value.

Table 1 .
Multi-year mean precipitation, temperature, ET, Q, and GPP averaged over each WRR or the entire CONUS during the baselineand the future period (2031-2060).