Hydrologic Response of Sierra Nevada Mixed-Conifer Headwater Catchments to Vegetation Treatments and Wildfire in a Warming Climate

We used a hydro-ecologic model (RHESSys) constrained by measurements of stream discharge, and spatially distributed snow and soil moisture, to simulate the impacts of operational forest treatments, historical wildfire and climate warming on productive mixed-conifer forests. We compared the response of two headwater catchments at the rain-snow-transition elevation in the wetter central Sierra and more water-limited southern Sierra. The variability of precipitation exerted a greater influence on annual evapotranspiration and runoff than vegetation changes from operational fuels treatment or historical wildfire. The short-term impacts of vegetation changes associated with wildfire, however, did have a greater effect on evapotranspiration and runoff than temperature increases in a warming climate. The average central-Sierra headwater response of evapotranspiration and runoff to fuels treatments (−12%, +12%, respectively) and wildfire (−43%, +46%) were greater than the projected responses to a 4.5°C temperature increase (+2 and −7%). The response in the southern Sierra was limited by lower annual precipitation and showed no response to fuels treatments; but the catchment showed respective changes of −11 and +17% in evapotranspiration and runoff for wildfire, versus +9 and −3% to a 4.5°C temperature increase. These results suggest that in the central Sierra, reductions in vegetation from either fuels treatments or historical wildfire can, temporarily, offset reductions in streamflow from a warming climate. In the southern Sierra, impacts of fuels treatments were small, and only more-extensive vegetation removal as would occur with wildfire, results in significant changes in hydrologic fluxes. Further research is needed to investigate how initial hydrologic changes and climate effects evolve as vegetation adapts and regrows following disturbance.


INTRODUCTION
It has long been known that mountain watersheds such as those on the western slopes of the Sierra Nevada, and the critical water supplies originating in these areas, are sensitive to climate warming (Pupacko, 1993;Jeton et al., 1996). Many of the winter storms that provide the deep seasonal snowpack in the western Sierra occur at temperatures above −3 • C (Bales et al., 2006), making precipitation vulnerable to a transition toward a higher rainfall fraction and reduced snowpack storage in a warmer climate (Knowles, 2002;Miller et al., 2003). The frequency of wildfires is also increasing in western U.S. forests as the temperatures warm, and fuel loads remain high, the legacy of a century-long history of suppressing the frequent low-intensity fires that previously kept vegetation densities low (Westerling et al., 2006). Climate warming is changing the timing, amount and quality of mountain water supplies, as demand pressures grow and as policy requires balancing statewide water supply and demand.
Climate projections for California point to a 3.1-5.0 • C increase in temperatures by the year 2100, with annual precipitation increasing or decreasing as much as 15% (Pierce et al., 2018). Greater changes in Sierra Nevada precipitation are possible by 2050-2100, as indicated under a highemissions scenario (Garfin et al., 2013). While confidence in precipitation projections is low to medium, taken together the projected temperature and precipitation reinforce the need to consider climate scenarios that are hotter and drier than even the most-severe droughts of the past 1100 years (Griffin and Anchukaitis, 2014).
One of the impacts from climate warming is an increase in the fraction of precipitation falling as rain versus snow and effects on runoff will depend on both climate and land use attributes (Bales et al., 2018). However, Berghuijs et al. (2014) suggested that catchments with a higher fraction of snowfall have higher streamflow than would otherwise be expected from precipitation and potential evaporation. A number of studies focusing on watershed response to changes in climate have been completed for the western slope of the Sierra Nevada, which has been identified as more sensitive to changes in temperature than the eastern slope due to the larger area of lower elevation (Pupacko, 1993). For projected temperature increases of 2-5 • C in the American and Merced river basins, and no change in precipitation amount or timing, Dettinger et al. (2004) reported that average-annual runoff generally remained constant, despite changes in the fraction of precipitation falling as rain and earlier snowmelt. These trends that have already been observed in the last half of the twentieth century (Stewart et al., 2004). An analysis of historical streamflow in the Sacramento River basin found that the interannual variability in precipitation explained 95% of differences in annual streamflow volumes while only 3% was explained by temperature (Risbey and Entekhabi, 1996), consistent with the minimum 80% of streamflow explained by precipitation reported by Duell (1994).
Wildfire is one of the key risks to North American ecosystems from climate change (Romero-Lankao et al., 2014). Westerling et al. (2006) discuss the competing influences of climate and forest management on increasing wildfire occurrence across the western U.S., suggesting that although recent climate change was the primary driver in Northern California, fire exclusion is also an important contributing factor in this region. In response to these changing conditions, Millar et al. (2007) encourage a proactive planning approach for forest management. Fuels treatments are an effective forest-management tool for mitigating wildfire risk in Sierra Nevada forests (Stephens, 1998;Collins et al., 2011;Stephens et al., 2013), and include selective thinning and prescribed burning for promoting fire-resilient landscapes (Agee and Skinner, 2005).
Forest-vegetation density and structure impact the interception of precipitation (Storck et al., 2002;Moeser et al., 2015), evapotranspiration amounts (Dore et al., 2010(Dore et al., , 2012Hawthorne et al., 2013), and the surface energy balance for snowmelt (Essery et al., 2008;Ellis et al., 2011;Mahat and Tarboton, 2012;Lundquist et al., 2013). Molotch and Meromy (2014) found elevation, temperature and precipitation were more influential than vegetation, using regression-tree analysis to rank relative physiographic and climatic influence on snow cover for the major Sierra Nevada watersheds. A modeling study using the Distributed Hydrology Soil Vegetation Model (DHSVM; Wigmosta et al., 1994) suggested that both forest cover and temperature increases will have significant, non-linear effects on snowpack and streamflow in the upper Tuolumne (Cristea et al., 2014). The relative effects of temperature and vegetation may then depend on the specific montane elevation range, vegetation type, and annual precipitation received in a watershed, requiring more-localized analyses to determine the dominant influences on evapotranspiration and runoff.
The specific aim of this study was to project the interacting effects of climate warming with forest treatments and disturbance on the annual water balance of productive Sierra Nevada mixedconifer forests in the elevation range that transitions from rainto snow-dominated precipitation (1500-2500 m). We focus on the extent to which reductions in vegetation that are consistent with relatively light thinning prescriptions and historical wildfire versus increasing temperatures will affect the partitioning of precipitation between evapotranspiration versus runoff.

MATERIALS AND METHODS
This study used a hydro-ecologic model, the Regional Hydro-Ecologic Simulation System (RHESSys version 5.14.7; Tague and Band, 2004), to integrate distributed-snow, distributedsoil-moisture and stream-discharge measurements, and to project water-balance response of two Sierra Nevada mixedconifer headwater catchments to temperature and vegetation perturbations. Extensive catchment and model descriptions are available in Saksa et al. (2017), with the most-relevant information provided here. Two types of vegetation change were simulated, a forest-thinning treatment implemented in 2012, and impacts of wildfires modeled with and without the thinning. Projected temperature increases from two climate scenarios were also simulated, Representative Concentration Pathways (RCP) 4.5 and 8.5, at 2050 and 2100. The climate and vegetation scenarios were then simulated together to determine the dominant factors controlling evapotranspiration and runoff, assessed over the range of dry-to-wet precipitation conditions observed during a 4-year period (2010-2013) for which field measurements were carried out.

Study Sites
Two headwater catchments at different latitudes along the western slope of the Sierra Nevada were monitored for climate, stream discharge, distributed snow depth and soil moisture during water years 2010-2013 (Figure 1). Bear Trap Creek (1.4 km 2 , 1560-1826 m elev) is located in the headwaters of the American River basin, in the central Sierra, and Big Sandy Creek (2.2 km 2 , 1776-2475 m elev) is located in the Merced River Basin, in the southern Sierra. Observed discharge was calculated from stream-level data recorded every 15 min, and a stage-discharge relationship developed for each stream. The catchments are dominated by mixed-conifer forests, a forest type covering 13-14% (∼52,500-56,500 km 2 ) of California (Barbour and Minnich, 2000), and have well-drained soils, classified as loamy-sand in Bear Trap and sandy or sandy-loam in the Big Sandy catchment. The headwaters receive a mix of rain and snow precipitation at the elevation of the catchment outlet, but are snow dominated at the upper elevations, suggesting these basins would be sensitive to lower snowfall and higher rainfall from increases in temperature.

Model Scenarios
Four vegetation scenarios were combined with climate projections to determine dominant influences in forest hydrology: no-treatment, thinning-treatment, no-treatment with wildfire, and thinning-treatment with wildfire. Strategically Placed Landscape Treatments (SPLATs; Finney, 2001), a fuels treatment strategy to lower the risk of high-severity fire by treating part of the landscape, were implemented at the fireshed scale during the summer of 2012. As part of SPLAT implementation on a larger fireshed, the mixed-conifer forest was selectively thinned in 95% of the Bear Trap catchment and in 32% of the Big Sandy catchment. It should be noted that these treatments reflect operational decisions of the local forest managers, which are constrained by topography, wildlife habitat, public input, budgets and other factors (North et al., 2015;Lydersen et al., 2019). LiDAR data, described in Kelly and Guo (2015), was used to determine vegetation-community type, and forest-plot measurements were performed before and after SPLAT implementation. Forest plot data were used to impute forest-structure characteristics into the individual vegetation areas (Su et al., 2016), to capture vegetation changes in both horizontal (canopy cover) and vertical (Leaf Area Index) (Saksa et al., 2017).
Understory-vegetation cover was also estimated using a linear equation developed from forest-plot data relating basal area and canopy cover with shrub-cover fraction (Hopkinson and Battles, 2015). The vegetation densities before and after fuels treatments were used to run the Forest Vegetation Simulator (FVS; Dixon, 2002) with the Fire and Fuels Extension (FFE; Reinhardt and Crookston, 2003) and the Fire Area Simulator (FARSITE; Finney, 2004) to project treatment impacts on wildfire severity and vegetation mortality (Fry et al., 2015). As the fire model was calibrated with recent observations, the post-fire vegetation scenarios reflect fire behavior under current climate conditions. FVS was then used to estimate forest-vegetation regrowth for 10 years after the simulated fire events. Following the 10 years of regrowth, overstory canopy cover was transferred directly from FVS, and LAI was calculated from tree lists using allometric estimates.
The model did not consider changes in surface characteristics such as soil hydrophobicity, reduced soil-infiltration capacity, and diminished litter cover that can occur immediately after fire. While these can be important in more semi-arid regions, the impacts are likely to be small in this region where hydraulic conductivity is quite high in the loamy-sand to sandy-soil textures, and infiltration excess runoff production is relatively rare. Wildfire simulations were based on current 95th percentile weather and fuel moisture conditions, as more extreme weather and wildfire events are expected with climate warming, we consider this a conservative scenario.
The climate scenarios were based on changes projected in the minimum and maximum daily temperatures for RCP 4.5, defined as a 4.5 W m −2 increase in radiative forcing relative to pre-industrialization with stabilization by 2100; and for RCP 8.5, which represents an 8.5 W m −2 radiative increase by 2100 that continues to rise (van Vuuren et al., 2011). Mean-annual minimum and maximum temperature anomalies were calculated using the 4-year annual mean of 2010-2013 as a baseline. The 4-year baseline period was + 0.4 • C (−0.2 • C to +0.9 • C) and + 0.1 • C (−0.8 • C to + 0.9 • C) of the long-term climate mean for minimum and maximum temperature in the Sierra climate region, respectively (California Climate Tracker; Abatzoglou et al., 2009). A 4-year trailing ensemble mean was calculated using the Coupled Model Intercomparison Project Phase 5 (CMIP5) output for the Community Climate System Model version 4 (CCSM4) and the Model for Interdisciplinary Research on Climate version 5 (MIROC5) (Figure 2). These models were chosen because they showed low error in bias analyses (Kattsov et al., 2013), and both were archived in the CMIP5 database with the required daily minimum and maximum temperatures. Using the temperature anomalies to produce a uniform offset in minimum and maximum temperatures on our observed-temperature data set aside the need for downscaling climate projections. Vapor pressure and relative humidity were derived using standard air temperature relationships in the model simulations (Tague and Band, 2004), and respond accordingly with increasing temperatures. Impacts on snowpack from vapor pressure and relative humidity are minor, with greater differences during the ablation period in higher elevations and wetter years (Roche et al., 2018a).

Water-Balance Model
RHESSys was used to simulate the hydrologic response to vegetation and climate scenarios. Model calibration was completed for the pre-thinning water years of 2010-2012, during which annual precipitation varied from drier than normal (−39%) to wetter than normal (+60%) conditions in the Sierra Nevada (Saksa et al., 2017). Drainage and subsurface-storage processes were calibrated by comparing simulated and observed streamflow at a daily time step, with precipitation not attributed to evapotranspiration or runoff considered subsurface bypass flow, a comprehensive term to account for all subsurface storage and routing, as in Saksa et al. (2017). Garcia et al. (2013) provides additional details on RHESSys storage and drainage parameters and standard calibration. Monte-Carlo style calibrations were completed by running 5000 sets of random parameters and FIGURE 2 | Climate-scenario temperature anomalies based on the 2010-2013 4-year mean. Each line connects data points for annual-mean temperature anomalies calculated as a 4-year trailing mean from 1950 to 2100. Shaded areas note the range of daily minimum and maximum temperature anomalies from the two climate scenarios used to calculate the mean.
selecting the parameter sets that conformed to a Nash-Sutcliffe and log-transformed Nash-Sutcliffe statistic higher than 0.60, as well as annual and August streamflow rates within 25% of observed values. Six parameter sets met the criteria for the Bear Trap catchment and 17 sets were acceptable for the Big Sandy catchment, providing a range of modeled responses to temperature and vegetation perturbations.
Increasing temperatures at the elevation range of these catchments will have a significant impact on the amount of precipitation received as snow versus rain, and the persistence of snowpack during the winter and into spring. RHESSys calibration was completed using a separate rain and snow precipitation input, and changes to initial snowfall rates were implemented using the linear transition of snowfall temperatures in the model, from −1 to 3 • C. Three inputs contribute to snowmelt rates: temperature, precipitation falling as rain, and radiation. All three components are affected by air temperature, and the relevant toplevel equations are listed below from Tague and Band (2004). Melt attributed to temperature is based on an empirical relationship to sensible and latent heat, and is calculated as: where β MT is the temperature-index melt coefficient, calibrated to 0.0005 for Bear Trap and 0.001 for Big Sandy, T air is the temperature in Celsius and F is the fraction of forest cover. Snowmelt from advection due to rainfall is calculated as: where ρ water is the density of water, TF is net throughfall onto the snowpack, cp water is the heat capacity of water, and λ f is the latent heat of fusion. Lastly, melt due to radiation is calculated as: where K direct and K diffuse are direct and diffuse net shortwave radiation, and L is net longwave radiation. Melt only occurs when the snowpack is ripe, but snow loss also can occur by sublimation from radiation energy input, calculated by adding the latent heat of vaporization to the latent heat of fusion in Eq. 3 (λ f + λ v ). RHESSys was also used to estimate overstory and understory transpiration and evaporation of water intercepted by forest canopy and litter, as well as soil evaporation and snow sublimation. RHESSys simulates vertical infiltration and drainage between surface storage, plant rooting zone, and unsaturated and saturated zones. Lateral redistribution of water follows topography. Previous application of RHESSys in snowdominated mountain environments demonstrates that the model can capture the impact of climate variation on eco-hydrologic processes such as streamflow (Zierl et al., 2007;Garcia et al., 2013), the impact of increases in atmospheric carbon dioxide concentration on conifer net primary productivity and water use efficiency (Vicente-Serrano et al., 2015), and snowpack (Christensen et al., 2008;Godsey et al., 2013;Morán-Tejeda et al., 2014;Bart et al., 2016). Additional details on RHESSys process representation can be found in these studies, plus Tague and Band (2004), and the open-access code maintained online 1 .

Vegetation and Climate Changes
Mean-annual water-balance components of evapotranspiration and runoff were assessed in response to vegetation and climate scenarios over the complete observation period (water years 2010-2013). Mean annual precipitation for this period was 1990 mm in Bear Trap and 1300 mm in Big Sandy. Selective thinning implemented in Bear Trap reduced mean LAI (Canopy Cover) from 9.9 (0.51) to 9.1 (0.49), with reductions from modeled wildfire being 8.8 (0.37) with SPLATs and 7.7 (0.29) without SPLATs. Wildfire in the catchment had a mean flame length of 2.2 m and crowning in 42% of the area without SPLATs, reduced to a mean flame length of 1.2 m and crowning in 19% of the area with SPLATs.
Mean LAI in Big Sandy was 5.0 (0.55 Canopy Cover), and the limited commercial thinning did not reduce the mean catchment LAI substantially (change of < 0.1). A small section of Big Sandy was thinned, with LAI being reduced by as much as 4.0; but the minor amount of area thinned combined with incremental increases in growth elsewhere led to the small changes in basinscale LAI. Wildfire in Big Sandy reduced LAI to 3.8 (0.47 Canopy 1 https://github.com/RHESSys/RHESSys Cover), with thinning prior to wildfire having an insignificant effect. Because of the limited catchment-scale thinning impacts on Big Sandy, results are only reported for the no-treatment and post-fire vegetation change. Wildfire in the catchment had a mean flame length of 1.5 m and crowning in 22% of the area.
Observed mean-daily winter temperatures during the months of heaviest precipitation (Nov-Apr), were 4.3 • C in both Bear Trap and Big Sandy catchments. Projected increases in meanannual temperature with RCP 4.5 were 1.2 • C by 2050 and 1.6 • C by 2100 in Bear Trap, with slightly smaller increases at Big Sandy, 1.0 • C by 2050 and 1.4 • C by 2100. In the RCP 8.5 projections, temperature increases by 2050 and 2100 are 1.8 and 4.7 • C in Bear Trap and 1.6 and 4.4 • C in Big Sandy, respectively.

Water-Balance Simulations
Simulations showed that vegetation changes had much greater effects on runoff and evapotranspiration than did the changes in temperature (Figure 3). Ninety-five percent confidence intervals were calculated for the 6 Bear Trap and 17 Big Sandy model calibrations, with runoff and evapotranspiration responses reported as fractions of precipitation. Confidence intervals were small for both pre-treatment and post-fire scenarios in Big Sandy. Confidence intervals in Bear Trap increased with decreasing vegetation and resulted in higher uncertainty of water-balance response with greater vegetation disturbance. In Bear Trap, the scenario of greatest vegetation change (no treatment, fire) increased the mean runoff fraction by 0.20 (from 0.44 to 0.64) and decreased the mean fraction of evapotranspiration by 0.20 (from 0.47 to 0.27), in response to the 22% LAI decrease (from 9.9 to 7.7) and 42% canopy decrease (0.51 to 0.29). This is equivalent to a drop in ET of about 398 mm yr −1 , from 935 to 537 mm yr −1 (Figure 4). In comparison, the climate scenario of greatest temperature increases (RCP 8.5, 2100) with no change in vegetation resulted in a smaller reduction in runoff, from 0.44 to 0.41, and smaller evapotranspiration increase, from 0.47 to 0.48. Responses of mean runoff and evapotranspiration fractions in 2050 and 2100 to RCP 4.5 temperature increases were limited to less than 0.03, as was the response to RCP 8.5 in 2050.
In the Big Sandy catchment, the simulated response of evapotranspiration and runoff to the 26% LAI decrease (from 5.0 to 3.7) and 14% canopy decrease (0.55-0.47) from modeled wildfire was more limited than in Bear Trap, increasing the runoff fraction by 0.06 (from 0.35 to 0.41) and decreasing the evapotranspiration fraction by 0.06 (from 0.54 to 0.48). This is equivalent to a drop in ET of about 78 mm yr −1 , from 702 to 624 mm yr −1 (Figure 4). A similar response in the reduction of evapotranspiration of 0.05 (from 0.54 to 0.49) was simulated due to the RCP 8.5 temperature increase by 2100, but the response of runoff was not significant, only increasing from 0.35 to 0.36. Precipitation not accounted for in runoff or evapotranspiration is routed to deeper groundwater storage and thus changes in evapotranspiration do not necessarily translate directly into runoff change in the model. Response to all other climate scenarios also resulted in changes of less than 0.03 in evapotranspiration and runoff fractions, similar to Bear Trap, even though Big Sandy has a higher overall evapotranspiration fraction due to the lower precipitation.

Precipitation Variability
Runoff in Bear Trap varied between wet, average, and dry years, and runoff response also varied across vegetation scenarios (fire and treatment). The difference in runoff between wet and dry years exceeded 750 mm while the runoff difference between no treatment/fire runoff and runoff from the most-substantial vegetation disturbance was limited to less than 500 mm (Figure 4). Trends of runoff and evapotranspiration response with increasing temperatures were the same in all climate scenarios except for RCP 8.5 in 2100, where evapotranspiration increases were greater in the wet years (<50 mm). In Big Sandy, precipitation variability between wet and dry years had a greater effect on evapotranspiration and runoff than the reductions in vegetation from wildfire (Figure 4). The differences in annual precipitation resulted in evapotranspiration differences near 400 mm without treatment or fire, and runoff differences close to 700 mm with post-fire vegetation. Water-balance response to reductions in LAI from wildfire were smaller, with evapotranspiration decreases and runoff increases of <200 mm. Elevated temperatures in 2050 and 2100 increased ET during mean to high precipitation years, offsetting some of the ET decline due to fire.

Hydrologic Timing and Storage
Precipitation falling as rain or snow, along with the accumulation and melt of the seasonal snowpack, determines the timing of soil infiltration, runoff and availability of water in the soil for use by vegetation. Projected temperature increases affected both precipitation phase and melt rate, and changes in vegetation density impacted snowmelt by modifying the surface-energy balance. Simulations of temperature and vegetation impacts on hydrologic storage and timing were assessed for 2010, a mean precipitation year (Figure 5).
In Bear Trap, post-fire vegetation losses advanced the snowpack melt-out date by about 3 weeks, while temperature increases by 2100 in RCP 4.5 advanced the melt-out date by about 4 weeks (Figures 5A,B). Increases in temperature by 2100 in RCP 8.5 showed no persistent snowpack, with the fraction of precipitation falling as snow decreasing from 0.40 to 0.10 ( Table 1). The reduced vegetation after fire increased soil-water storage during the dry season from minimal storage to within 150 mm of saturated winter conditions (∼350 mm). Increased temperatures resulted in an earlier start of the soil-water storage recession from wet winter to dry summer conditions by approximately 1.5 weeks, associated with the reduced snowpack. Evapotranspiration with increased temperatures becomes more limited in the early summer because of the earlier soil drying. The reduced post-fire vegetation resulted in evapotranspiration reductions during all seasons and was not limited by soilwater storage. Runoff timing was accelerated with both scenarios of increased temperatures and reduced vegetation, with peak runoff about 8 weeks earlier. Post-fire vegetation resulted in 67% higher annual runoff, with increased peak (+54%) and FIGURE 4 | Evapotranspiration (upper 6 panels) and runoff (lower 6 panels) for forest-treatment and disturbance scenarios during dry (2012), mean (all years), and wet (2011) precipitation conditions for current and projected temperatures Vertical bars indicate the 95% confidence interval based on the multiple parameter sets for each catchment. Left column of panels is Bear Trap and right column of panels is Big Sandy. Big Sandy mean basin LAI following treatments were not different than having no treatment and are not shown.
Frontiers in Forests and Global Change | www.frontiersin.org  frequency (5-14 events) of high-flow events, where increasedtemperature scenarios resulted in lower peak flow (−35%) at similar frequency. In Big Sandy, post-fire vegetation resulted in the loss of a persistent winter snowpack, similar to temperature increases by 2100 with RCP 8.5 (Figures 5F,G). Increases in temperature by 2100 with RCP 4.5 and 8.5 reduced the snowfall fraction from 0.60 to 0.53 and 0.29 of precipitation, respectively. Soilwater storage from infiltration increased earlier in the winter and became more saturated with temperature increases (300 mm) than in the control (230 mm), but the dry season recession curve also started about 4 weeks earlier than in the control scenario. Evapotranspiration response to the vegetation reductions was muted, but evapotranspiration was reduced in all vegetation and temperature simulations from the earlier drawdown in soil-water storage. Peak runoff timing was shifted about 12 weeks, from early June to early March, in the post-fire and RCP 8.5 scenarios, with higher peak flow (+29%) in the post-fire vegetation simulation than with projected increases in temperature (−9% to +4%). The earlier peak runoff in the RCP 4.5 scenario was only shifted about 8 weeks, to mid-April.

DISCUSSION
These are highly productive ecosystems compared to other temperate conifer forests (Millar, 1996), where vegetation change from wildfire, and to a lesser extent forest management, impacted evapotranspiration and streamflow more than projected temperature increases from climate warming. On the other hand, impacts from inter-annual climate variation between relatively wet and dry years were greater than the impacts of vegetation change or climate warming. The hydrologic impacts of wildfire were greater than those associated with fuel treatments and the magnitude of vegetation change impacts were greater in the wetter central-Sierra study site. These results are based on two relatively small Sierra Nevada watersheds that enable the incorporation of substantial observation data for constraining a physically based model, but that also limit the ability to fully capture the wide range of geoclimatic variation within the Sierra. Nonetheless, our findings demonstrate how climate and vegetation change effects can interact to produce significant, short term hydrologic changes, but also show why it can be challenging to meaningfully generalize about these impacts in the context of climate and site differences.
Estimated hydrologic responses to vegetation change in this study were relatively high, particularly for the wetter Bear Trap catchment. Boisramé et al. (2019) found streamflow increases of about 20-40 mm as managed fires were allowed to return to the Illilouette Creek Basin, which is lower than the ∼75 mm change in the climatically similar Big Sandy Creek catchment. The difference can likely be attributed to a greater vegetation change from a reduction in forest cover in the Big Sandy headwater simulated wildfire, a single event that burned the entire catchment during 95th percentile weather conditions (e.g., high temperature, low humidity) in which high-severity wildfires are more likely to occur. Vegetation change in the Illilouette Creek catchment is a result of the ongoing managed wildfires that burned portions of the larger Illilouette Creek catchment in various weather conditions over 40 years, and incorporated conversion of forested areas to shrubs and meadows. Bart et al. (2016) showed the potential for water balance changes in the southern Sierra region in excess of 100 mm with reduced vegetation, and changes in vegetation type from forest to shrubs. The greater 400 mm yr −1 decline in ET with high-intensity wildfire in the high-precipitation, mixed-conifer region of the central-Sierra American River basin was approximately 70 mm higher than the findings of Roche et al. (2018b) and Roche et al. (2020). For example, the nearby 2014 King Fire resulted in a 330 mm (±80 mm) ET reduction in areas burned by highseverity fire, averaged over 4 years after the fire event. Saksa et al. (2020) also demonstrated the higher response of the water balance to horizontal forest structure changes in RHESSys, as such the 42% reduction in canopy cover contributed substantially to the change in ET.

Water-Balance Simulations
Ficklin and Barnhart (2014) suggest using multiple parameter sets and General Circulation Model outputs, as significant differences in hydrologic projections will occur from uncertainty in model parameterization and climate scenario. In this study, we used the ensemble means from two climate models and two emissions scenarios, with 6 (Bear Trap, American River Basin) and 17 (Big Sandy, Merced River Basin) parameter sets to incorporate some accounting of uncertainty. Additional uncertainties can originate from climate downscaling, which was not done in this study, and model structure (Wilby and Harris, 2006). We simulated a uniform increase in the minimum and maximum daily temperatures, using projected temperature anomalies for 2050 and 2100. Studies often use a single mean-temperature adjustment to project effects of climate warming on Sierra Nevada watersheds (e.g., Young et al., 2009;Meyers et al., 2010), but minimum and maximum temperatures may have periods of non-linear increases, which can modify the diurnal temperature range (Easterling, 1997;Vose et al., 2005).
In the RCP 4.5 projections, minimum and maximum temperatures increased an average of 0.017 and 0.020 • C yr −1 in the American and 0.016 and 0.019 • C yr −1 in the Merced between 2013 and 2100, respectively. RCP 8.5 projections of minimum and maximum temperatures resulted in mean increases of 0.044 and 0.055 • C yr −1 in the American and 0.045 and 0.053 • C yr −1 in the Merced between 2013 and 2100, respectively. The rate differences of minimum and maximum temperature will impact snow accumulation and melt estimates, along with ecological estimates in RHESSys, such as the estimation of vapor-pressure deficit and the Jarvis-based calculations of stomatal conductance (Jarvis, 1976) used to estimate transpiration in RHESSys, which incorporates functions of mean and minimum daily temperatures to limit maximum conductance.
The simulation results of temperature perturbations generally agree with the findings of Dettinger et al. (2004) that without statistically significant changes in precipitation, annual volumes of streamflow will generally remain steady in the American and Merced River basin areas with increasing temperatures. Dettinger et al. (2004) did find a general trend of +8% precipitation per year with climate warming, but the effect was small compared to the interannual variation, which was also found to have a large effect on runoff in this study (Figure 4). Previous work in the Merced River basin (Christensen et al., 2008) also showed the sensitivity of transpiration to precipitation at an elevation range similar to Big Sandy, with their model extending into higher elevations where transpiration became increasingly sensitive to changes in temperature. The sensitivity to transpiration at higher elevations may also be why Bales et al. (2018) estimated a greater increase in evapotranspiration (+12-15%) with a +1 • C temperature increase in the Kings River Basin during the recent California drought, leading to a −5% change in runoff over the basin, calculated as precipitation minus evapotranspiration. Null et al. (2010) used the Water Evaluation and Planning (WEAP21) waterbalance model that includes supply, demand, and policy scenario capabilities to project progressively decreasing mean annual flows in all major Sierra Nevada basins with increasing temperature, including a respective −5.6 and −6.3% for the Merced and American with +4 • C. Temperature increases in this study were similar for the RCP 8.5 scenario in 2100, resulting in a similar response for projected mean runoff in Bear Trap (American) of −5.3%, but a lesser runoff response in Big Sandy (Merced) of −1.2%. Changes in Water Use Efficiency and evapotranspiration with increasing atmospheric CO 2 concentrations were not considered in this study (De Kauwe et al., 2013).

Precipitation Variability
The strong response of evapotranspiration to vegetation density in Bear Trap, and to annual precipitation in Big Sandy, suggests that within a Budyko (1974) framework of competing water and energy limitations of transpiration, the Big Sandy region tends to be more water limited (Figure 6). The difference in energy and FIGURE 6 | Annual evapotranspiration as a function of dry to wet precipitation conditions in Bear Trap and Big Sandy. Forest vegetation was constant (pre-treatment), with fitted lines highlighting the different trends between the two catchments. water limitations will affect the magnitude of the water-balance response to changes in vegetation and temperature, consistent with Zhang et al. (2001), who showed the potential for increased response of evapotranspiration with reduced forest cover in regions with higher precipitation. In both sites, yearly runoff was influenced by interannual precipitation variability more than temperature, similar to results from previous work in the western Sierra (Duell, 1994;Risbey and Entekhabi, 1996).
The individual years of 2011 and 2012 were selected from the 4 years of simulation (2010)(2011)(2012)(2013) to provide a spectrum of response to climate and vegetation perturbations during a wet and dry year, respectively. Antecedent-moisture conditions can modify watershed response to disturbance, so the progression of dry to wet years may be important. Precipitation in 2010 was close to the long-term mean for both regions, followed by the wet year of 2011, and dry years of 2012 and 2013. Shallow subsurface water-storage capacity and the rates of soil drainage versus evapotranspiration will impact the magnitude of summer-baseflow response in low-precipitation years and where temperature increases lead to earlier snowmelt (Jefferson et al., 2008;Tague and Grant, 2009;Huntington and Niswonger, 2012). The uncertainty associated with the calibrated model parameter sets of subsurface flow increased in Bear Trap Creek with simulated reductions in vegetation (Figures 3, 4). Improved characterization of subsurface properties in Sierra watersheds are needed to enhance our understanding and predictive capability of runoff response to climate (Shaw et al., 2014), which could be used to further constrain the simulation uncertainty of vegetation scenarios in this study. The total subsurface storage capacity, and plant-available water storage (Garcia and Tague, 2015;Tague et al., 2019), could limit the evapotranspiration increases seen in warmer Sierra Nevada watersheds that have been used as a proxy for temperature increases with climate change (Goulden and Bales, 2014).

Hydrologic Timing and Storage
Both reduced vegetation and increased temperatures resulted in more energy being absorbed by the snowpack, with persistent snow cover eliminated by 2100 in RCP 8.5 (Figures 5B,H). In contrast to results from Berghuijs et al. (2014) that indicated streamflow will generally decrease with decreasing snowfall, these results are consistent with studies such as Tague and Peng (2013) that argue that snowmelt inputs can in some cases support higher rates of evapotranpiration than winter precipitation. In this study, the increased winter (November-April) snowmelt led to higher water storage in the soil and increased streamflow, while evapotranspiration generally remained unchanged. The earlier melt reduced soil-water storage, evapotranspiration and streamflow in May through July in both Bear Trap and Big Sandy (Figure 5). These results agree with Tague and Peng (2013), showing that evapotranspiration response to temperature increases will depend on timing of snowmelt. Forrest et al. (2018) also project higher winter flows by 2050 for watersheds supplying California hydropower facilities, a significant source of renewable-power-generation -but noted that due to increased spillage events, higher flows did not necessarily yield additional power generation.
The response of snowmelt to vegetation loss from wildfire in Big Sandy was much greater than in Bear Trap. LAI drives radiation attenuation (Varhola and Coops, 2013), and in RHESSys, radiation transferred through the canopy to the snowpack follows a Beer's Law type of exponential curve. In Bear Trap, LAI was reduced from 9.9 to 7.7 in the post-fire scenario, is within the saturated range of radiation absorption by the canopy, and is similar to previous reported LAI values for ponderosa pine forests in the area (Goldstein and Hultman, 2000;Gersonde et al., 2004;Campbell et al., 2009). In Big Sandy, the LAI was reduced from 5.0 to 3.7, which is similar to previous LiDAR LAI estimates for the Sierra National Forest (Zhao et al., 2012;Tang et al., 2014). The lower LAI in Big Sandy is within the exponential increase of radiation with changes in vegetation, resulting in snowpack patterns more similar to RCP 8.5 in 2100. The combination of reduced vegetation and elevated snowmelt led to the highest peak runoff in the post-fire scenario, compared to the control or climate projections, but will be further influenced by forest gap size and slope orientation (Ellis et al., 2013).
Changes in the timing of snowmelt and runoff also have implications for fish habitat and other aquatic species that rely on summer baseflow. Using the end of the 2010 water year (Sep 30) as an indicator of low flow in this study showed in baseflow reduced more following wildfire (−19%) than temperature increases (−4 to −8%) in Bear Trap Creek. In Big Sandy Creek, however, baseflow was reduced more following temperature increases (−29 to −33%) than wildfire (−25%). Godsey et al. (2013) note than low flows in this region depend on both Snow Water Equivalent level of both current and previous years, and impacts are influenced by subsurface storage capacity and underlying bedrock. Meyers et al. (2010) further extend their temperature warming assessment to winter and spring flooding in the region, which may negatively impact both brook and rainbow trout, but the greater proportion of winter flooding may impact brook trout more severely.

CONCLUSION
The water-balance responses to simulations of temperature and vegetation perturbations in productive mixed-conifer forests with wildfire disturbances showed that vegetation changes from operational fuel treatments and historical wildfire exerted a greater influence on annual evapotranspiration and runoff than did projected temperature increases in a warming climate. However, inter-annual variation in precipitation had a greater influence on runoff than did effects due to either fuels treatments or wildfire. Hydrologic impacts associated with historic wildfire were generally greater than those associated with operational fuel treatments. In the wetter central Sierra, headwater evapotranspiration decreased, and runoff increased 40-50% for a simulated wildfire event. In the more waterlimited southern Sierra, the headwaters response was constrained to a respective −11 and +17% change in evapotranspiration and runoff following a simulated wildfire event. In contrast, evapotranspiration increases and runoff decreases to a 4.4 • C temperature increase were less than 10% of current values for all headwaters.
Climate warming will eliminate the persistent seasonal snowpack at these elevations (1500-2500 m) by 2100 in the RCP 8.5 projections, becoming rain dominated as the amount of precipitation falling as snow was reduced from the current 40-60% down to 10-29%. The early snowmelt has cascading effects on the rate and timing of soil-water storage, and thus on evapotranspiration and runoff. Increases in temperatures resulted in peak streamflow occurring up to 12 weeks earlier, but post-fire vegetation conditions also increased peak runoff, especially in the American River headwaters of Bear Trap Creek. These results suggest that in the central Sierra, reductions in vegetation from either light thinning and fuels treatment or historical wildfire generate increases in flow, however, only wildfire or more-extensive fuels treatments provide significant increases in the southern Sierra. In both cases, these reductions are premised on vegetation changes that can be sustained. Further work to examine the dynamics of vegetation regrowth following fuels treatments, wildfire, and vegetation adaptation to a warming climate are needed to determine longer-term impacts of vegetation change on runoff and evapotranspiration.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
PS contributed to the study design, data processing, modeling, and drafted the manuscript. MC and RB contributed to the study design, modeling approach, and manuscript revisions. CT contributed to the modeling approach and execution, and manuscript revisions. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We would like to thank California Department of Water Resources, United States Forest Service, Sierra Nevada Research Institute, and Sierra Nevada Adaptive Management Project personnel for their support and assistance. We thank Patrick Womble and Sarah Martin for their significant work in field data collection at UC Merced, and to Norm Miller at UC Berkeley for assistance with the climate datasets. We thank John Battles and Danny Fry at UC Berkeley for the forest vegetation and wildfire modeling components incorporated in this study. We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups for producing and making available their model output. For CMIP the U.S. Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals.