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

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 and Cadol 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 and Cadol 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 ecoregion-dependent 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 vegetation—intended 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.

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 post-wildfire 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.1. Experimental design

2.2. Vegetation recovery

We assessed post-wildfire vegetation recovery, allowing this recovery to vary by ecoregion division (Bailey 1995). Where possible we used existing post-wildfire vegetation recovery curves, including in the Mediterranean Division (Hope et al 2007) and the Mediterranean Regime Mountains (McMichael et al 2004), or existing post-wildfire vegetation recovery datasets (Wine and Cadol 2016). Where these were unavailable we extracted peak annual time series of normalized difference vegetation index (NDVI) from remotely sensed Landsat-derived Level I at surface reflectance imagery (Masek et al 2006). NDVI is an effective means of assessing plant leaf biomass in that this band ratio attains a maximum value when chlorophyll absorption of red light peaks and healthy mesophyll tissue reflection of infrared light decreases (Campbell and Wynne 2011).

Typically, we considered recovery within all non-intersecting 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_i) of the ith year post-fire as:

$$\begin{equation} \rm NVR_{\it i} = \frac{{V{I_{\it i}} - V{I_{min}}}}{{V{I_{pre}} - V{I_{min}}}} \end{equation} \tag{ 1 }$$

where VI_i is the vegetation index during the ith post-fire 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:

$$\begin{equation} 1 - {{{\rm e}}^{ - \frac{{{t}}}{{{{{l}}_{{j}}}}}}} \end{equation} \tag{ 2 }$$

where t is the number of full years that have elapsed since the fire and l_j 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:

$$\begin{equation} {\rm{WI}}{{\rm{I}}_{{\rm{\it ij}}}} = \mathop \sum \limits_{{{t}} = 1984}^{{{t}} = {{i}} - 1} {\rm{B}}{{\rm{A}}_{{t}}} \cdot {{{\rm e}}^{ - \frac{{{{i}} - 1 - {{t}}}}{{{{{l}}_{{j}}}}}}}, \end{equation}$$

where BA_t is the proportion of the watershed burned during year t.

2.4. Watershed precipitation

2.4. Statistical analysis

3.1. 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).

Table 2. Streamflow models of increasing complexity, where y_ij is streamflow from the ith year of the jth watershed, as a function of precipitation (x₀), random intercepts by watershed (u_0j), random slopes by watershed (u_1j), ecoregion division (α_k(j)), and wildfire impact index (x₁). 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.

Model #	Model form	BIC	RMSE
1	${{{y}}_{{{ij}}}} = {{{\beta }}_{{o}}} + {{{\beta }}_1} \cdot {{{x}}_0}_{{{ij}}} + {{{u}}_0}_{{j}} + {{\rm{e}}_{{{ij}}}}$	5485	110.63
2	${{{y}}_{{{ij}}}} = {{{\beta }}_{{o}}} + {{{\beta }}_1} \cdot {{{x}}_0}_{{{ij}}} + {{{u}}_0}_{{j}} + {{{u}}_1}_{{j}} \cdot {{{x}}_0}_{{{ij}}} + {{\rm{e}}_{{{ij}}}}$	5225	94.15
3	${{{y}}_{{{ij}}}} = {{{\beta }}_{{o}}} + {{{\beta }}_1} \cdot {{{x}}_0}_{{{ij}}} + {{{\alpha }}_{{{k}}\left( {{j}} \right)}} + {{{u}}_0}_{{j}} + {{{u}}_1}_{{j}} \cdot {{{x}}_0}_{{{ij}}} + {{\rm{e}}_{{{ij}}}}$	5217	94.07
4	${{{y}}_{{{ij}}}} = {{{\beta }}_{{o}}} + {{{\beta }}_1} \cdot {{{x}}_0}_{{{ij}}} + {{{\beta }}_2} \cdot {{{x}}_1}_{{{ij}}} + {{{\alpha }}_{{{k}}\left( {{j}} \right)}} + {{{u}}_0}_{{j}} + {{{u}}_1}_{{j}} \cdot {{{x}}_0}_{{{ij}}} + {{{\rm e}}_{{{ij}}}}$	5140	91.09

Zoom In Zoom Out Reset image size

3.2. 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_j). Thus, the two key parameters governing WII are l_j 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.

Zoom In Zoom Out Reset image size

3.3. 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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 availability—if 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.

Zoom In Zoom Out Reset image size

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 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 change- to 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 water-limited 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 and Yoder 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 fire—have 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.

In modeling wildfire impacts on water resources, vegetation recovery dynamics—modeled here as a non-parametric 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, Bart 2016, Wine and Cadol 2016). As important as burn severity is in influencing post-wildfire 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 and Cadol 2016); further work developing a model that better accounts for non-linear interactions between ecological disturbances and SPAC dynamics is underway.

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² 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.

We wish to acknowledge funding from New Mexico EPSCoR WC-WAVE NSF award 1329470, a New Mexico Geological Society grant-in aid, Geological Society of America, the New Mexico Water Resources Research Institute, and the United States–Israel Educational Foundation—the Fulbright Commission in Israel. The lead author is currently supported by a post-doctoral fellowship at Ben Gurion University. We also gratefully acknowledge Josh Malinow's contribution in designing figure 1. Finally, we thank two anonymous reviewers for comments that greatly improved the initial manuscript.