CHARACTERIZATION OF POST-FIRE STREAMFLOW RESPONSE ACROSS WESTERN U . S . WATERSHEDS

This research investigates the impact of wildfires on watershed flow regimes, specifically focusing on evaluation of fire events within specified hydroclimatic regions in the western United States. Information on fire events and watershed characteristics were collected through federal and state-level databases and streamflow data were collected from U.S. Geological Survey stream gages. Eighty two watersheds were identified with at least ten years of continuous pre-fire daily streamflow records and five years of continuous post-fire daily flow records. For each watershed, percent change in annual runoff ratio, 5 low-flows, high-flows, peak flows, number of zero flow days, baseflow index, and Richards-Baker flashiness index were calculated using preand post-fire periods. The gathered watersheds were divided into nine regions or clusters through k-means clustering and regression models were produced for watersheds grouped by total area burned. The coefficient of determination (R2) was used to determine the accuracy of the resulting models. Results show that low flows, high flows, and peak flows increase significantly in the first two years following a wildfire and decrease over time. Relative response was utilized to scale 10 response variables with respective percent area of watershed burned in order to compare regional differences in watershed response. Watersheds in Cluster 9 (eastern CA, western NV, OR) typically demonstrate a negative relative post-fire response, in that when scaling response to area burned, a slight negative response is observed in flow regimes. Most other watersheds show a positive mean relative response. In addition, regression models show limited correlation between percent watershed burned and streamflow response, implying that other watershed factors strongly influence response. 15


Introduction
The number of wildfires in the western United States (US) is increasing annually, on average costing federal agencies billions of dollars a year in suppression efforts (Whitlock, 2004) and causing increases in flood events destructive to both life and infrastructure in many parts of the world (Daniel G Neary, 2003;Juli G. Pausas, 2008).Westerling et. al. (2006) showed that the western fire regime exhibited a significant transition from infrequent and short-duration events to higher frequency, longer duration regimes during the mid-1980's.The greatest increases in fire frequency were found to occur in mid-elevation forests, most commonly in the Northern Rockies, Sierra Nevada, southern Cascades, and western Coast Ranges in northern California and southern Oregon (Littell et al., 2009).This marked change is strongly correlated with climate change impacts, such as warmer springs and longer dry seasons, commonly in occurrence with reduced winter precipitation rates and earlier spring snowmelt.Overall, Westerling et. al. 2006 determine that, though land-use history may be a significant factor in the spatial distribution of wildfires within specific forest types, changes in fire regimes in the western US can most likely be attributable to recent changes in climate.Other notable research has also provided significant correlatory evidence between climate change and wildfire occurrences (Littell et al., 2009;Moritz et al., 2010).
Though wildfires are a part of the natural process of vegetation dynamics, they cause wide-ranging changes to ecosystems (Daniel G Neary, 2003;Santos et al., 2015) depending on numerous factors, most importantly burn severity.Studies examining the effects of wildfires on a small-scale, such as in plot-sized and laboratory experiments, show high fire temperatures can result in the combustion of organic matter within soils and cause permanent alteration to the chemical structure of local clays, decreasing soil stability (Shakesby and Doerr, 2006).Water-repellent soil layers can be created in a discrete layer on or below the soil surface through chemical bonding of the combusted organic matter to mineral particles, potentially increasing overall topsoil erosion rates in burned regions (Wilkinson et al., 2009), though this hydrophobicity is highly variable depending on fire behavior, burn severity and soil properties (DeBano, 2000).
At larger scales, such as entire watersheds or multiple watershed systems, studies of post-fire erosion rates have shown incompatible conclusions (Moody and Martin, 2001;Owens et al., 2013;Smith et al., 2011), though this is most likely due to the variability of precipitation events and general climate patterns (Moody et al., 2013).In terms of water quality, contaminant levels can be dramatically increased for many years after a wildfire in both soil (Burke et al., 2010) and stream systems (Emelko et al., 2011;Stein et al., 2012;Burke et al. 2010), increasing the workload on source water protection organizations in communities reliant upon burned watersheds for drinking and farm water.Furthermore, wildfires are readily attributed as the cause of substantial increases in debris flows (Benavides-Solorio and MacDonald, 2001;Cannon et al., 2001;Meyer et al., 2001).
Studies evaluating post-fire water yield change are highly disparate owing to the transient nature of climate patterns, variations in basin geomorphology, and vegetation recovery patterns, and the resulting complex interactions (Moody et al., 2013).
For example, studies in rangeland regions of the US found moderate increases in flow, infiltration, and erosion rates after major wildfires, with trends continuing for as long as 15 years (Emmerich and Cox, 1994;Frederick B. Pierson, 2009;Hester et al., 1997).Fires in chaparral environments, such as in southern California, exhibited increased flows up to as much as two orders of magnitude, with much of this occurring in the dry season (Coombs and Melack, 2013;Kinoshita and Hogue, 2015;Loáiciga et al., 2001).Fires in other chaparral environments were found to also yield flow increases, such as in South Africa (Lindley et al., 1988;Scott, 1993), Cyprus (Hessling, 1999), and France (Lavabre et al., 1993).Additional increases to post-fire flow regimes were found in temperate, forested catchments as well (Neary et al., 2005;Watson et al., 2001).A concise summary of historic changes in US post-fire stream systems is found in Neary et al. (2005), documenting changes in 1st year runoff and peak flows, encompassing a range of ecological regions.Conversely, several studies found limited or no significant changes to hydrologic systems post fire, or attributed fluctuations to natural annual variability (Aronica et al., 2002;Bart and Hope, 2010;Britton, 1991;Townsend and Douglas, 2000).
These discrepancies in post-fire flow response lead to the question of which watershed characteristics have the greatest influence over the observed response?Moody et al. (2013) provide a succinct summary of soil-related theories, such as reduced infiltration due to increases in soil-water repellency, increased overland flow velocities due to increased bare ground, and reduced infiltration caused by soil-sealing.Theories commonly found in literature attribute flow changes to a wider range of factors, including reduction in interception and evapotranspiration (Lavabre et al., 1993;Scott, 1993) and increased hydrophobicity of soils (Neary et al., 2005).In regards to altered peak flows, conflicting evidence is found regarding the importance of burned watershed areas with some studies finding an inverse correlation between peak flows and watershed size (Biggio and Cannon, 2001;Neary et al., 2005) and others finding no relation at all (Bart and Hope, 2010).
The current study builds upon prior studies but undertakes a more comprehensive assessment of post-fire streamflow changes in the US by examining burned watersheds that encompass a wide spectrum of climatological and geophysical parameters.A variety of flow parameters are also examined which describe changes to flow regimes at several levels.Furthermore, the variability in response by distinct regions is investigated, anticipating distinct differences influenced by regional climate.With downstream communities at risk for flooding, and also relying on catchment runoff for water supply, investigating alterations in post-fire discharge over large scales will provide critical information for regional managers on post-fire runoff mitigation.In addition, understanding factors controlling discharge response will help inform development and calibration of surface water models used for post-fire streamflow predictions.

Study Areas
A total of 82 burned watersheds were utilized (Fig. 1), encompassing a wide range of spatial, temporal, climatological, and topographic factors (Fig. 2).These watersheds were exclusively limited to those with significant wildfires (burned > 5%) and adequate (continuous 15 years daily flow) discharge records available in the USGS streamflow database (U.S. Geological Survey, 2014) identified through the GAGES-II database (Falcone, 2011).The majority of available watersheds are overwhelmingly found in the western US, predominantly in California, Oregon, and Idaho, with several located in the North East, Florida, and Kentucky.Due to discharge and burn severity data limitations, the fires in this study cover a temporal range from water years 1984 through 2010.Average percent area burned ranges from 5-97%, with a mean of 25%, over a range of watershed areas from 4.6-9209 km2.The wide spatial distribution of the studied watersheds results in mean elevations and burn-area slopes varying from 13-2760 m and 0.11-16% respectively (Fig. 2).
The most important difference between many of these watersheds is the variation in climate, the values of which were collected from the GAGES-II dataset (Falcone, 2011), which catalogs all watersheds in the US monitored by the USGS and have at least one twenty year period with continuous daily flow records.Average basin precipitation ranges from 29-220 mm/yr, with a mean of 72 mm/yr, and average temperature ranges from 1.4 -23 °C, with a mean of 10 °C.Important for identifying snow dominated regions is the percent of precipitation (PPT) that falls as snow (%Snow/PPT), which ranges from 0-72%.
Relative humidity ranges from 39-73%, with a mean of 55, and potential evapotranspiration ranges from 400-1200 mm/yr, with a mean of 633 mm/yr.
Watershed vegetation types vary across the watersheds as well.Evergreen forest and shrub vegetation are the overwhelmingly dominant land cover type over all watersheds used in this study (Fig. 3).The high proportion of evergreen is due to the

Data Collection
Watershed parameters utilized in this study were chosen to encompass the variability in geophysical parameters found throughout the watersheds used in this study (Table 2.1).To identify spatial trends in post-fire response, watersheds were first grouped through k-means clustering based on geographic and climatological data (described in section 3.3 below).

Watershed selection
Watersheds were selected based upon the continuity of USGS mean daily flow records available before and after the fire event.
Watersheds were required to have continuous daily flow records (>95% of daily flow records accounted for in each year) for a minimum of 10 years pre-fire and 5 years post-fire.Using the approximately 9,000 watersheds in the GAGES-II dataset (Falcone, 2011), delineated watersheds were spatially cross-referenced with the MTBS database of historic wildfires (2009).
The results were again cross-referenced with USGS daily flow records (U.S. Geological Survey, 2014) to identify watersheds with the required flow records, resulting in 263 unique watersheds in the US with greater than 5% total burn area in a single water year.Of these, 23 contained 2-3 wildfires within the same year burning over 5% of the total area.The remainder contained only a single significant fire in the year of interest.Further exclusion of watersheds was based on the presence of major dams within the watershed flow regimes extracted from the GAGES-II database (Falcone, 2011), resulting in a final collection of 82 watersheds.

Streamflow and Precipitation Data
Daily flow and peak flow data were obtained from the USGS National Water Information System (U.S. Geological Survey, 2014), matched for the range of 10 years pre-fire and 5 years post-fire for each selected watershed.Monthly precipitation data was collected from the PRISM database (PRISM Climate Group, 2004), a nation-wide 4 km resolution gridded monthly dataset that extrapolates station climate measurements over unmonitored areas using a complex topographic-and climatebased algorithm.Monthly national precipitation rasters were averaged for each watershed for all months within the flow record period.

Climatological data
Watershed climatological parameters used in this study included percent of precipitation that falls as snow (%Snow/PPT) and the aridity index (AI).The %Snow/PPT for each watershed was available through the GAGES-II dataset (Falcone, 2011) and the aridity index was calculated for each watershed as: where P avg is average precipitation and PET avg is average potential evapotranspiration, both of which were available in the GAGES-II dataset.

Response Variables
A range of response values were selected to quantify post-fire flow changes across a variety of regimes, including flows relating to dry seasons (low flows, base flows) and wet seasons (high flows, peaks flows).

Low, high, and peak flows
Pre-fire low-flow (LF) and high-flow (LF) metrics were calculated for each of the ten years prior to the fire water year and averaged to produce a single value.Low flows (LFs) were defined as the average of mean daily flows with a 90% exceedance within a single water year.To reduce calculation bias due to zero flow days commonly found in ephemeral stream systems, zero flow days were eliminated from exceedance value calculations.High flows (HFs) were defined similarly, with a 10% exceedance threshold used to isolate larger volume flows (Kinoshita and Hogue, 2015).Changes in LFs and HFs were calculated as the post-fire percent change from the average 10 water years pre-fire.Post-fire values were calculated for the 1st year (LF.one,HF.one), the 2nd year (LF.two,HF.two), and the 5 year mean (LF.five,HF.five).
Peak flows (PFs) were defined as the largest mean daily flow measurement each water year and post-fire changes in PF were calculated as the percent change of the 1st year (PF.one), 2nd year (PF.two), and 5 year mean (PF.five) peak flow measurements from the pre-fire ten year mean.Percent changes in the number of zero flow days were calculated similarly (Nzero.one,Nzero.two,Nzero.five).

Runoff Ratios
The runoff ratio (RO) is defined as the fraction of total annual runoff depth over total annual precipitation: where P tot is total annual precipitation, Q tot is total annual runoff depth, and A ws is watershed area.RO was calculated for the ten years pre-fire and 5 years post-fire using PRISM precipitation and USGS mean daily flow data.Post-fire RO response was calculated as the percent change between the 1st year, 2nd year, and average 5 year values post-fire (RO.one,RO.two, RO.five) and the pre-fire 10 year mean.

Base flow and Richards-Baker indices
Base flow index (BFI), defined as the fraction of total streamflow that is baseflow (Baker et al., 2004), was calculated for each water year through the R package 'hydrostats' (Bond, Nick, 2015), that applies the Lyne-Hollick filter(V.D. Lyne, 1979).BFI response was calculated as the percent change of the 1st, 2nd, and 5 years BFI post-fire from the mean of the 10 years pre-fire BFI (BFI.one,BFI.two, BFI.five).
The Richards-Baker index (RB) quantifies the frequency and rapidity of short-term changes in streamflow (flashiness) based on daily flow data through the equation: where q is mean daily flow, t is time, and q is daily flow (Baker et al., 2004).RB response was calculated as the percent change from the average RB over the ten years pre-fire to the RB of 1 year, 2 years, and average of 5 years post-fire (RB.one,RB.two, RB.five).

k-means Clustering
To better understand regionalize differences in post-fire flow response and create region-specific regression models, watersheds were classified into unique regions through k-means clustering (MacQueen, 1967) based on all GAGES II watersheds.This approach partitions an N-dimensional population of observations into clusters with minimal variation, allowing for relatively simple similarity grouping.For the current study, the ideal ensemble of clusters was one that produced easily recognizable regions with unique climatological characteristics.Large-scale clustering methods have been applied in prior watershed classification studies, but utilized more complex streamflow and ecological indices as parameters (McManamay et al., 2014;Poff, 1996).
Wildfires in this study were typically found in western evergreen and shrub environments, so clustering by only these watersheds would likely produce regions biased by fire occurrence.To limit this impact, we applied the mclust package in R Information Critera (BIC) to determine the most accurate model and cluster count (Schwarz, 1978).Various groupings of simple parameters were used for clustering including watershed latitude and longitude, elevation, AI, %Snow/PPT, and mean monthly and seasonal flow statistics.

Clustering
The k-means clustering performed on the GAGES-II watershed set yielded 9 clusters or regions (Fig. 4).The most important clusters for the current study are 6 through 9, which assemble 77 of the 82 watersheds into unique regions (Fig. 5).These four clusters have unique characteristics (Fig. 6).Watersheds in cluster 6, on average, have the highest Kfact, though almost all are burned less than 20%.They have relatively moderate %Snow/PPT values and the lowest AI values.Watersheds in cluster 7 have the highest average area and elevation, and accordingly the highest %Snow/PPT and the lowest NDVI.Cluster 8 contains watersheds with the widest range of relative fire sizes, including watersheds burned from as little as 10% to as great as 97%.
Watersheds in cluster 8 also have the lowest average elevations and areas, as well as the smallest %Snow/PPT and low AI.The percent of the burn area rated as high burn severity is also the greatest on average in cluster 8 watersheds.Cluster 9 watersheds have the lowest Kfact and highest elevations.These watersheds also have high %Snow/PPT and AI.

Response Variable Distribution and Analyses
Calculated flow response variables indicate an extremely wide range of post-fire system responses (Fig. 7).The greatest ranges occur within variables representing changes in low flows, such as Nzero.one(st.dev = 243%), BFI.one (236%) and Nzero.five(201%).The tightest ranges are typically found within mean five year variables where extreme changes are muted, such as RB.five (28%), RO.five (39%), and BFI.five (48%).Response variable means range from as low as 0.92% (RB.five) to as great as 115% (Nzero.one).Due to the nature of response calculations (section 3.2.1),Nzero values were limited.Nzero.one and Nzero.two were found for 22 watersheds and Nzero.five was found for 27 watersheds.

Trend analysis
Comparison of 1st year response variables to 5 year mean variables generally produce line-of-best-fit slopes greater than one, with a mean slope of 1.9 (Fig. 8).The greatest slopes are found in LF.one, HF.one, and BFI.one.BFI.one is unique in the magnitude of its slope versus BFI.five, which is twice as large as that of LF.one.Only Nzero.oneversus Nzero.fiveyields a slope less than one.Lines-of-best-fit of 2nd year response variables versus 5 year mean variables produce significantly lower slopes, with a mean of 0.78.Only LF.two and RB.two yield greater values.We can infer from this that the greatest increases in these response variables in the five water years following a fire are found in the first year.Additionally, only LFs exhibit greater response in the second year than in the following three years.
Comparing response variables to percent area burned, lines-of-best-fit shows the greatest slopes in 1st year responses for LF, HF, and PF variables (Fig. 9).In the case of RO, lines exhibit similar increases between RO.one and RO.two.Typically, 2nd year slopes tend to be steeper than 5 year mean values.The exception to this is Nzero.two,where 2nd year values decrease significantly with increased fire size.Overall, Nzero is found to increase post-fire, though due to both a small sample size and a short time period, results are most likely uncertain.BFI increases with increasing burn size, though to a lesser extent than the previously mentioned variables.Only RB indicates little linear correlation to burn area, with marginal 1st and 2nd year slopes.In fact, RB.five decreases with increasing percent area burned.Overall, these findings confirm previous smaller-scale studies in which the greatest flow responses occur immediately after fire events and decrease with time.

Response variability by cluster
To simplify an analysis of response variables, boxplots are provided comparing responses across the four significant clusters noted above (Fig. 10).In this instance and that of the CONUS plots, all variables are scaled by dividing the variable by the Watersheds in cluster 8 (found in CA) differ in overall magnitude with respect to both the standard deviations and the mean values of the other clusters.Standard deviation within relative response variables is significantly higher than that of clusters 7 and 9, though lower than that of cluster 6. Notable instances of this are observed in the cases of PF.one (st.dev = 11.7%) and BFI.one (16.4%).Mean values are the second largest among the clusters, close in magnitude to cluster 6.  Variability in response variables also generally decreases with increasing percent watershed burned.Linear regression modeling of response variables by percent watershed burned yields the error statistics found in Figure 11

Conclusions
Post-fire changes in streamflow are found to be highly variable across regions of the western U.S. and some trends can be difficult to discern.In general, flow response was found to be greatest in both magnitude and variability within the first year following wildfire and was shown to decrease over a five year period.The most important trends in post-fire response were found within low flow, high flow, and peak flow variables.This study found that these variables increase significantly in the first and second years following wildfire.On a regional basis, the following trends were identified among relative response variables: -Cluster 6 (AZ/NM/UT/CO/MT/WY/SD) -greatest variability and lowest magnitude response.Positive response in first year and five year means.
-Clusters 7 (UT/NV/WY/ID/MT/OR/WA) -lowest variability and a positive mean response.Positive response in second year and five year means.
-Cluster 8 (CA) -second-largest mean relative response.Positive response in first year and five year means.
-Cluster 9 (OR/CA/NV) -lowest variability and a negative mean relative response.Almost entirely negative response across the temporal scale.
-When examining a collection of burned watersheds of this scale, five year mean values consistently show a positive response across clusters 6-8.
Importantly, correlation between percent of a watershed burned and flow response was only significant among some of the first year response variables, and only when including watersheds burned greater than 50%.Much of the response variability seemed to be produced within watersheds burned less than 40%.Specific watershed properties also have a significant influence on post-fire streamflow response, though some of the variability is undoubtedly due to transient climate patterns.A complementary paper is forthcoming that will incorporate multiple regression models to identify how various geophysical parameters control flow response.The authors predict that burn severity, vegetation, and slope will have a strong influence on streamflow response.Results from this study will assist water resources organizations in developing post-fire water budgets, as well as highlight regions where emergency services may be required to increase their funding and preparedness.

Figure 1 .
Figure 1.CONUS map of the locations of the 82 watersheds utilized in this study.

Figure 2 .
Figure 2. Boxplots summarizing the range of watershed areas, elevations, percent area burned, percent of precipitation that falls as snow, and aridity index.

Figure 4 .
Figure 4. CONUS map of the results of k-means clustering of the GAGES-II watershed set.

Figure 5 .
Figure 5. CONUS map of the results of watershed clustering.

Figure 6 .
Figure 6.Variations in explanatory variables by cluster.

Figure 7 .
Figure 7. Boxplots of the response variables utilized in this study.

Figure 8 .
Figure 8. Scatterplots of 1st and 2nd year response variables versus 5 year mean response variables.

Figure 9 .
Figure 9. Scatterplots of all response variables versus respective watershed percent area burned.Simple regression lines provided to show correlations. 10

Figure 10 .
Figure 10.Boxplots of relative response variables by cluster.

Figure 11 .
Figure 11.Scatterplots of the R 2 and p-value results of simple linear regression modeling of response variables versus an increasing percent burn limitation.Squares are R 2 and circles are p-values.Filled squares are R 2 >-0.50, filled circles are p <= 0.05.The blue dashed line is R 2 = 0.50, the red dotted line is p = 0.05.
(i.e.decreasing sample size with increasing percent area burned).Figures include the adjusted R 2 and p-value significance tests (alpha = 0.05).Generally, linear modeling of 1st year response variables increases in accuracy as included values are limited by increasing percent burn area.LF.one,HF.one, PF.one, and RB.one show substantial increases in adjusted R 2 once included watersheds are limited to those exceeding a 50% burn area (n=12).Included p-values indicate that several of the LF.one and PF.one models are statistically significant.Applying the same methods to 2nd year values produces dissimilar results, with adjusted R 2 exceeding 0.5 in only a single instance (HF.two), and few significant p-values.However, in the cases of LF.two and RB.two, R 2 values increase with increasing percent burn threshold.Unsurprisingly, simple regression of 5 year mean values versus percent area burned produce mixed results with only LF.five and PF.five allowing for adjusted R 2 values greater than 0.5, few of which are statistically significant.BFI.five shows a single instance of a high adjusted R 2 value but is statistically insignificant.Simple regression modeling was also performed on response variables by limiting included watersheds by decreasing percent area burned (i.e.decreasing sample size with decreasing percent area burned) and results demonstrated zero significant adjusted R 2 values.