In ecoregions across western USA streamflow increases during post-wildfire recovery

Continued growth of the human population on Earth will increase pressure on already stressed terrestrial water resources required for drinking water, agriculture, and industry. This stress demands improved understanding of critical controls on water resource availability, particularly in water-limited regions. Mechanistic predictions of future water resource availability are needed because non-stationary conditions exist in the form of changing climatic conditions, land management paradigms, and ecological disturbance regimes. While historically ecological disturbances have been small and could be neglected relative to climatic effects, evidence is accumulating that ecological disturbances, particularly wildfire, can increase regional water availability. However, wildfire hydrologic impacts are typically estimated locally and at small spatial scales, via disparate measurement methods and analysis techniques, and outside the context of climate change projections. Consequently, the relative importance of climate change driven versus wildfire driven impacts on streamflow remains unknown across the western USA. Here we show that considering wildfire in modeling streamflow significantly improves model predictions. Mixed effects modeling attributed 2%−14% of long-term annual streamflow to wildfire effects. The importance of this wildfire-linked streamflow relative to predicted climate change-induced streamflow reductions ranged from 20%−370% of the streamflow decrease predicted to occur by 2050. The rate of post-wildfire vegetation recovery and the proportion of watershed area burned controlled the wildfire effect. Our results demonstrate that in large areas of the western USA affected by wildfire, regional predictions of future water availability are subject to greater structural uncertainty than previously thought. These results suggest that future streamflows may be underestimated in areas affected by increased prevalence of hydrologically relevant ecological disturbances such as wildfire.


Introduction
Since the 1950s accelerated human activities have altered key Earth cycles and systems in ways that have profound, often indirect, future effects on renewable water resources (Karl and Trenberth 2003) and have made the assumption of hydrologic cycle stationarity untenable (Milly et al 2008). Acknowledging that our climate is non-stationary, future water resource availability relative to the present is typically projected on the basis of future greenhouse gas emissions scenarios and their consequent impacts on precipitation and evaporative demand, as determined by hydrometeorologic models (Milly et al 2005). Implicit in these highly impactful projections of climate change effects on water resources is an assumption that anthropogenic climate change is essentially the only component of global-scale human-induced changes relevant to future water resource availability. Indeed the accuracy of projections of future water resource availability is assumed to be high, bolstered by assessment of an ensemble of climate models (Milly et al 2005). These projections are in turn considered and heeded by western USA water managers as they formulate long-term plans to cope with large reductions in the availability of renewable freshwater resources. Rising temperatures, earlier snowmelt, and changing land management paradigms-including reduced grazing rates and a transition from suppressing to managing wildfire-have contributed to increased wildfire frequency in western North America. Following wildfire hydrologic fluxes (blue arrows) are modified to reflect reduced transpiration and enhanced groundwater recharge.
However, in addition to effecting modifications to climate forcings used in conventional physically based hydrologic predictions, climate change also effects ecological disturbances; these climatic controls on ecological disturbance may in turn be amplified by other components of global change (Dale et al 2001, Scheffer et al 2001. Widespread upland global change-induced ecological disturbances include woody encroachment (Van Auken 2000, Wine et al 2015), tree mortality due to pest-pathogen complexes (Kurz et al 2008), and wildfire (Westerling et al 2006). Assuming that these ecological disturbances can substantially influence Earth's future water balance, they represent structural uncertainties in the current approach to projecting future water resource availability. Though there has been progress in including wildfire in dynamic vegetation models that simulate hydrological processes (Krinner et al 2005), the long-term impact of wildfires on ungauged basins throughout the western USA remains unknown. To date there is evidence of increased water yields due to pest-pathogen complexes (Bearup et al 2014) and decreased water yields due to woody encroachment (Huxman et al 2005). However, the extent of these ecological disturbances in the western USA is not yet sufficiently widespread to estimate their future effects on regional water availability from retrospective environmental data.
In contrast, evidence is mounting that increases in western USA wildfire activity since the 1980s (Westerling et al 2006) significantly enhance regional water availability (Kinoshita and Hogue 2015, Bart 2016, Wine andCadol 2016). This increase in wildfire activity in the western USA is occurring in response to changing climatic conditions and land management paradigms (Parisien et al 2016) (figure 1). Climatic changes increasing the frequency of large wildfires include warmer temperatures and earlier snowmelt (Westerling et al 2006). Land management paradigms contributing to an increase in large wildfire frequency include reduced grazing rates (with resultant increases in fine fuel loads) and an outstanding fire deficit due to fire suppression since the 1940s (Marlon et al 2012).
Hydrologic consequences of wildfires have been widely investigated, but typically at a local or regional scale, in investigations that implement disparate methodologies. While enhanced infiltration-excess overland flow and associated flooding are well-known consequences of wildfires (Wine and Cadol 2016), these fires also reduce transpiration (Dore et al 2010), increase soil water storage (Cardenas and Kanarek 2014), and enhance groundwater recharge (Ebel 2013) and accompanying baseflow (Stoof et al 2012, Kinoshita and Hogue 2015, Bart and Tague 2017. Enhanced post-wildfire infiltration excess overland flow-a consequence of fire impacts on soil hydraulic properties-is often erroneously attributed to soil hydrophobicity, when in fact surface cover, which prevents soil sealing, exerts a more important control on enhanced post-wildfire infiltration excess overland flow (Larsen et al 2009). In contrast, the principle mechanism underlying post-wildfire increases in groundwater recharge immediately following wildfires is reduced (or eliminated) transpiration, with reductions in precipitation interception ranking as secondary (Ebel 2013). Following the initial disturbance, the formerly forested system may be colonized by herbaceous plants having different functional traits than coniferous forest stands. During this vegetation recovery phase-conventionally referred to as secondary succession-shallower rooting depth and shorter growing season duration (in that order of importance) of herbaceous plants relative to evergreen trees may increase groundwater recharge, depending on interactions with transient climate (Wine et al 2015). Additionally, enhanced post-wildfire preferential flow may serve as another mechanism increasing groundwater recharge following wildfires (Stoof et al 2014). In snow-dominated regions, earlier soil thaw and longer periods of soil drainage may increase groundwater recharge following wildfire (Ebel et al 2012). Long-term watershed-scale response to wildfire integrates-over both space and time-many of these factors to determine the total wildfire impact on water yield. This work focuses on watershed-scale response to wildfire and assumes that hydrologic impacts of wildfires depend on-or are correlated with-vegetation recovery.
Despite recent evidence that wildfire increases regional water availability (Kinoshita and Hogue 2015, Bart 2016, Wine andCadol 2016), this hypothesis remains contrarian. Prevailing belief remains that land cover changes or ecological disturbances are largest at small spatial scales-plant, stand, and catchment-and are vanishingly small at larger spatial scales (Wilcox et al 2006). This phenomenon might be explained by a threshold effect in which a disturbance of some magnitude-commonly of at least 20% of the watershed area-is required to elicit a measurable hydrologic impact (Bosch and Hewlett 1982, Stednick 1996, Brown et al 2005, Wine and Cadol 2016, Buma and Livneh 2017, though the relative importance of measurement error versus a buffering effect remains unknown. In either case, should disturbances become more widespread the chances of being able to detect them should increase, even at larger spatial scales, contradicting conventional thought. However, the hypothesis that wildfires increase regional water yield in watersheds across the diverse ecoregions comprising the western USA remains untested. The goal of this work was to improve understanding of the relative importance of different global change components in influencing water yield. To advance this goal we tested the contrarian hypothesis that ecoregiondependent wildfire impact on long-term streamflow is larger relative to anthropogenic climate change impacts than previously thought-even at the large watershed scale (figure 2). We propose a simplified framework to test this hypothesis consisting of the following objectives: (1) Determine post-wildfire vegetationintended as an indicator of transpiration-recovery for each ecoregion division within the western USA. (2) Develop and validate a statistical model to estimate wildfire impacts on streamflow across all gauged wildfire impacted basins in the western USA. (3) Compare the relative magnitudes of wildfire and climate change impacts on streamflow, to test the hypothesis that the magnitude of wildfire enhancement suggests that wildfires have the potential to substantially offset climate change-induced reductions in streamflow.

Methods and study area
We developed a single framework to directly interrogate the importance of wildfires as a structural uncertainty affecting the accuracy of future water resource availability projections at the scale of the western USA. This framework integrates geospatial climate (Daly et al 1994), wildfire (Eidenshink et al 2007), streamflow, and moderate spatial resolution remotely sensed postwildfire vegetation recovery (Masek et al 2006) datasets from 39 watersheds across seven ecoregion divisions in the western USA (figures 3 and 4). We modeled the rate of asymptotic post-wildfire recovery of plant leaf biomass by ecoregion (table 1). We then used these ecoregion-specific vegetation recovery curves in a linear mixed-effects (LME) modeling framework to estimate the proportion of streamflow attributed to wildfires in all long-term gauged basins across the western USA, after controlling for climatic effects. Finally, we compare the wildfire impacts on streamflow estimated in this study to existing projections of climate change impacts on streamflow.

Experimental design
The objective of this study was to assess the hydrologic effects of wildfires across gauged basins of the western USA. We first selected all USGS reference watersheds 20% or more of which were burned within  Table 1. Ecoregion divisions included in this study. l is the characteristic recovery time scale parameter; larger values suggest slower post-wildfire vegetation recovery. Values in parentheses (where applicable) represent the number of burn scars from which post-wildfire vegetation recovery was derived from satellite imagery. Burned area refers the average proportion of a watershed burned (1986−2015); this value may exceed 100% implying that certain areas were burned repeatedly. The mean wildfire impact index over the study period is denoted as WII.  . These 39 watersheds range in size from small catchments to the large 6784 km 2 Yellowstone River watershed, thereby making the analysis directly relevant to regional water availability. (The larger watersheds may contain several distinct ecoregions and biomes due to spatial variation in environmental conditions.) Next, we compiled existing post-wildfire vegetation recovery curves and create new curves as necessary such that post-wildfire vegetation recovery curves will be available for each of seven western USA ecoregion divisions. From these recovery curves, we develop a wildfire impact index based on wildfire events and recovery curves. We then estimated annual precipitation within each watershed at a water year time scale from PRISM (Precipitation-elevation Regressions on Independent Slopes Model) (Daly et al 1994). At this point, all of the parameters required to fit a series of mixed effects statistical models are available, and these models are fit, validated, and used to estimate wildfire impact on streamflow. Lastly, we compare wildfire impact on streamflow to projections of change in streamflow by 2050 considering climate change alone.

Vegetation recovery
We assessed post-wildfire vegetation recovery, allowing this recovery to vary by ecoregion division (Bailey 1995 Typically, we considered recovery within all nonintersecting wildfire perimeters in Landsat scenes intersecting the 39 study watersheds. However, in the Marine Regime Mountains due to a high frequency of cloud-covered conditions, we considered post-wildfire recovery only within those wildfire perimeters that intersected study watersheds. During vegetation recovery variability may occur in vegetation index data due to climatic conditions-especially wet years that enhance plant growth or dry years that inhibit growth-thereby masking the overall recovery trend. Furthermore, the peak vegetation index is affected by satellite imagery availability. (In this analysis all data from Wine and Cadol (2016) are included as well as all available Level 1 at surface reflectance imagery from Landsat 5 (1984−2012) and Landsat 7 prior to scan line corrector failure (1999−2003). This amounts to analysis of ∼10 000 Landsat images altogether.) This availability in turn is influenced by the number of satellites that are in orbit and fully functional as well as by cloud cover.
Within each burn scar we calculated normalized vegetation recovery (NVR ) of the ith year post-fire as: where VI is the vegetation index during the ith postfire year, VI min is the lowest post-fire vegetation index value measured, and VI pre is pre-fire vegetation index. We iteratively fit normalized peak annual vegetation index data by burn scar to equation (2) and selected the median value for each ecoregion division, discarding those burn scars where departures from asymptotic recovery are particularly severe and numerical convergence is not achieved. We fit all remotely sensed recovery datasets with a Gauss-Newton algorithm (R Core Team 2017) to the simplified asymptotic function: where t is the number of full years that have elapsed since the fire and l is the characteristic recovery time scale parameter, which varies by ecoregion. (This simplified function masks certain vegetation dynamics such as tree mortality that lag wildfire events.) Finally, we calculated a wildfire impact index (WII) of year i in ecoregion j as: where BA is the proportion of the watershed burned during year t.

Watershed precipitation
To statistically control for precipitation effectswith the aim of estimating wildfire effects more accurately-we derive water year precipitation time series from PRISM, a leading measurement based spatial prediction approach well-suited to mountainous watersheds. For most analyses, we used a contiguous USA Albers equal area conic projection. All geospatial analyses were performed using ArcGIS 10.5; where automation was required the arcpy site package was used.

Statistical analysis
To test the hypothesis that wildfires significantly increase streamflow across the western USA, we created a series of increasingly complex mixed effects models (Pinheiro et al 2017) whose terms were specified a priori (table 2). To check that the model met relevant statistical assumptions we made a series of plots to check the residuals for normality, trends relative to fitted values, and trends relative to model terms after applying square root transformations to the precipitation and streamflow terms. The model that minimizes the Bayesian information criterion (BIC) is assumed to be optimal, though we also iteratively validate each model via n-fold cross-validation and present both sets of results.

Vegetation recovery
Over the course of the three-decade span of this study 320 wildfires burned forested areas within the study watersheds. Availability of pre-fire and seven consecutive years of post-fire imagery ranged from two burn scars in MaRM, which is commonly cloudy, to 273 in TDD. Where complete recovery trajectories were available, normalized vegetation recovery values range from 0-minimum vegetation index attained during the post-wildfire recovery period-to 2.5-representing enhanced vegetation vigor relative to pre-fire conditions. (Pre-fire vegetation vigor is assigned a value of unity.) In keeping with an assumption of asymptotic vegetation recovery, the lowest median normalized recovery values occur in the year immediately following wildfire events in TSuRM, TDD, and TStRM (figure 5). In contrast, in TSSD (and perhaps MaRM) recovery departs from monotonic-even at early time-with lowest recovery values in the second year following the wildfire event, consistent with delayed tree mortality. While absent inter-annual climatic variability, monotonic vegetation recovery might be expected, such variability-enhanced by the El Niño Southern Oscillation and Pacific Decadal Oscillation-causes departure from monotonic vegetation recovery. Thus, it is readily apparent that the satellite-derived vegetation recovery analyzed here represents a post-wildfire vegetation recovery signal superimposed on a dynamic climate signal. Fitting to this composite signal yields variability in recovery time scale parameters. The most rapid recovery occurs in TSSD, MD, and MeRM, whereas slowest rates of recovery occur in MaRM, TStRM, and TSuRM (table 1). is streamflow from the ith year of the jth watershed, as a function of precipitation (x 0 ), random intercepts by watershed ( 0 ), random slopes by watershed ( 1 ), ecoregion division ( ( ) ), and wildfire impact index (x 1 ). Models are assessed according to the Bayesian information criterion (BIC) and root mean square error (RMSE) calculated from n-fold (leave one out) cross-validation.

Wildfire impact index (WII)
In the absence of wildfire, the WII takes on a value of zero. In the year immediately following a fire this index increases by an amount equal to the proportion of the watershed burned. Following fires, WII decreases as an asymptotic function (figure 6) governed by the characteristic recovery time scale parameter (l ). Thus, the two key parameters governing WII are l and burned area (table 1). Burned area-as a proportion of watershed area-is largest in TSuRM, MaRM, and MD study watersheds and smallest in TDD, TStRM, and MeRM watersheds. In the framework presented here long-term hydrologic impact of wildfire is directly related to the area under the WII curve. The highest mean WII is associated with remarkably high burned area in TSuRM, whereas the second highest mean WII (MaRM) is associated with the largest recovery time scale (table 1). In contrast, the lowest WII occurs where recovery time scale is short or burned area is small, in TSSD and TDD, respectively.

Mixed effects modeling
To test the hypothesis that western USA wildfires significantly increase water availability across this region, we use this spatiotemporal variability in post-wildfire vegetation recovery as a quantitative predictor of streamflow. In doing so we fit and validate four sequentially more complex linear mixed effects models (table 2). Our results show that wildfires significantly improve prediction of streamflow across the western USA (p < 0.001), after controlling for precipitation and rainfall runoff relationships that vary by watershed ( figure 7, table 2). Model 4, the only model to include wildfire impacts yields the lowest BIC as well as the lowest n-fold (leave one out) cross-validation derived RMSE ( figure 8, table 2). The proportion of streamflow attributed to wildfire impacts varies substantially by ecoregion division (figure 9) due to spatial variation in wildfire frequency and post-wildfire recovery rate and ranges from 2%-14%. The largest wildfire contribution to streamflow occurred in the Tropical/Subtropical regime mountains ecoregion division, where dry lightning strike density is anomalously high (Abatzoglou et al 2016).
As a point of comparison climate change-based ensemble hydrologic predictions suggest a 1%-18% reduction in streamflow in western North America by 2050 (Milly et al 2005). This magnitude of wildfire contribution to streamflow-typically neglected in solely climate-based predictions of future water availabilityif considered would increase the structural uncertainty of predicted streamflow to 20%-370% (figure 10) of the magnitude of the value predicted by conventional climatically based approaches (Milly et al 2005). Hence, our results demonstrate that predicting future streamflow on the basis of climatic change alone rather than on the basis of global change-the integral of climatic and ecological changes-results in overly uncertain predictions of future water availability.

Discussion
While the results of this retrospective study have important implications for predictions of future water resource availability, the results presented here may differ from future contributions of ecological disturbances  1986−2015) to precipitation as modeled by PRISM (Daly et al 1994). In years affected by recent wildfires the difference between the orange and gray triangles indicates the fitted effect on streamflow attributed to wildfire influence.
to water resource availability. There is widespread agreement that wildfires will increase in certain regions of the western USA on the basis of diverse lines of reasoning including paleofire data (Marlon et al 2012), climatic trends that influence conditions conducive to wildfire (Jolly et al 2015), and integrated modeling of key controls on wildfire-climate, vegetation biomass, and ignition sources (Krawchuk et al 2009). Similarly, while pest-pathogen-complex-induced tree mortality currently elicits limited watershed-scale hydrologic impact (Penn et al 2016, Slinski et al 2016, anticipated future increases in the conditions conducive to die-off (Breshears et al 2009) could make this disturbance relevant to regional hydrology in the future. Conversely 'unparalleled' woody encroachment-estimated as high as 300 million ha in North America-would have an opposing influence on water yields relative to wildfire and die-off (Van Auken 2009). Consideration of ensembles of possible future disturbance scenarios and their hydrologic consequences in addition to changing climate is needed to transition from climate changeto global change-based hydrologic predictions, with the aim of reducing structural uncertainties inherent in model predictions. Inclusion of ecological disturbances may require modeling not just ecological dynamics, but also the economics and public policy that influence Anthropocene ecological disturbances. In addition to likely future increases in ecological disturbance frequency, the death of climatic stationarity (Milly et al 2008) also implies the death of the theory of ecological succession. Ecological succession is defined as 'the process by which a plant or animal community successively gives way to another until a stable climax is reached (Stevenson and Lindberg 2010)'. In the presence of a now non-stationary climate this theory of post-disturbance succession is moot. As a consequence, the simplified model of vegetation recovery fit in this study (figure 5) may tend to underestimate disturbance impacts by incorrectly assuming that plant communities will revert to the pre-fire climax community. Climax species in waterlimited regions are well adapted to maximize extraction of water from the subsurface (Seyfried et al 2005) consistent with the vegetal equilibrium or ecological optimality hypothesis (Eagleson 1978, Eagleson 1982, Eagleson and Tellers 1982, Milly and Eagleson 1987, Milly 1994, Milly and Dunne 1994, Kochendorfer and Ramirez 2010, O'Grady et al 2011, despite critiques of this hypothesis (Kerkhoff et al 2004). Following a disturbance under post-successional conditions, we speculate the series of plant communities that colonizes a site is likely to be less well adapted to the site due to the non-stationary climatic conditions, thereby potentially enhancing groundwater recharge.
Though the present work focuses on post-wildfire recovery and hydrologic impacts on the order of years to decades, these effects are superimposed on longer term forest dynamics-on the order of decades to centuries, which are particularly important from a water yield perspective in non-water limited watersheds, such as those in the Pacific Northwest. Long-term work on eucalypts in Australia observed relatively high water yield immediately following forest establishment, with decreasing water yields for 25 years followed by water yield increases until 200 years as stand leaf area index (LAI) decreased (Vertessy et al 2001). One mechanism behind this decrease in LAI and transpiration as a stand ages is now known to be hydraulic limitation, which increases as trees grow taller and branches grow longer, increasing hydraulic resistance (Murty et al 1996, Ryan andYoder 1997). Replacement of trees with low stomatal resistance by trees having high stomatal resistance has also been observed to increase streamflow as a forest ages post-disturbance under energy-limited conditions (Hornbeck et al 1997). Thus, the asymptotic vegetation recovery assumption employed here constitutes a source of uncertainty. Finally, an additional source of hydrologic complexity in post-wildfire recovery involves changes in transpiration rates of surviving trees. Thinning experiments-a good analog for low canopy burn severity in which a subset of trees survive the firehave demonstrated that remaining trees use more water (Whitehead et al 1984), thereby buffering the expected transpiration response. In the case of the Australian experience with wildfire-affected eucalypts, Nolan et al (2014) found that leaves resprouted along the trunk and branches, resulting in lower hydraulic resistance and hence higher transpiration per unit leaf area by nearly a factor of two. Therefore, despite the widespread practice of taking vegetation index as an index of transpiration (Nagler et al 2005, Wine and Hendrickx 2013, Cadol and Wine 2017, this approach is a simplification in areas in which the vegetation index-transpiration relationship changes following wildfire.

Future directions
In modeling wildfire impacts on water resources, vegetation recovery dynamics-modeled here as a nonparametric function-present a major source of uncertainty. Despite the evident importance of burn severity in influencing post-fire hydrologic response, hydrologic impacts of varying burn severity have often been neglected (Ebel 2013 As important as burn severity is in influencing postwildfire vegetation recovery, it interacts with climatic dynamics, land management strategies, and vegetation type. As all of these dynamics change in a global change context, future advances toward a more mechanistic approach that account for details of these processes will be necessary for accurate future predictions. Further research is needed to improve understanding of how western North American wildfires interact with the soil-plant-atmosphere continuum (SPAC) theory. According to this theory transpiration is limited by evaporative demand, plant canopy conductance, and soil water availability. Consequently, the hydrologic impacts of any disturbance will be mediated by these factors (Feikema et al 2013, Wine et al 2015, Wine andCadol 2016); further work developing a model that better accounts for non-linear interactions between ecological disturbances and SPAC dynamics is underway.

Conclusion
The results of this study suggest a need to revisit the role of land cover change and ecological disturbances in influencing streamflow. Ecological disturbances resulting in land-cover change are conventionally hypothesized to exert negligible influence on the water balance at increasingly large spatial scale (Blöschl et al 2007). However, this study, which includes five large watersheds that are each over 1000 km 2 in area demonstrates that the unprecedented scale of Anthropocene ecological disturbances has made these disturbances substantial and significant predictors of present and future regional water yields in ecoregions across western North America (figure 2). This study and past work (Wine and Cadol 2016) demonstrate that at small spatial scales wildfire contribution to streamflow is intermittent, whereas if a threshold of aerial extent continues to be exceeded at larger spatial scales this patchwork of wildfire histories integrates to sustained increases in streamflow due to wildfire effects that may rival or exceed projected climate change impacts.