Abstract
In recent years climate change and historic fire suppression have increased the frequency of large wildfires in the southwestern USA, motivating study of the hydrological consequences of these wildfires at point and watershed scales, typically over short periods of time. These studies have revealed that reduced soil infiltration capacity and reduced transpiration due to tree canopy combustion increase streamflow at the watershed scale. However, the degree to which these local increases in runoff propagate to larger scales—relevant to urban and agricultural water supply—remains largely unknown, particularly in semi-arid mountainous watersheds co-dominated by winter snowmelt and the North American monsoon. To address this question, we selected three New Mexico watersheds—the Jemez (1223 km2), Mogollon (191 km2), and Gila (4807 km2)—that together have been affected by over 100 wildfires since 1982. We then applied climate-driven linear models to test for effects of fire on streamflow metrics after controlling for climatic variability. Here we show that, after controlling for climatic and snowpack variability, significantly more streamflow discharged from the Gila watershed for three to five years following wildfires, consistent with increased regional water yield due to enhanced infiltration-excess overland flow and groundwater recharge at the large watershed scale. In contrast, we observed no such increase in discharge from the Jemez watershed following wildfires. Fire regimes represent a key difference between the contrasting responses of the Jemez and Gila watersheds with the latter experiencing more frequent wildfires, many caused by lightning strikes. While hydrologic dynamics at the scale of large watersheds were previously thought to be climatically dominated, these results suggest that if one fifth or more of a large watershed has been burned in the previous three to five years, significant increases in water yield can be expected.
Export citation and abstract BibTeX RIS
1. Introduction
In the mid-1980s large wildfires in western North American forests increased markedly in spatial extent, duration, frequency, and severity in association with higher spring and summer temperatures, as well as reduced winter precipitation (Westerling et al 2006, Miller et al 2008, Williams et al 2010, Jolly et al 2015). This regional increase in large wildfires occurred in association with an unprecedented multi-year drought that may have been a consequence of climate warming and a harbinger of a prolonged intensification of aridity in this region (Cook et al 2004, Seager et al 2007). In addition to prolonged drought, another factor contributing to this increase in large wildfires is a history of fire suppression in the western United States that has resulted in a 'fire deficit' relative to long-term patterns (Marlon et al 2012). Examination of historical periods of climate warming implies that conditions conducive to large wildfires will continue as Earth's climate warms further (Calder et al 2015).
In the intermountain region of the western US wildfires interact with the North American monsoon. This phenomenon brings thunderstorms accompanied by high intensity rainfall during mid to late summer (Underwood and Schultz 2004). Such high intensity convective precipitation may interact with post-wildfire conditions to produce flash flooding (Underwood and Schultz 2004). Wildfires combust live plant canopy cover and leaf litter and increase the severity of soil water repellency (Mataix-Solera and Doerr 2004, Granged et al 2010, Versini et al 2012). Canopy removal tends to increase throughfall, which—after combustion of the leaf litter layer—promotes surface sealing, and reduces infiltration (Shakesby et al 2000, Savadogo et al 2007, Larsen et al 2009), thereby potentially modifying monsoon season hydrologic response by increasing runoff ratio, peak discharge, and flood frequency (Mataix-Solera and Doerr 2004, Granged et al 2010, Versini et al 2012). Commonly the extent to which post-wildfire infiltration-excess overland flow increases is similar across burn severities (Cawson et al 2013, Vieira et al 2015).
In addition to effects of fire on soil hydraulic properties that influence partitioning of high intensity monsoonal precipitation between infiltration into the subsurface and infiltration-excess overland flow, fire also may reduce transpiration (Kinoshita and Hogue 2011, Ebel 2013, Cardenas and Kanarek 2014, Zhou et al 2015). Kinoshita and Hogue (2015) examined dry season water yield in 14–51 km2 southern California watersheds following wildfires that burned 87%–95% of the watershed area. They observed 118%–1090% increases in low flow volumes following the 2003 Old Fire and attributed this increased discharge to basin-wide reduction in transpiration resulting from plant canopy removal. Based on these results they suggest that such increases in baseflow may supplement regional water supplies that are often stressed in semi-arid regions.
The hydrologic effects of wildfire are, however, modified by spatiotemporal scale. At the point, hillslope and small watershed scales many effects of fire are well known. In contrast, at the scale of a large watershed, fire may be only one process modifying watershed function, particularly if only a portion of the watershed has been burned. Fire effects on fluxes measured at the watershed outlet may be modified by such processes as climatic anomalies, 'extreme' changes in geomorphology following the fire due to debris flows (Pelletier and Orem 2014, Orem and Pelletier 2015), channel routing, overbank flooding, flood wave attenuation, and the degree to which convective precipitation events within large watersheds coincide with burned areas. With regard to temporal scale, over long-term periods the hydrologic effects of wildfires represent only one potential influence on hydrological processes. Others include natural climate variability (Cook et al 2007, Wine et al 2012b, 2015), anthropogenic climate change (Cayan et al 2010), and land-cover change (Wilcox et al 2008, Wine and Zou 2012). Commonly hydrologic studies of wildfire effects adopt short time horizons and small spatial extents (Shakesby and Doerr 2006). Examples of such studies range from simulating rainfall on runoff plots for several minutes (Granged et al 2010) to measuring post-fire discharge from watersheds up to 100's of square kilometers in area following individual fires for periods of up to seven years (Kinoshita and Hogue 2015, Bart 2016).
Therefore, this study aims to improve understanding of the importance of fire effects on hydrologic response of large semi-arid mountainous watersheds in the southwestern USA relative to long-term climate variability. Our long-term goal is to fulfill this aim through the use of coupled surface-subsurface, physically based hydrologic models. However, many commonly used hydrologic models are ill-suited to modeling post-fire hydrology in which infiltration rates may initially increase at early time during a rainfall event—as amphiphilic molecules within the topsoil transition from unwettable to wettable orientations—and subsequently decrease with time—as the soil approaches saturation (Imeson et al 1992, Wine et al 2012a). In this study we therefore adopt a statistically based approach to interrogate the relative importance of candidate controls on hydrologic response of large watersheds subjected to wildfires.
2. Methods
2.1. Site selection
To elucidate the relationship between wildfires and streamflow, we obtained geospatial wildfire occurrence data developed by the Monitoring Trends in Burn Severity (MTBS) project (Eidenshink et al 2007) and available from the launch of Landsat 5 through the present (1984–2014). Since MTBS only extends back to 1984, we supplemented the MTBS data with federal wildland fire occurrence data compiled by the USGS starting in 1980. We then selected three watersheds in New Mexico that each experienced two or more large wildfires, were unaffected by manmade hydraulic control structures, and had continuously operating stream gages over the period of interest from 1982 to 2014 (figure 1).
Figure 1. Geospatial analysis revealed that the Jemez (top right), Mogollon (bottom right), and Gila (bottom right) watersheds are well-suited for assessment of the hydrologic relations associated with wildfires incident on southwestern USA watersheds because these three watersheds are natural—unaffected by major channel flow control structures, have long-term stream gauges at their respective outlets, and have been profoundly impacted by large wildfires. Burn severity fields from the MTBS project represent wildfires that burned between 1984 and 2014; those burn severity fields derived from post-scan-line-corrector-failure Landsat 7 imagery exhibit striping that is most apparent in the Jemez watershed.
Download figure:
Standard image High-resolution image2.2. Site characteristics
The Jemez watershed intersects two ecoregion divisions—tropical and subtropical regime mountains to the west and steppe regime mountains to the east—and two ecoregion provinces—Arizona-New Mexico mountains semi-desert-open woodland-coniferous forest-alpine meadow province to the west and southern Rocky Mountain steppe-open woodland-coniferous forest-alpine meadow province to the east (Omernik 1987). Spatially weighted mean monthly precipitation estimated from PRISM (1982–2014) typically reaches a low in June (20 mm) throughout the watershed, then peaks in July (75 mm) and August (90 mm). Of total annual precipitation (table 1), near the base of the Jemez Mountains 36% falls in early-afternoon convective events during July and August (Bowen 1996). Monthly mean temperatures averaged throughout the watershed typically start the year as low as −4 °C in January and attain a peak of 17 °C in July. Averaging over the watershed, temperatures typically remain below freezing in December through February. However, the coldest (high elevation) locations typically remain below freezing from November through March. At the Quemazon SNOTEL site located on the eastern rim of Valles Caldera snow may begin to accumulate at high elevations as early as October. Snow water equivalent (SWE) stored in the mountain snowpack typically peaks in March and most snow melts before May. Median monthly streamflow from the perennial Jemez River peaks during spring snowmelt in April (13 mm) and May (10 mm) and drops to 1.5 mm by July.
Table 1. Watershed characteristics.
| Watershed | Area | Elevationa range | Precipitationb range | Mean annual precipitation | Mean annual streamflow | Maximum temperatureb range | Minimum temperature range |
|---|---|---|---|---|---|---|---|
| km2 | m | mm yr−1 | mm | mm | °C | °C | |
| Jemez | 1223 | 1720–3439 | 310–1040 | 630 | 54 | 9–21 | −5–3 |
| Mogollon | 191 | 1660–3239 | 490–1130 | 720 | 144 | 11–22 | −1–5 |
| Gila | 4807 | 1423–3329 | 380–1210 | 520 | 34 | 9–24 | −2–6 |
aElevation is derived from 30 m shuttle radar topographic mission digital elevation models. bPrecipitation (annual) and temperature (monthly mean) are estimated from PRISM.
Rocky Mountain Spruce-Fir Forest occurs at the highest elevations (2750–3430 m), dominated by Engelmann spruce (Picea engelmannii) and corkbark fir (Abies lasiocarpa). Partly below this ecosystem (2630–3110 m), lies moist Rocky Mountain Aspen Forest and Woodlands, dominated by quaking aspen (Populus tremuloides). Rocky Mountain Mixed Conifer Forest lies in a similar elevation range (2600–3050 m), dominated by Douglas-fir (Pseudotsuga menziesii), white fir (Abies concolor), blue spruce (Picea pungens), southwestern white pine (Pinus strobiformis), limber pine (Pinus flexilis), and ponderosa pine (Pinus ponderosa). Generally below aspen and mixed conifer forest ecosystems lies Rocky Mountain Ponderosa Pine Forest and Woodland (2450–2840 m), dominated by ponderosa pine. Typically, below these forest and woodland ecosystems lies shrubland. Rocky Mountain Montane Shrubland (2540–2870 m) is dominated by Gambel oak (Quercus gambelii) and New Mexico locust (Robinia neomexicana). Rocky Mountain Montane Riparian Shrubland occurs in the same elevation range along streams and is dominated by thinleaf alder (Alnus tenuifolia). Within a similar elevation range (2560–2870 m) Rocky Mountain Montane Grasslands dominated by upland bunch grasses occur (Muldavin et al 2006).
The Gila and Mogollon watersheds follow similar climatic and ecological patterns to the Jemez, but with higher annual precipitation and temperature (table 1), reflecting their more southerly positions closer to the gulfs of California and Mexico. Both the Gila and Mogollon watersheds are located within the Arizona-New Mexico mountains semi-desert-open woodland-coniferous forest-alpine meadow province. The Gila region is characterized by unusually high occurrences of lightning-ignited typically low-severity fires, with five lightning-caused fires per 100 ha annually (Rollins et al 1999, 2002). In this region hillslopes with northeastern topographic aspect tend to experience less water-limitation, greater productivity, greater continuity of fine fuels, and as a consequence higher fire frequency relative to other aspects (Rollins et al 2002). Historically, ponderosa pine and Douglas fir potential vegetation types have been most likely to burn relative to other land-cover types (Rollins et al 2002). (Potential vegetation types classify sites based on the climax community predicted for a site (Keane et al 1999).)
2.3. Post-wildfire vegetation recovery
To assess the duration over which transpiration would be reduced following a wildfire, we assessed recovery of enhanced vegetation index (EVI) over the six years following 53 non-intersecting wildfires that affected the Jemez, Mogollon, or Gila watersheds. Mesophyll tissue in healthy plant leaves strongly reflects infrared radiation whereas chlorophyll absorbs red light (Campbell and Wynne 2011). These characteristics facilitate remote sensing of plant canopy vigor and post-wildfire vegetation recovery. The EVI signal is related to the difference between reflectance in the near infrared and red regions of the electromagnetic spectrum, which is associated with the photosynthetic activity of healthy vegetation. We obtained all Landsat-derived cloud-free at-surface-reflectance-based (Masek et al 2006) EVI scenes from 1984 to 2014 that were not affected by scan-line corrector striping. To avoid issues related to overlapping fires we modified the MTBS wildfire perimeter polygons to remove any overlapping wildfire areas. With overlapping areas removed, we extracted EVI for each wildfire polygon from each EVI image using the arcpy site package to yield an annual time series of peak EVI before and after each wildfire. Altogether there were 53 wildfires for which EVI values were available during the year before the wildfire as well as the following six years. We determined the number of years post-wildfire during which transpiration was reduced using a repeated measures ANOVA followed by a post hoc test consisting of paired t-tests subjected to a Bonferroni adjustment.
2.4. Wildfire hydrologic assessment
Following selection of watersheds with natural flow regimes subjected to repeated catastrophic wildfires, regression analyses required additional information about the characteristics of fires in each watershed as they relate to watershed-scale hydrologic response. To this end we created two different predictors of wildfire impacts on hydrologic response for analysis of the monsoon season (June–September). The first predictor (INF3) interrogates the combined influence of increases in infiltration excess overland flow and decreases in transpiration. A meta-analysis by Vieira et al (2015) suggested that significant increases in infiltration-excess overland flow are limited to the three years following a wildfire, but do not differ significantly with burn severity. Thus the first predictor quantifies the total watershed area burned in the previous three years. The second predictor (TRANSP6) interrogates the influence of reduced post-wildfire transpiration on streamflow. This predictor quantifies the total watershed area burned in the previous six years, the duration over which transpiration is reduced as found in this study. At an annual time step, we followed a similar scheme to determine which water years experience hydrologic effects associated with wildfires, but include only wildfires that burned during the previous two (INF2) and five years (TRANSP5), respectively. This decrement of the period over which fire areas were integrated reflects the fact that large wildfires typically burn after spring runoff, which dominates annual streamflow in these watersheds.
2.5. Hydrologic time series derivation
Assessing the hydrologic relations associated with large wildfires first required isolating these effects by removing climatic influences on discharge. To this end we assessed the capacity of candidate predictor variables—precipitation, temperature, and SWE—to predict discharge metrics. We used Python functionality available through the arcpy site package at version 10.3 of ArcGIS to derive from existing monthly 4 km PRISM raster fields (Daly et al 1994) a time series of annual descriptive statistics for temperature and precipitation within each watershed by water year. PRISM is an algorithm specifically developed for spatial prediction of precipitation and temperature in mountainous areas subject to orographic precipitation and rain shadow effects that may not be appropriately captured in many conventional geostatistical approaches. To determine snowpack dynamics we calculated peak SWE at a representative SNOTEL (SNOw TELemetry; operated by the Natural Resources Conservation Service) station proximal to each watershed (figure 1). Finally, we obtained discharge measured by the USGS at the outlet of each watershed.
2.6. Statistical analysis
The objective of the statistical analysis was to assess controls on annual total streamflow (mm), peak flow (mm d−1), and low flow (mm d−1) at both the annual time scale and during the monsoon season (June–September). Peak and low flow were determined as the highest and lowest flows of each year, respectively, at a daily time step for the perennial Jemez and Gila watersheds. To minimize zero-flow data in the intermittent Mogollon watershed, low flows were represented by the median daily flow during the monsoon season and the first quartile of daily flow at the annual time scale. For each continuous variable we assessed the normality of its statistical distribution using the Shapiro–Wilk test (Royston 1982) and transformed non-normally distributed variables to normality using the 'Ladder of Powers' concept (Helsel and Hirsch 2002) with the understanding that hydrological variables are commonly positively skewed due to large, rare, extreme events. To enhance model realism we also allowed for interactions between terms such as between mean temperature and peak SWE or precipitation and wildfire. Following any necessary transformations, we implemented stepwise regression to select a subset of candidate variables that minimize the Akaike information criteria (Akaike 1974). (To promote significance of terms at 
we increased the multiple of the number of degrees of freedom used to calculate the penalty for adding additional terms to the model to 2.7.) To ensure that the regression models complied with applicable assumptions we assessed normality in the regression residuals via the Shapiro–Wilk normality test, the presence of outliers via the Bonferroni outlier test, constancy of error via the score test, autocorrelation of residuals via the Durbin–Watson test, and the acceptability of variance inflation factors. We required that climatic predictor variables have a known and expected physical relationship, either direct or inverse, with the response variables. Where wildfire terms emerged as significant we quantified the marginal contribution of the wildfire term (Gromping 2006). We used R 3.2.4 to implement all statistical analyses (R Core Team 2016).
3. Results
3.1. Post-wildfire vegetation recovery
Prior to wildfire, annual peak EVI averaged 
(standard error) in the 53 tested burn perimeters (figure 2). EVI decreased significantly in the year following fire (p < 0.001) to 0.22 ± 0.007. Thereafter, EVI gradually recovered through the fifth year following wildfire when it remained lower than pre-fire conditions (p < 0.001) with an annual peak value of 0.23 ± 0.005. Finally, six years after the initial fire recovery to 0.25 ± 0.006 was achieved as indicated by a spectral signature statistically indistinguishable (p = 1.0) from pre-fire conditions.
Figure 2. Evolution of enhanced vegetation index—an indicator of transpiration—from pre-wildfire conditions through post-wildfire recovery within 53 non-overlapping wildfire perimeters in the Jemez, Mogollon, and Gila watersheds of New Mexico. Error bars indicate a 95% confidence interval around the mean.
Download figure:
Standard image High-resolution image3.2. Jemez
From 1982 to 2014, seven wildfires burned 225 km2 (18%) of the Jemez watershed. The largest wildfires were the 2011 Las Conchas and 2013 Thompson Ridge wildfires, which burned 103 and 74 km2, respectively, within the Jemez watershed. Despite the large area burned during these two wildfires, the cumulative proportion of watershed area burned in the Jemez was the lowest of the three watersheds considered. After controlling for climatic and snowpack variability, wildfires were not correlated with any significant increases in total discharge or low flows during the monsoon season (table 2, figure 3). However, the interaction between precipitation and the proportion of the watershed burned in the previous three years (INF3) did aid in predicting peak flows (tables 2 and 4, figure 3). This significant increase in peak flows was associated with a positive precipitation anomaly in 2013 that coincided with a peak in the cumulative area burned in the previous three monsoon seasons.
Table 2. Statistically significant predictors and associated p-values of monsoon season total streamflow (q), peak flow (qpeak), and low flow (qmin).
| Site | Response | Yeara | P | Pmin | Pmax | SWEpeak | SWEMay | Tavg | Tmax | Tmin | INF3 | Pmax · INF3b | P · INF3b | R2 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Jemez | q | <0.001c | 0.02 | <0.001 | −0.001 | 84.0 | ||||||||
| qpeak | 0.142 | 0.063 | <0.001 | −0.009 | 0.011 | 66.8 | ||||||||
| qmin | 0.046 | 0.004 | 0.003 | −0.098 | 75.1 | |||||||||
| Mogollon | q | <0.001 | 0.022 | 64.6 | ||||||||||
| qpeak | <0.001 | 0.023 | 48.9 | |||||||||||
| qmin | <0.001 | 0.004 | −0.028 | −0.140 | 0.045 | 81.4 | ||||||||
| Gila | q | 0.006 | 0.005 | 0.300 | 0.047 | 79.8 | ||||||||
| qpeak | <0.001 | 0.032 | 69.3 | |||||||||||
| qmin | <0.001 | −0.002 | 53.3 |
aPredictors include total annual spatially-averaged precipitation (P), minimum monthly spatially averaged precipitation (Pmin), maximum monthly spatially averaged precipitation (Pmax), maximum spatially averaged temperature (Tmax), peak snow water equivalent (SWEpeak), May 1st snow water equivalent (SWEMay), minimum spatially averaged temperature (Tmin), and mean annual spatially average temperature (Tavg). bBullets denote interaction terms. cTabulated values indicate the p-value of each regression term included in each model; p-values associated with negative regression coefficients are assigned corresponding signs.
Figure 3. Monsoon season measured (bars) and modeled (points) precipitation, streamflow, peak flow, and low flow anomalies from 1982 to 2014 for the Jemez (left), Mogollon (center) and Gila (right) watersheds. The orange line in the top panel indicates the cumulative proportion of the watershed burned in the previous three years.
Download figure:
Standard image High-resolution imageAt the annual time step, significantly lower annual low flows from the Jemez watershed occurred in association with INF2 (table 3). These anomalously low flows occurred in late June (table 3; figure 4). Two of these low flows coincided with unusually low precipitation during June 2012 and 2014. The third low flow occurred while the Thompson Ridge fire actively burned in 2013.
Table 3. Statistically significant predictors of annual total streamflow (q), peak flow (qpeak), and low flow (qmin). Predictor abbreviations are given in table 2.
| Site | Response | Year | P | Pmin | Pmax | SWE | Tmax | Tmin |
|
|
TRANSP5 | INF2 | R2 (%) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Jemez | q | 0.005 | <0.001 | −0.034 | 88.0 | ||||||||
| qpeak | <0.001 | −0.009 | 82.6 | ||||||||||
| qmin | −0.078 | 0.001 | −0.005 | −0.053 | 76.9 | ||||||||
| Mogollon | q | 0.041 | <0.001 | 0.002 | <0.001 | 83.0 | |||||||
| qpeak | 0.021 | 0.003 | 32.3 | ||||||||||
| qmin | 0.009 | 0.071 | 53.5 | ||||||||||
| Gila | q | <0.001 | 0.009 | 0.006 | 81.8 | ||||||||
| qpeak | 0.808 | −0.019 | 0.005 | 0.054 | −0.627 | 0.052 | 0.034 | <0.001 | 80.4 | ||||
| qmin | <0.001 | 49.6 |
Figure 4. Annual measured (bars) and modeled (points) precipitation, streamflow, peak flow, and low flow anomalies from 1982 to 2014 for the Jemez (left), Mogollon (center) and Gila (right) watersheds. Blue and red bars refer to positive and negative anomalies respectively. The orange line in the top panel indicates the cumulative proportion of the watershed burned in the previous two years; x's indicate the proportion of the watershed burned in the previous five years.
Download figure:
Standard image High-resolution image
3.3. Mogollon
From 1982 to 2014, eight wildfires burned 352 km2 within the Mogollon watershed. (The total area of this watershed burned exceeds the watershed area because certain areas were burned more than once.) The largest wildfires were the 2012 Whitewater-Baldy complex, 1996 Lookout, and 2003 Dry Lake complex wildfires that burned 175, 56, and 46 km2, respectively, within the Mogollon watershed). Of the three watersheds considered, the Mogollon watershed experienced the greatest proportion of area burned cumulatively. After controlling for climatic and snowpack variability, total, peak, and low flows from the Mogollon watershed during the monsoon season all increased significantly in association with INF3 (tables 2 and 4; figure 3). Total, peak, and low flows each peaked in 2013 when a positive precipitation anomaly followed the Whitewater-Baldy complex fire, which burned 92% of the watershed in the previous year. The model predicting monsoon season flow from the Mogollon watershed explained 64.6% of the total streamflow variance, of which 7% was attributed to wildfire effects (table 4). The model suggests that of the total 950 mm of measured discharge from the Mogollon watershed from 1982 to 2014 during the monsoon season, 200 mm of streamflow may have discharged due to the influence of wildfire. Of this 200 mm of monsoonal discharge associated with wildfire impacts 120 mm discharged during 2013. Following the Whitewater-Baldy complex wildfire runoff coefficients were 48% and 36% during the 2013 and 2014 monsoon seasons, respectively. These fire-affected runoff coefficient are 12–16 times higher than mean runoff coefficient during monsoon seasons unaffected by wildfire (3%).
Table 4. Relative importance—proportion of variance explained—of wildfire term in regression model.
| Site | Response | Period | Importance (%) |
|---|---|---|---|
| Jemez | q | Monsoon | —a |
| Jemez | qpeak | Monsoon | 17.1 |
| Jemez | qlow | Monsoon | — |
| Mogollon | q | Monsoon | 10.9 |
| Mogollon | qpeak | Monsoon | 20.6 |
| Mogollon | qlow | Monsoon | 5.0 |
| Gila | q | Monsoon | 7.6 |
| Gila | qpeak | Monsoon | 6.9 |
| Gila | qlow | Monsoon | — |
| Mogollon | q | Annual | 21.7 |
| Mogollon | qpeak | Annual | 63.3 |
| Mogollon | qlow | Annual | — |
| Gila | q | Annual | 11.6 |
| Gila | qpeak | Annual | 9.4 |
| Gila | qlow | Annual | — |
| Jemez | q | Annual | — |
| Jemez | qpeak | Annual | — |
| Jemez | qlow | Annual | — |
aIndicates that wildfire term was not selected for best subset regression even when permitted.
At an annual timescale both total streamflow and peak flows increased significantly in association with INF2 (tables 3 and 4; figure 4). The model predicting streamflow explained 83% of the total variance of which 17% was associated with wildfire influence. Of the total annual discharge from the Mogollon watershed from 1982 to 2014 (4740 mm), 550 mm were associated with wildfire. Of this 170 and 100 mm discharged in 2013 and 2014, respectively. Annual runoff coefficients in 2013 and 2014 were 39% and 29%, respectively. The mean runoff coefficient in years unaffected by wildfire was 10%.
3.4. Gila
From 1982 to 2014 97 wildfires burned 3998 km2 on the Gila watershed. The largest wildfires were the 2012 Whitewater-Baldy complex, 2003 Dry Lake complex, and 2011 Miller complex wildfires, which burned 559, 410, and 341 km2, respectively within the Gila watershed. Of the three watersheds considered the greatest total cumulative area was burned within the Gila watershed. From 2003 to 2005 and in 2013 over 20% of the watershed area had been burned in the previous three years. After controlling for climatic and snowpack variability, total and peak monsoonal flows significantly increased in association with INF3 (table 2, figure 3). Record total and peak flows occurred in 2013 when a positive precipitation anomaly followed both the aforementioned Whitewater-Baldy and Miller complex fires, which together yielded INF3 of 20%. The monsoon season streamflow model predicted 79.8% of total streamflow variance of which 6% was attributed to wildfire (table 4). The model suggests that of the 290 mm of total monsoon season discharge from the Gila watershed from 1982 to 2014, 112 mm were associated with wildfires. Of this only 30 mm of wildfire-related flows occurred in 2013, reflecting a less episodic wildfire influence than on smaller watersheds, such as Mogollon.
At an annual timescale total discharge and peak flows increased significantly in association with TRANSP5 and INF2, respectively (tables 3 and 4; figure 4). The streamflow model explained 81.8% of total annual streamflow, of which 9% was related to wildfires. Of total annual discharge (1982–2014) of 1130 mm, 250 mm was associated with wildfire.
4. Discussion
The present study provides early evidence of certain amplified post-wildfire hydrologic responses at the large watershed scale that are consistent with increased regional water supply. Results from the Mogollon watershed in the present study as well as past work have extensively documented increased groundwater recharge and peak flows at small scales. Past work suggests that hydrologic effects of wildfires may be undetectable on large watersheds such as the Gila because at the scale of large watersheds large-scale climatic patterns dominate hydrologic response (Blöschl et al 2007). However, the present study demonstrates that in watersheds such as the Gila with high frequencies of large wildfires—associated with unusually high frequencies of lightning strikes—up to a third of the watershed may be affected by recent wildfires. In such large watersheds we demonstrate a direct relationship between wildfires and water yields.
Whereas at the scale of small watersheds we observed large short-lived streamflow increases following wildfires, at the large watershed scale we observed proportionally smaller increases in streamflow that are sustained over time as new sites are regularly burned while previously burned sites recover. The smaller increases in streamflow at the scale of the Gila watershed relative to the Mogollon watershed reflect a smaller proportion of the watershed burned. The proportion of the watershed burned interacts with precipitation anomalies, especially anomalies during the monsoon season, to determine the marginal increase in streamflow each year.
Past work in Mediterranean climates has suggested increased regional water yield as a result of wildfires (Kinoshita and Hogue 2015, Bart 2016). Bart (2016) suggested that post-wildfire streamflow responses remain elevated during the seven years following a wildfire and that the magnitude of this effect is related to the proportion of area burned and precipitation. This California work was limited to a Mediterranean climate unaffected by snow and in which precipitation is out of phase with evaporative demand. Thus, the present study extends this work to demonstrate increases in regional water yield associated with wildfires in a climate co-dominated by spring snowmelt and the North American monsoon during the summer, and in which precipitation and evaporative demand are in phase. Furthermore, the largest watershed examined by Bart (2016) was the 625 km2 Arroyo Seco (North) watershed, whereas the present study demonstrates increases in post-wildfire water yield at yet a larger spatial scale, that of the 4807 km2 Gila watershed.
The absence of a significant increase in water yields in the Jemez contrasts with the significant increase observed in the Gila and Mogollon watersheds. In fact this result is more consistent with our understanding of scale based on the literature. There exists an extensive body of hydrologic literature showing that hydrologic effects measured at the point, hillslope, and small watershed scales are often negligible, undetectable, or masked at larger spatial scales relevant to water resources management (Wilcox et al 2006, Wilcox and Huang 2010) because the importance of particular processes in controlling a hydrologic response is itself a function of scale (Dooge 1997, Blöschl 2006). Bosch and Hewlett (1982) suggest that vegetation treatments of less than 20% do not yield measurable changes in streamflow. This apparent 20% threshold could reflect some combination of limitations on measurement accuracy, increased transpiration by remaining plants, and increased evaporation. In any case, the 20% threshold suggested by Bosch and Hewlett (1982) reflected field to watershed scale measurements. The present study is consistent with a similar 20% threshold required to elicit hydrologic impacts of a disturbance at the large watershed scale.
Given that this study considered over 100 wildfires, measuring the small-scale hydrologic response of each would not likely be feasible. However, while the precise small-scale response of each is not known, the qualitative characteristics of post-wildfire impacts on the hydrologic cycle, including flooding have been documented following certain wildfires such as the Las Conchas wildfire (Orem and Pelletier 2015) and qualitatively documented in the Jemez in cases in which flumes were washed out following wildfires. Furthermore, inclusion of the Mogollon watershed, 92% of which was burned by the Whitewater-Baldy complex wildfire, provides a template for the expected hydrologic responses caused by wildfire. This template includes reduced transpiration (Dore et al 2010) and increased infiltration-excess overland flow following high intensity monsoonal storms (McLin et al 2001) resulting in increases in streamflow (tables 2–4, figures 2–4). The lack of small-scale measurements does mean that, for example, transmission losses between small and large scales are not known precisely and the precise streamflow contributions of areas burned at different severities are unknown.
Further research is needed in the form of a longer-term study that capitalizes on long-term wildfire atlases that span important changes in wildfire and land management paradigms. Such a study could test the hypothesis that increases in regional water yield are a peculiarity of the Anthropocene increase in wildfire frequency and were not detectable in earlier decades. Also, further research is necessary to better understand the 20% threshold suggested by Bosch and Hewlett, and supported by the present work. Specifically, numerical modeling is needed to better understand the extent to which the 20% threshold is related to measurement error versus increased ET in undisturbed areas. Finally, given that the present study considered burned area as area within wildfire perimeters, further research is necessary to address the role of burn severity in influencing changes in hydrologic response over time.
Acknowledgments
We wish to acknowledge funding from New Mexico EPSCoR WC-WAVE NSF award 1329470 and receipt of a 2015 New Mexico Geological Society grant-in aid. We also acknowledge valuable critique of a draft of this manuscript by Dr Thomas W Swetnam and an anonymous reviewer.