The magnitude and spatial patterns of historical and future hydrologic change in California ’ s watersheds

.


INTRODUCTION
Projections of future climates from General Circulation Models (GCMs) diverge more widely for precipitation than for temperature (IPCC 2007).Ecologically, precipitation is highly relev www.esajournals.org 1 February 2015 v Volume 6(2) v Article 24 vant, as moisture availability is a critical determinant of terrestrial ecosystem type (Milly et al. 2005, Trnka et al. 2012).There is a need, therefore, to address hydrologic dynamics under climate change as part of climate impact assessments (Girvetz et al. 2009, Marcarelli et al. 2010).This has led to attempts to couple hydrologic models with climate projections, with the goal of simulating future hydroclimatic dynamics.Some hydrologic models use daily or hourly climate data to drive watershed-scale hydrologic projections (e.g., Leavesley et al. 1992, Wood et al. 2002, Hay et al. 2011).Such models are informative for flood and drought forecasting.However, characterizing the landscape for ecological applications requires fine spatial scale, which can result in computational limitations over large regions.
Simpler, process-based hydrologic models have a number of advantages for climate change studies.First, the most common and accessible resolution for relevant data are monthly observations, which are generated on a year-by-year basis by the GCM groups who develop forecasts using IPCC emissions scenarios (IPCC 2007).Hydrologic models that can use these monthly climate projections as data inputs can render temporally simpler, but informative assessments of watershed sensitivities to changing climate (Flint and Flint 2012a, b).Second, for effective comparison of ecosystem sensitivities to climate change, hydrologic models should address entire landscapes, rather than single watersheds more typical of temporally intensive hydrologic model applications.This broader scope provides insight into multi-watershed or landscape-scale ecosystem dynamics, as well as being able to calculate hydrologic response with the spatial resolution to capture habitat-relevant processes (Band et al. 1996).Process-based, monthly time step models can achieve these goals with considerably lower computational requirements than finer temporal scale models.Third, if a process-based model is calibrated to historical conditions it can be useful in climate change studies because it can potentially better represent local and regional future responses to change than statistical or rule-based models (Cuddington et al. 2013).Relevant to all hydrological models, natural resource managers often require finer spatial scale data than are generated by GCMs, therefore requiring downscaling of GCM outputs for use in resource planning (Littell et al. 2012).
We used the Basin Characterization Model (BCM; Flint andFlint 2007, Flint et al. 2013) to examine change in hydrologic conditions between historical and current time periods : 1951-1980 (historical) and 1981-2010 (current), and from current time to two futures (2071-2100) that represent drier and slightly wetter future climate conditions.The BCM is a process-based model with several modules (Fig. 1A) that balances the hydrologic cycle for each grid cell in the area modeled.Modules in the BCM simulate four parts of the hydrologic cycle: snowpack, climatic water deficit, groundwater recharge (as net infiltration below the root zone), and surface runoff.
The BCM differs from many watershed models by its ability to spatially characterize net infiltration below the root zone corresponding to variation in shallow bedrock, which may become groundwater recharge (Flint and Flint 2012a, b).The BCM uses mapped bedrock type and estimates of bedrock permeability to partition available water into spatially distributed recharge and runoff, thus providing an additional layer of local constraint to estimates of recharge that are typically determined as everything that is not runoff, and only represented at a basin scale.The model's inclusion of both runoff and recharge across large areas for multiple projections is a capacity we sought for this study, because the representation of recharge is important under future climates, as water availability becomes more uncertain.
In addition, running the BCM at fine spatial resolution can quantify hydrologic dynamics at scales relevant to vegetation and ecological processes.This ability permits the use of the BCM outputs for a large range of natural resource management issues, and its outputs have already been used in climate-related studies of scale sensitivity of species distribution models, forest die-off, and fire severity in California and the western United States (for example Franklin et al. 2012, Millar et al. 2012, van Mantgem et al. 2013).Development of spatially downscaled climatic inputs to 270-m grid cells supported the fine-scale application of the BCM in this study (Flint and Flint 2012a, Thorne et al. 2012, Flint et al. 2013).The model is by design sufficiently parsimonious to enable application over large areas, only including parameterization needed to provide meaningful output.We used monthly time steps because of the future downscaled projections available for analysis, and because monthly time steps enable the evaluation of long time periods necessary for analyses of future climates over large areas.This study assesses the hydroclimatic dynamics of climate change across 5135 watersheds comprising California, a Mediterranean climate region where mountain snowpack has declined under historical warming (Mote et al. 2006), melting of annual snowpack is occurring earlier (Knowles et al. 2006), and where decline in available water from precipitation can have considerable impacts at ecological, ecosystem and economic levels (e.g., Hayhoe et al. 2004, Westerling et al. 2006, Tanaka et al. 2006, and Howitt et al. 2014).
This study used the BCM output maps of historical and two future hydroclimatic condi-tions, representing wetter and drier conditions, to compare the hydrological change among 5135 watersheds in the California hydrologic region (a region encompassing California and all the watersheds that flow into it, Fig. 1B, here called California), in a comparative assessment of watershed-scale response to changing climate.Focusing on four BCM outputs-snowpack, climatic water deficit (CWD), recharge, and runoff-we asked: 1. What is the magnitude of historical and projected future change in the hydrologic cycle? 2. Are watersheds that show the greatest historical change the same as those predicted to have the greatest change under future conditions? 3. To what degree are the same watersheds most sensitive to hydrologic stresses under both wetter and drier futures?

METHODS
Model outputs for snowpack, CWD, recharge, and runoff at a 270-m grid scale were spatially averaged by watershed.For each time period (i.e., past, present, and future) and variable analyzed, we used 30-year monthly and yearly means and their standard deviations (SD), calculated for each grid cell.The 30-year time slices are based on the recommendation to use 1951-1980 as a baseline, during which time the North American climate changed relatively little (Hansen et al. 2012).The current and future 30year times used are extensions of that reasoning.We used the standard deviation (SD) of the current time period (month or year) as a measure of the return interval against which to evaluate the change in hydrologic conditions for the past and projected future climatic conditions.For comparison of the historical period to current, we spatially downscaled the historical monthly and water-year climate maps of California which are portrayed by 800-m grids from the Parameter-elevation Relationships on Independent Slopes Model (PRISM; Daly et al. 2008).We summarized change by watershed and ecoregion in units of millimeters (mm) of water gained or lost for each variable, and as measured by the number of SDs from the current 30-year mean for each variable.We assume that for some ecological questions, change as measured in variability could provide more insight into impacts than changes indicated by the difference in mean conditions between periods (Landres et al. 1999).We also quantified change from current to future projections in 2071-2100, for two GCMs, PCM and GFDL using the A2 emissions scenario, that represent respectively a wetter and a drier future in California.We developed watershed-specific values for 5135 watersheds (USGS HUC 12 watersheds, 2011 version; http://nhd.usgs.gov/index.html)that comprise the region, including CWD, April 1st snowpack, runoff, and recharge.

Input data
Historical climate maps (minimum and maximum air temperature [Tmin, Tmax] and precipitation [PPT]) were derived from the empiricallybased PRISM monthly precipitation and air temperature database and maps that are available at 800-m spatial resolution (Daly et al. 2008).
We used a spatially downscaled version of the PRISM data (270-m grids) that was developed using a gradient-inverse-distance squared approach (Flint and Flint 2012a) that integrates location and elevation with monthly climate values.The future climate projections used were statistically downscaled to 12-km using constructed analogues (Hidalgo et al. 2008, Flint andFlint 2012a), and spatially downscaled to 800-m to perform a bias correction to adjust projections to match the mean and standard deviation of 1950-2000 climate (Flint andFlint 2012a, Thorne et al. 2012).This correction provides the ability to compare historical conditions and future projections without introducing differences in means or variation.Future projections were then spatially downscaled to 270-m for model application.
We selected the A2 emissions scenarios from two GCMs: GFDL (the National Oceanic and Atmospheric Administration's Geophysical Fluid Dynamics Laboratory CM2.1 model; Delworth et al. 2006, Stouffer et al. 2006) and PCM (the Parallel Climate Model developed by National Center for Atmospheric Research and the U.S. Department of Energy; Meehl et al. 2003); based on analyses done by Cayan et al. (2008Cayan et al. ( , 2009)), who used several criteria for selection of models to downscale for California.The GCMs selected were required to produce a realistic simulation of aspects of California's recent historical climate, particularly the distribution of monthly temperatures and the strong seasonal cycle of precipitation.They were required to contain realistic representations of regional features, such as the spatial structure and variability of precipitation.In addition, the models selected were to have differing levels of sensitivity to greenhouse gas forcing.Emissions scenarios are based on assumptions that range from highly mitigated global emissions to business-as-usual, and current CO 2 measurements reflect levels that exceed those used in any of these models (Carbon Dioxide Information Analysis Center; www.cdiac.ornl.gov).Therefore, we only apply projections using the medium to high CO 2 emissions reflected in the A2 emissions scenario (Special Report on Emissions Scenarios, IPCC 2007).
These scenarios were compared by Krawchuk and Moritz (2012) to an ensemble of 16 CMIP3 global climate models (Meehl et al. 2007) who found that by 2010-2039, GFDL in , GFDL in California is in the mid-range for temperature increases and low range for precipitation; while the PCM is in the low range for temperature and precipitation.By 2070-2099, the GFDL is still in the mid-range among GCMs for temperature but among the lowest of the 16 GCMs for precipitation, thus supporting its use in representing a warmer, drier future.The PCM falls in the low range for temperature and low-mid range for precipitation; annual precipitation shows little change, but is higher than the GFDL.Although not ranking high among precipitation projections within the 16 GCM ensemble, in comparison to GFDL the PCM projects a wetter future.

Running the BCM
Each module of the BCM represents a process (Fig. 1A), meaning the model as a whole is mechanistic, rather than statistical.The assumption of this approach is that once the model is calibrated, the relationships between deterministic factors will hold steady as climate changes, and permit more accurate estimates of local environmental conditions, particularly when comparing the relative hydrologic dynamics between basins, where the influence of site factors may cause variation in hydrologic outcomes.Map-based variables required to run the BCM in addition to temperature and precipitation are the bedrock type with estimates of bedrock permeability, depth and hydraulic properties of soil, and the slope, aspect, and elevation from a digital elevation model.These were assembled (Thorne et al. 2012), and the model was calibrated by using historical stream gauge data from 159 relatively unimpaired watersheds around California (Flint et al. 2013).
The calculation sequence in the BCM starts with potential evapotranspiration (PET; Flint and Childs 1991) calculated on the basis of solar radiation (slope, aspect, elevation, and topographic shading; Flint and Childs 1987), and is a modified version of the Priestley-Taylor PET algorithm (Priestley and Taylor 1972).Precipitation (PPT) and temperature are used to calculate the proportion of PPT that falls as snow, and the accumulation, sublimation, and melt of snowpack (Fig. 1A).Available water calculated from PPT and PET (the energy loading), and the associated snow processes, interact with soil values in each grid cell to calculate soil moisture, which permits assessment of the hydrologic dynamics of the soil.From these, actual evapotranspiration (AET) and climatic water deficit (CWD; CWD ¼ PET À AET) can be calculated (Fig. 1A).When soil moisture capacity is exceeded by available water, the BCM allocates the surplus to runoff and infiltration of water below the root zone (here called recharge).These two are linked in the BCM, and permit a calibration of calculated recharge or runoff by varying the permeability of the underlying bedrock to partition excess water between these two fates.To calculate basin discharge, the recharge and runoff grid cell values are summed for all locations within a watershed above where a stream gauge has recorded streamflow, and goodness-of-fit is achieved to measured data by post-processing equations derived to represent subsurface processes and the gains and losses in streams, including seepage, base flow, and deep groundwater recharge.This study's model calibrations were developed from 159 historical streamgage measurements across California (Thorne et al. 2012, Flint et al. 2013).Averaged goodness-of-fit statistics for those basins resulted in r 2 of 0.73 and Nash-Sutcliffe Efficiency statistic of 0.67, indicating relatively good matches of modeled to measured streamflow, including basins that varied geologically and hydrologically across the entire region with different amounts of impairments due to urbanization, diversions, impoundment, or agriculture (Flint et al. 2013).
We ran the BCM using monthly climate input data, summed the monthly outputs to calculate water year values (for California's Mediterranean climate, water year was defined as October-September), and summarized 30-year means for climatic and hydrologic output variables.The times used were 1951-1980;1981-2010;2010-2039;2040-2069;and 2070and -2099and . The 1951and -1980 time period was used as baseline against which to measure historical change, and the current time period  was the baseline to assess future changes.For the current time, we calculated the average and SD values for each variable in each grid cell.For all other times, we calculated the mean values over 30 years.All outputs were summarized by watershed.These results were used to assess the change per hydrologic variable by watershed, as measured in units of current time standard deviations as follows: Change ðin current SD unitsÞ ¼ T 2 30-year mean À T 1 30-year mean current 30-year SD where T 1 is 1951-1980 for historical analyses and 1981-2010 for future analyses, and T 2 is 1981-2010 for historical analyses and 2070-2099 for future analyses.The standard deviation analysis permits a perspective how of much change has occurred, or will occur, relative to current variability.
Results were summarized by watershed (n ¼ 5135) and for 10 ecoregions in California (Fig. 1B; modified from Hickman 1993).For the snow module, we present snow accumulated by April 1st, the date also used by water managers in California to assess yearly water supply held in snow, as a measure of anticipated springtime runoff.This measure combines the snow accumulated during California's precipitation season with the effects of warming air temperatures that influence the seasonality of snowmelt (Stewart et al. 2004).To reduce the volume of results presented, only future change to 2070-2099 is presented, and this paper focuses on change in the four hydrologic output variables: annual CWD, recharge, runoff, and snowpack on April 1st.Results from additional variables are presented in Appendix: Tables A1 and A2, and Appendix: Figs.A1-A5.

BCM output analysis
The four variables in millimeter (mm) units, were either accumulated over the water year, or accumulated by April 1st, and the yearly values were used to produce the 30-year means.For each variable, we then sorted the 5135 watersheds by amount of change between time periods (increase in CWD, and decrease in recharge, runoff, and snowpack) and identified the 514 watersheds (10%) with the greatest change towards drier and more ecologically stressful conditions between time periods, here called the ''area pattern analysis'' to evaluate where within California change occurred.We assigned a value of 1 to each selected watershed and added selected watersheds from the four analyses to see which watersheds had high levels of change for more than one hydrologic cycle variable.
We also selected watersheds where the change in each variable ranked in the top 20% of all change (magnitude of change, not a set number of watersheds as above).This cutoff identifies a subset of the watersheds selected by the area pattern analysis.The results are presented together in maps that show both the pattern for the 514 watersheds in the area pattern analysis, and the 20% of watersheds with the greatest amount of change.The selection of 10% by area and 20% by value as cutoffs is arbitrary, but permits identification of enough watersheds to visualize ecoregional patterns of change.
Finally, we normalized the change in each variable for the historical and future analyses.For variables that showed both positive and negative change, we used the absolute value of change.This produced a 0-1 value for change per variable per watershed.The values of the four variables could then be averaged per watershed, permitting a rank-ordering of watersheds with an index of overall hydrologic change.This index was used to compare how similar or different are the spatial patterns of the most affected watersheds under the two future scenarios (spatial congruence; Seo et al. 2009).The index was also used to assess spatial congruence of historically high-ranking watersheds with the distribution of future highly ranked watersheds.We used spatial cutoffs of 514 (10%) and 1028 (20%) watersheds with the highest index values from the historical and future analyses to examine the spatial patterns of change among all time periods.Results from these analyses are shown to illustrate relationships among watersheds at ecoregional scales, and to highlight regions that will emerge as newly-stressed under future projections.

BCM model outputs
Current yearly CWD watershed values range from 2 to 1519 mm (Fig. 2A), with standard deviations (SD) between 0 and 182 mm (Fig. 2B).Since 1951Since -1980, annual , annual CWD in California has increased by 17 mm, with the greatest ecoregional increase (36 mm) east of the Sierra Nevada (Fig. 2A; all ecoregional results in Appendix: Table A1).Change from historical to current time (positive or negative) is not .1 SD of current variability for any location, and mostly ,0.5 SD (Fig. 2B).Under the GFDL projection, CWD increases throughout California (140 mm), with particularly large increases in the Sierra Nevada, East of Sierra Nevada, and the Modoc Plateau.Many watersheds in the Sierra Nevada increase by between 1.5 and .2.0 SD of current CWD.Under the PCM projection, CWD increases in all but a few watersheds, although by a lesser amount than under the GFDL projection (Fig. 2C-E; 61 mm average for all California).These increases range from 0.5 to 1.5 SD of current values for the Sierra Nevada, Central western CA, and Northwestern CA ecoregions, with smaller increases elsewhere (Fig. 2F-H).Under both GFDL and PCM projections, the greatest increase in CWD is in the southern Sierra Nevada.
Current April 1st snowpack ranges from 0 to 1089 mm (0-511 mm SD; Fig. 3A, B).With the Fig. 2. Climatic Water Deficit (CWD) averaged over the current 30 years and summarized by watershed (A), and its standard deviation (B).CWD represents unmet water demand by plants, which can be a measure of stress to plants.Historical and future changes in CWD are shown, as differences from the current 30-year mean value (mm) and by the number of current 30-year standard deviations (C-E and F-H, respectively).Historical change in CWD is typically less than 0.5 SD, while in projected futures, CWD for GFDL ranges from 0.6 to 2 SD, and PCM from 0 to 1.5 SD.The legends for C-E and for F-H are distributed across each set of panels.v www.esajournals.orgexception of very high elevation increases in snowpack, confirmed by observation (Mote 2006, Andrews 2012), spring snowpack has declined historically by 10 mm for California, with the greatest declines in the Cascade Ranges (À33 mm) and the Sierra Nevada (À29 mm), and in many watersheds this change is .1 SD (Fig. 3C, F).Future snowpack declines under both GFDL (À45 mm) and PCM projections (À35 mm) with the greatest declines for both in the Sierra Nevada, Cascade Ranges, and Modoc Plateau (Fig. 3D, E).
Current runoff values range from 0 to 2800 mm, with standard deviations of 0-837 mm (Fig. 4A, B).From the 1951-1980 time period, changes in runoff are patchy, with almost no change across the whole state, but a decrease of 22 mm in the Northwestern CA ecoregion, and an increase of 12 mm in the Central western CA ecoregion (Fig. 4C).For nearly all locations the changes represent less than 0.5SD from current values (Fig. 4F).Under the GFDL projection, California Fig. 3. April 1st snowpack averaged over the current 30 years and shown by watershed (A), and its standard deviation (B).Historical and future changes in April 1st snowpack are shown, as differences from the current 30year mean value and by the number of current 30-year standard deviations (C-E, F-H, respectively).Historical change ranges from 0.5 to .2SD, while future projections show decreases for GFDL ranging from 0.5 to 2 SD, and decreases under PCM from 0 to 1 SD.The legends for C-E and F-H are distributed across each set of panels.v www.esajournals.orgrunoff decreases À36 mm, with the biggest decreases in the Cascade Ranges (À123 mm), Northwestern CA (À95 mm), and the Sierra Nevada (À82 mm) (Fig. 4D).Under the PCM projection, California runoff increases by 26 mm, with the greatest increases in the Sierra Nevada (68 mm) and Northwestern CA (61 mm) (Fig. 4E).Runoff under the GFDL projection is mostly within 1 SD of current values (Fig. 4G).There are some locations in very dry parts of the state that show larger SD change values, but occur where changes of a few mm are large relative to the very low amount of runoff to start with.Runoff predominantly increases around the state under the PCM projection and many parts of the Great Valley (91 watersheds out of 555) and desert Current California recharge values range from 0 to 1634 mm, with standard deviation values ranging from 0 to 399 mm, and the California average is 127 mm (Fig. 5A, B).Change from historical recharge is mostly very low, and below 0.5 SD of current time (Fig. 5C, F).Under the GFDL projection, California recharge decreases by a mean of 30 mm, with the greatest decrease in the Northwestern CA ecoregion (À103 mm) (Fig. 5D).Under the PCM projection, California recharge increases by 2 mm.The greatest increases are in the Central Valley (þ14 mm) and Central Coast (þ13 mm) while the Northwestern ecoregion sees a decrease of À19 mm (Fig. 5E).As with runoff, very dry regions show great change in standard deviation units, but

BCM output analysis
According to the area pattern analysis (selection of the top 10% of watersheds, n ¼ 514), the largest historical change for recharge and runoff is concentrated in the Northwestern CA ecoregion, while changes in CWD and April 1st snowpack were concentrated in the Sierra Nevada and Cascade Ranges (Fig. 6A-D).Historical changes in CWD ranged from 40 to 87 mm, with associated decreases between 42 and 275 mm for April 1st snowpack, 14 and 125 mm for recharge, and 12 and 300 mm for runoff (see Appendix: Table A2 for full area pattern analysis results).Forty-four watersheds are in the top 20% change by value with historical increases in CWD between 63 and 87 mm; 12 watersheds had decreases in snowpack between 180 and 275 mm, and fewer numbers of watersheds for historical change in recharge and runoff.When overlaying the four variables for the area pattern analysis, the Sierra Nevada and Northwestern CA ecoregions emerge as the areas with the greatest historical change (Fig. 6E).
Under the area pattern analysis for GFDL, the top 514 watersheds experience increases in CWD between 192 and 336 mm, declines in April 1st snowpack from 178 to 798 mm, declines in v www.esajournals.orgrecharge from 105 to 328 mm, and declines in runoff from 117 to 389 mm.The Sierra Nevada Mountains are the location of the most change for all variables, while the Northwest CA ecoregion also shows high levels of decline for recharge and runoff (Fig. 7A-D).These numbers are amplified in the 20% most impacted watersheds by value: CWD in 72 watersheds increases above 269 mm, April 1st snowpack decreases more than 638 mm in 20 watersheds, recharge decreases more than 256 mm (13 watersheds) and runoff decreases more than 302 mm (14 watersheds).When the area pattern analysis results are overlaid, the Sierra Nevada Mountains emerge as the ecoregion with the most watersheds that will be negatively impacted in two or more of the hydrologic categories (Fig. 7E).
The area pattern analysis for the PCM projection shows a similar intensification of CWD and decrease in snowpack compared to change under the GFDL projection, primarily in the Sierra Nevada and Northwest ecoregions (Fig. 8A-D).However, runoff and recharge actually increase slightly in most mountainous regions.Decreases in runoff identified in the area pattern analysis are in the Modoc Plateau (Northeast CA), and in Northwest CA for recharge (Fig. 8A, B).Overall, fewer watersheds are projected to be negatively impacted by multiple variables under the PCM projection than the GFDL projection (Figs.7B v www.esajournals.organd 8B): for CWD, change in the top 20% of watersheds is projected to be 34% less than under GFDL; while for April 1st snowpack depletion it is 24% more; and for recharge it is 62% less.
The index of normalized values of hydrologic change identifies an 82% overlap among the top 514 watersheds with the most change under GFDL and PCM projections.These watersheds are located primarily in the Sierra Nevada ecoregion, (263 of 716 watersheds; Fig. 9A).Some differentiation appeared using 20% of all watersheds (n ¼ 1028), with 78% overlap for most affected watersheds in the Sierra Nevada (n ¼ 378) and Northwestern CA ecoregions (n ¼ 213 of 744 watersheds; Fig. 9B).Watersheds in the Central western CA (n ¼ 60 of 477) and Southwestern CA (n ¼ 43 of 361) ecoregions respond to increasing water availability under the wet PCM projection, but not the dry GFDL projection (Fig. 9B), while an additional 151 Northwestern CA watersheds show high levels of water depletion under the GFDL projection but not the PCM projection.
When comparing historical and future hydrologic change by index for the top 514 watersheds, 46% had already undergone change historically.When considering the top 1028 watersheds with the greatest future change, the proportion already showing change rises to 56%.These watersheds are concentrated along the highest elevations in the Sierra Nevada ecoregion, as well as in the Northwestern CA ecoregion (Fig. 10A, B).Many of the watersheds in this index not among the top 514 historically, but active under both GFDL and PCM projections, occur in three regions: the highest elevation watersheds of the Sierra Nevada ecoregion (in the south), in lower elevations in the northern Sierra Nevada ecoregion, and west of the Mount Shasta region (Fig. 10B).Areas that are among the 514 most active historically that are not projected to be so under the future include portions of the Southwestern CA ecoregion (San Bernardino Mountains and Big Bear Lake; Fig. 10A), and portions of the desert lands east of the Sierra Nevada ecoregion.By contrast, under the 1028 watershed analysis, lower watersheds on the western side of the Sierra Nevada ecoregion and mountain ranges in the Southwestern CA and Central western CA ecoregions begin to show high levels of hydrologic change under one or both of the future projections (Fig. 10B).

DISCUSSION
The use of the BCM, a process-based model, to simulate hydrologic dynamics from climate variables allowed for comparisons of hydrologic change among watersheds under projected future climates.The 5135 watersheds showed considerable variability related to their site characteristics.We identified watersheds that have already experienced and/or are projected to undergo large hydrologic change, for both individual parts of the hydrologic cycle, and as measured by an index that integrates the component parts.Historically, the Sierra Nevada ecoregion has already shown changes in snowpack and CWD, while recharge and runoff changes have been concentrated in the Northwestern CA ecoregion.These trends continue and are amplified under the GCM future Fig. 9. Combined rank-ordered change in combined hydrological outputs for watersheds by 2099.These images portray spatial congruence between the BCM model outputs for GFDL and PCM using the BCM normalized projections of change for the four variables.When 10% of all watersheds are considered (A), relatively few watersheds that show the most response are seen (in red or light blue).However, when 20% of all watersheds are considered (B), differentiation is more apparent, with many watersheds' hydrologic activity decreasing under the drier GFDL model in the Northwestern ecoregion (in red) and increasing in the central and southwestern ecoregions under the wetter PCM model (in light blue).
v www.esajournals.orgprojections we examined, with additional watersheds in varying locations entering active change in their water cycle depending on projected future conditions and the assessment criteria used.This discussion first addresses the integrated hydrologic variables, then the individual variables analyzed, followed by a few comments on model limitations and next steps.

Normalized change
Normalizing the change in each BCM variable allowed watersheds to be summarized and rankordered, and also permitted the use of cutoffs, as measured by the amount of change or by a selection of some proportion of all watersheds, to assess change through time.The normalized values could then be combined into a single index value, which permits an integrated ranking of watersheds.It also permitted the tracking of when a watershed's hydrologic cycle will show high response to climate change, and when it may cease to do so, because there is too little water left (Figs. 9 and 10).Spatial and temporal differentiation of maximum watershed response to climate change was particularly evident when considering the top 20% of watersheds showing change (Fig. 10B).This area-based approach to identifying sensitive watershed regions on a landscape permits a view of the spatial differences across large regions that arise when considering climate projections that diverge in terms of precipitation.
For example, watersheds in pale green in the Southwestern CA ecoregion, and dry watersheds east of the Sierra Nevada ecoregion, have been active historically, but do not rank highly in the future (Fig. 10B).This indicates where water availability becomes so low as to not permit much future change.Many watersheds in the coastal ranges show change under either the wetter or drier futures (Fig. 9B) but have not had change historically (Fig. 10B).These include watersheds that show high levels of response under both futures, and pale blue watersheds that indicate an increase in degree of change (increasing runoff ) under the PCM projection.It Fig. 10.Combined rank-ordered watersheds for historical and future projections using the BCM normalized projections of change for four variables: CWD, snowpack, recharge and runoff.Watersheds that show the most change historically and under both futures are clustered predominantly in the Sierra Nevada Mountains when 10% of all watersheds are considered (A).Many watersheds in the Northwestern ecoregion also appear under a 20% analysis (B).Watersheds in bright red show where relatively little change has occurred historically, but are predicted to experience high levels of change under both futures, while pale tan (drier) and blue (wetter) show the watersheds that will show more change under one, but not both futures.
v www.esajournals.org is interesting to note that the red watersheds include some of the highest elevation watersheds in the Sierra Nevada ecoregion, owing to large reductions in snowpack by the end of the century, but also some lower, foothill watersheds, where changes reflect shifts in the snow/rain transition zone, as well as a large number of locations around the Mt.Shasta region (Fig. 10A, B).
This view of when and where watersheds will show maximum hydrological response to climate change can be helpful in watershed management, and the rank-ordered approach allows water managers to visualize the potential levels of impact.For example, highly productive watersheds used for water delivery can be evaluated to assess when to expect their maximum hydrologic response to drying climate conditions, and their subsequent decline in water availability.Comparing between wetter and drier future projections can identify the watersheds expected to have high levels of change under both projections (Fig. 9A, B), which can help managers in assessing the level of risk to watersheds for which they are responsible.The converse can also be measured, to identify watersheds that have relatively low projected levels of change, although this paper focuses on watersheds at risk, rather than those that are more resilient.The measure of resilience would be to identify watersheds above some minimum level of precipitation that are also in the lowest levels of change using the BCM and normalized index approach.
Ecoregions can also be used to rank the entire study area and understand the region-wide hydrologic responses to climate change.For example, mountain watersheds in the Southwestern CA ecoregion, such as those in the northern San Bernardino Mountains and in the Big Bear Lake region, were identified during the historical change area pattern analysis (Fig. 10A), and cause the Southwestern CA ecoregion to have the fourth highest historical rate of change in CWD among the 10 ecoregions in the study area.However, in the future projections, the Southwestern CA ecoregion drops to seventh in terms of future change in CWD.Climatic water deficit in the Southwestern CA ecoregion continues to increase under future projections, but other ecoregions experience significantly higher rates of CWD increase.

Individual variable dynamics
The historical warming in California occurred in all ecoregions, and is about one-third to onefifth of projected warming under the future projections analyzed (Appendix: Tables A1 and  A2).Overall, precipitation has increased by 8 mm relative to the historical value, although the wettest ecoregion, Northwestern CA, decreased by 33.4 mm, and with a maximum decrease in one watershed of 529 mm.The change in precipitation is higher in the projected futures, increasing statewide by 49 mm under the PCM projection and decreasing by 95 mm under the GFDL projection (Appendix: Table A1).The BCM outputs integrate these contrasting trends to provide site-mediated simulations of the hydrologic cycle, represented by the variation among watersheds in any ecoregion, where aspect, topographic shading, soil properties, or geologic type might have a greater influence on moderating the loss of water in a watershed.The BCM outputs of historical runoff show a very slight increase across all of California, with a decrease under the GFDL projection, and increasing more rapidly than historical under the PCM projection.The BCM outputs for historical recharge decrease slightly across the study region.Under the GFDL projection, the decrease is amplified, while under PCM it changes to a slight increase.The trends in other BCM variables show no change in direction, although in many cases there is a strong amplification of historical rates.
The uncertainty of future regional precipitation drives much of the variation found in the BCM outputs.However, projected increases in air temperature interact with precipitation, particularly for driving PET, as well as by affecting the snow accumulation and melt.PET predominantly showed moderate increases over historical time (,80 mm in most locations), but it rises sharply in future climate models (Fig. 8; Appendix).Particularly under the GFDL projection, increases in PET impact ecosystem and ecological processes by increasing CWD and lowering water available for recharge and runoff.What effects this level of increased metabolic demand will have on the many plant species in the region is undetermined, but water deficit stress well beyond current 30-year variability is one way to identify places where we might expect physiognomic shifts in dominant vegetation types (Lenihan et al. 2008), in part due to increased post-fire mortality among stressed forest trees (van Mantgem et al. 2013).These dynamics may also correspond with shifts in the ranges of plant species (Loarie et al. 2008, Thorne et al. 2008, Loarie et al. 2009, Ackerly et al. 2010).
One of the largest changes in hydrologic processes resulted from increased air temperature that influences springtime snowmelt.The historical extent and thickness of April 1st snowpack diminished in all but a few highelevation locations, which has large implications for the water supply of California, where a slow springtime snowmelt helps sustain both human populations and natural resources through the dry summer months.The fine spatial-scale application of the BCM allowed detection of locations of historical April 1st snowpack increase that generally have not been captured by more coarse-scale climate model applications (e.g., 12-km grid cell studies; Knowles and Cayan 2002, Snyder et al. 2004, Maurer and Duffy 2005), but the phenomenon was captured in a 2-km grid cell model (Howat and Tulaczyk 2005).The BCM model outputs are also corroborated by empirical studies of historical snowpack dynamics (Regonda et al. 2004, Stewart et al. 2004, Mote 2006).Disagreement with interpretation of empirically measured trends in historical snowfall levels by Christy (2012) can potentially be resolved if one considers that it is change in temperature that primarily affects change in depth of spring snowpack.In the BCM outputs, historical change in snowpack influenced runoff and recharge differently across the state.For example, a decrease in snowpack in the Cascade Ranges (À33 mm) led is accompanied by a 7 mm increase in runoff, and a 2 mm decrease in recharge; while an increase of 6 mm of snowpack East of the Sierra Nevada occurred as well as an 8 mm increase in runoff and a 6 mm increase in recharge (Appendix: Table A1).
Modeled change in projected runoff corresponds to the uncertain direction of change in precipitation, and is mediated by the proportion of precipitation that falls as snow and the timing of snow melt.Runoff is projected to continue at similar or increasing levels under the PCM projection, especially in the mountainous regions, but to decrease under the drier GFDL projection.However, projected change in runoff is mostly under 0.5 SD for both the wetter and drier GCMs.The exception is in desert regions, where increases in runoff (particularly for the PCM projection) are projected to be far above current variability, because of the current very low levels of runoff, and the actual magnitude of change projected is very small.While warming is generally expected to cause earlier spring runoff, and to cause more extreme precipitation events (Flint and Flint 2012b) that are not represented in this study because of the 30-year annual means chosen for analysis, our results suggest that for many places, the amount of future annual runoff is similar to current amounts.This is remarkable because projected increases in PET reduce excess water available for runoff.However, our results do not moderate extreme precipitation events that may arise under both futures during which precipitation overwhelms soil water capacity for short durations of high precipitation, yielding greater annual runoff numbers than would otherwise be expected.
Similar to runoff, recharge historically decreased in Northwestern CA because of decreases in precipitation, and increased slightly elsewhere (Appendix: Table A1).Future changes in recharge are also similar to runoff but with less variability, and generally correspond to the wet and dry projections of precipitation by the two GCMs.However, review of changes in the ratio of recharge to runoff is informative.For some basins, recharge becomes more dominant than runoff under the drier GFDL projection.Recharge dominates over runoff when excess water does not overcome soil moisture storage capacity, and bedrock can accommodate the recharge.This generally occurs in locations with deep soils or permeable bedrock (Flint et al. 2013).As future statewide snow cover diminishes under increasing temperature, exposed areas potentially permit recharge to function throughout the precipitation season, rather than only after snowmelt.However, where soils are shallow, the snowmelt processes result in increased runoff for both GCMs.When snowpack melts in the spring, the water quickly overcomes the soil moisture storage and becomes runoff.Our projections indicate that if there is excess water afforded by increases in precipitation, it will become runoff rather than recharge by over an order of magnitude for ecoregions dominated by v www.esajournals.orglow permeability bedrock and shallow soils, such as the Sierra Nevada ecoregion.However, in locations such as the Modoc Plateau with permeable volcanic geology or the Great Valley with deep soils, recharge may prevail over runoff for both wet and dry scenarios (Figs. 4 and 5).
Climate has also driven historical increases in CWD, which proved to be most reliant on air temperature.This is observable where moderate historical increases in precipitation were either held in soils or became recharge or runoff early in the season (due to earlier melt of snowpack), and therefore resulted in increased annual CWD.Increased CWD occurred throughout the state, except for parts of the San Francisco North Bay, Central western CA, and Northwestern CA regions (which experienced less of an increase in air temperature) and the northern Great Valley (where the soils are thick enough to store excess water from precipitation).However, CWD becomes much more pronounced under future projections, with the eastern and northern parts of the state showing the greatest change under both projections (East of Sierra Nevada, Sierra Nevada, Modoc Plateau and Cascade Ranges ecoregions; Appendix: Table A1).This trend is due to increases in PET resulting from changes in air temperature, snowpack, and the limitation of the soils to hold additional water where, or if, precipitation increases.Thus projected increases in CWD, which is potentially a more direct predictor of species distributions (Stephenson 1998), agricultural demand, and other processes that rely on seasonal soil moisture, can be used to identify areas of much higher future stress.
These results suggest that as historical warming continues some watersheds that until now have been hydrologically stable will exhibit significant changes in the volume and timing of their hydrologic cycles.For example, mountainous regions that have been somewhat buffered from hydrologic impacts may experience more pronounced impacts related to water shortages and CWD.These impacts may include tree mortality due to drought, increased pathogen outbreaks, and increased fire risk (e.g., Breshears et al. 2005, Westerling et al. 2006, de la Giroday et al. 2012), all of which may lead to significant changes in timberland species composition and structure, including conversion to other vegetation types.Projected overall drying and lessening of snowpack depth raises interesting management challenges.In relatively arid regions like parts of California, the need to retain as much water as possible in watersheds for infiltration will likely lead to increased focus on ecosystem adaptation management, particularly for forests, riparian zones, and meadows and wetlands to capture water through natural processes (Millar et al. 2012, Capon et al. 2013), and to more intensive management of human-controlled systems such as underground storage and the timing and volume water transfers (Hanak and Lund 2012).
The future conditions of the major agricultural regions will depend on whether the future plays out more along the PCM or the GFDL projections.Under the PCM, CWD increases in most agricultural areas, but the increase is typically less than half a standard deviation of current variability.Climatic water deficit increases are much greater under the GFDL projections, suggesting farmers will need more adaptive measures to address increasing soil aridity (Mehta et al. 2013).Regardless of hydrologic conditions, increased temperatures will impact winter chilling hours necessary for bud set and fruit production, which is already observed (Baldocchi andWong 2008, Luedeling et al. 2009), as well as earlier warm temperatures that could result in mismatches in flowering and pollination (Hegland et al. 2009).For natural ecosystems, the Sierra Nevada appear to be under the greatest increase in stress from CWD under both projections, and the three western CA ecoregions experience varying levels of CWD increase, which are lower under the PCM projections while the GFDL projections produce similar levels of CWD as in the Sierra Nevada.These increases are likely to have both physiological impacts to plant species, and to change the background level of fire return in these ecosystems.

Model limitations and next steps
There are trade-offs that go with the representation of fine-scale processes associated with soils and potential evapotranspiration and the calculation of the water balance over long time periods and large spatial scales.For example, snowpack dynamics typically occur as a series of accumulation events, and a variable rate of melt.However, we only portray the estimated net snowpack accumulated at a single point in time (April 1st); albeit one used by resource managers to estimate the springtime reserve of moisture held by snow that is used for water delivery management during California's dry season.We used snowpack projections on this date to assess the relative changes that climate change will impose.As we are successful in matching the timing of the snowmelt peaks and can calibrate the basin outflow to match the measured outflow (Flint et al. 2013), we believe that use of a monthly snow model, summed for an annual accumulation date is informative.Tests of the BCM historical record have been done to assess how well simulated snowpack matches snow covered area, snow water equivalent at snow courses, and persistence of snowpack at mapped locations of glaciers (Flint andFlint 2007, Curtis et al. 2014).Additionally, exercises to illustrate a match to the monthly runoff peak for Sierra Nevada basins could add confidence to the simulation of snowpack and the projections of April 1st snowpack displayed here.However, it should be noted that BCM outputs are generally in agreement with other studies of future snowpack (Miller et al. 2003, Cayan et al. 2008).There is a limitation to the monthly model because we are currently using only one set of snow accumulation and melt coefficients for the entire state, and the northern Sierra Nevada snow dynamics differ from those in the southern high elevation basins.This could be addressed in the future by calibrating separate models for all the basins, to then develop different coefficients by region.
Another potential shortcoming to this modeling approach is that it does not account for hydrologic connectivity that may redistribute soil moisture and thereby influence patterns of evapotranspiration and recharge.While we acknowledge that this factor is at work, we feel the redistribution of soil moisture driven by these processes typically occurs on a scale smaller than 270 m, so this has not been considered a limitation to the BCM.Further studies on groundwater flow would help clarify this issue.
The use of the BCM outputs for fine-scale representation of climatic and hydrologic trends can be similar to the way in which GCMs are used, in that the model outputs provide a platform for inter-comparison of regions even if individual GCMs used may be more or less accurate when compared to ground-based measurements.Challenges for modelers of these physical and biophysical processes include determining how to incorporate a fine scale level of detail, developing methods for additional model calibration and validation, and developing means of projecting the BCM across larger areas.
Because of the modular nature of the BCM, it is possible to make two types of improvements.First, any particular module's calculations may be updated and improved (Fig. 1).For example, if improvement to the PET equation (Priestley and Taylor 1972) used in PET calculations (Sheffield et al. 2012) or identification of varying PET for different vegetation types could be calculated, these could be applied using an existing vegetation map to render more accuracy in plant-driven parts of the model.Second, input data maps may be updated and improved.Efforts to improve the BCM outputs in our study area are currently focused on improvement to the snow-driven module.Spatially distributed calibration of the snow accumulation and snowmelt calculations to better match measured snow covered area and snow-water equivalent from snow sensors and snow courses, and better define the soil zone exploited by plants for actual evapotranspiration could render more accurate outputs for this module.Despite these simplifications and the need to continue refining individual modules, the comparative approach described here is informative to regional water planning and management.Indeed, developing BCM models for the historical record and future projections for consistent representation across very large areas would be beneficial, and is being developed for the western USA.

Fig. 1 .
Fig. 1.Schematic showing the processes modeled in the Basin Characterization Model (A).The upper left is the snow module, the right-hand shows the steps in evaporative demand.The large box is where climate values interact with site characteristics, which outputs to the evapotranspiration module, and to runoff and recharge.Basin discharge and groundwater recharge represent post-model summaries that can be tallied for areas of different extents, such as watersheds or ecoregions.(B) The Study Area includes all watersheds that drain into California.Ecoregions discussed in the text are shown in blue, and the outlines of individual basins in grey.

Fig. 4 .
Fig. 4. Runoff averaged over the current 30 years and shown by watershed (A), and its standard deviation (B).Historical and future changes in runoff are shown, as differences from the current 30-year means and by the number of current 30-year standard deviations (C-E, F-H, respectively).Historical change is typically less than 0.5 SD, although a few basins in the Great Valley emerge with relatively high levels of change.Under future projections, decreases in GFDL generated runoff typically declines by ,1 SD, although some basins in the desert regions show an increase with higher levels of change.Under PCM much of the deserts and Great Valley ecoregions show increases ranging from 0.5 to .2SD.The legends for C-E and F-H are distributed across each set of panels.

Fig. 5 .
Fig. 5. Recharge averaged over the current 30 years and shown by watershed (A), and its standard deviation (B).Historical and future change in recharge are shown, as differences from the current 30-year means and by the number of current 30-year standard deviations (C-E, F-H, respectively).Historical change is typically less than 0.5 SD, although a few basins in the Great Valley emerge with higher levels of change.Under future projections, decreases in GFDL generated recharge declines are typically ,1.5 SD, although some basins in the desert regions show an increase ranging from 0.5 to 2 SD.Under PCM much of the deserts and Great Valley ecoregions show increases ranging from 0.5 to .2SD.The legends for C-E and F-H are distributed across each set of panels.

Fig. 6 .
Fig. 6.Historical change in BCM variables between 1951-1980 and 1981-2010.The watersheds in orange are the 514 (10%) of 5135 watersheds with the highest amount of change towards drier conditions (decreasing snowpack, recharge and runoff, and increasing CWD) over the historical time (A-D).Watersheds shown in red represent the watersheds whose changes were in the top 20% of change for each variable.The sum of the 514 watersheds that have had the most change from each of the four categories (E) is shown, to identify basins which were most impacted (most exposed) for two or more hydrologic parameters.

Fig. 7 .
Fig. 7. Future change in BCM outputs from 1981-2010 and 2071-2100.The watersheds in orange are the 514 (10%) of 5135 watersheds with the highest amount of change towards drier conditions (decreasing snowpack, recharge and runoff, and increasing CWD) projected under GFDL A2 (A-D).Watersheds shown in red represent the watersheds whose changes were in the top 20% of change for each variable.The sum of the 514 watersheds projected to have the most change from each of the four categories (E) is shown, to identify basins that are in the highest category (most exposed) for two or more hydrologic parameters.

Fig. 8 .
Fig. 8. Future change in BCM outputs from 1981-2010 and 2071-2100.The watersheds in orange are the 514 (10%) of 5135 watersheds with the highest amount of change towards drier conditions (decreasing snowpack, recharge and runoff, and increasing CWD) projected under PCM A2 (A-D).Watersheds shown in red represent the watersheds whose changes were in the top 20% of change for each variable.The sum of the 514 watersheds projected to have the most change in each of the four categories (E) is shown, to identify basins that are in the highest category (most exposed) for two or more hydrologic parameters.This image differs from the GFDL (Fig.7) image in that the northeastern region (Modoc Plateau) is identified for a decrease in runoff.This is due to a declining snowpack, which opens highly porous soils for ground water recharge during the wet season.

Fig. A1 .
Fig. A1.Minimum annual temperature (Tmin) for current time (1981-2010) and the variability over 30 years.Values are shown on a per-watershed basis.There are 5135 represented in this and subsequent images.Values represent the mean (A) or standard deviation (B) for each watershed.Historical and future change in Tmin, relative to the current 30-year annual mean or current 30-year annual standard deviation.The top row (C-E) represents the change in value per watershed, the bottom row represents the same change, as in the number of current time standard deviation (F-H).In these images, historical change is rarely over two SD, but in almost all watersheds the future change is over 2 SD of current time.The legends for C-E and F-H are distributed across each set of panels.

Fig. A2 .
Fig. A2.Maximum annual temperature (Tmax) for the current 30 years (A) and the variability over the current 30 years, portrayed in standard deviation units (B).Values are shown on a per-watershed basis.Values represent the mean (left) or standard deviation (right) for each watershed.Historical and future change in Tmax, relative to the current 30-year means or current 30-year standard deviations.The top row (C-E) represents the change in value per watershed, the bottom row (F-H) represents the same change, as in the number of current time standard deviation.In these images, historical change is rarely over 1 SD, but in all watersheds the future change is over 2 SD of current time.The legends for C-E and F-H are distributed across each set of panels.

Fig
Fig. A3.Precipitation (PPT) averaged over the current 30 years (A) and the variability over the current 30 years, portrayed in standard deviation units (B), shown on a per-watershed basis.Historical and future change in PPT.The top row (C-E) represents the change in value per watershed, the bottom row (F-H) represents the same change, as measured in the number of current time standard deviations.In these images, historical change is never over 0.5 SD, and future change is almost always below 1SD of current time.The legends for C-E are distributed across each set of panels.

Fig. A4 .
Fig. A4.Potential evapotranspiration (PET) averaged over the current 30 years (A) and the variability of PET (B), shown on a per-watershed basis.(B) Historical and future change in PET, relative to the current 30-year means or current 30-year standard deviations.The top row (C-E) represents the change in value per watershed, the bottom row (F-H) represents the same change, measured in the number of current time standard deviations.In these images, historical change ranges up to 2 SD of current time.In the future projections, GFDL is almost always above 2 SD of current time, while PCM shows variation between 1 and .2SD.The legends for C-E and F-H are distributed across each set of panels.

Fig. A5 .
Fig. A5.Actual evapotranspiration (AET) averaged over the current 30 years (A) and the variability of AET (B), shown on a per-watershed basis.Actual evapotranspiration represents the potential water demand, as limited by available water, that is conducted by plants.Historical and future change in AET, relative to the current 30-year means or current 30-year standard deviations.The top row (C-E) represents the change in value per watershed, the bottom row (F-H) represents the same change, measured in the number of current time standard deviations.In these images, historical change is typically less than 0.5 SD, while in projected futures, both GFDL and PCM are typically below 1 SD of current time.The legends for C-E and F-H are distributed across each set of panels.