Climate drives inter-annual variability in probability of high severity fire occurrence in the western United States

Our spatial domain is a 1/8th degree latitude/longitude modeling grid (∼12 km resolution) west of –102.0625 longitude encompassing eleven states: WA, OR, CA, ID, UT, NV, MT, WY, CO, NM, AZ. The spatial resolution of data varies from 30 m (e.g. observed high severity burned area, vegetation) to ∼12 km (e.g. climate, topography). For higher resolution variables, we calculated fractional area of each variable within each modeling pixel.

The temporal domain was limited by Landsat image-derived burn severity data (1984–2014). Using land surface and hydroclimate predictors and monthly fire discovery events, we model gridded high severity fire occurrence at a monthly time step. Individual fire records contain the discovery date, which often approximates the ignition date and conditions likely conducive to rapid fire spread (Westerling et al 2006).

The California Energy Commission specified 1961–1990 as a reference period for the four California state climate and vulnerability assessments to date, to reflect average climate conditions prior to significant climate change impacts. We retain the 1961–1990 normal period for consistency and comparability of results with prior work by the authors.

We acquired burn severity data from the Monitoring Trends in Burn Severity database (MTBS Data Access 2009, www.mtbs.gov; accessed 12/2008 and 9/2016), resulting in a dataset with a total of 6808 fires (4493 through 2006) coded by discovery date (month, year), using thematic burn severity images comparable across space and time. We used ESRI Arc Macro Language (ESRI 1999) to intersect burn severity data with a 1/8th degree grid, assigning each fire to the grid cell where a majority of area burned. We then calculated fractional area burned in six severity classes (unburned to low severity, low severity, moderate severity, high severity, increased greenness, or unclassified) by voxel (latitude, longitude, year × month) (Eidenshink et al 2007).

Potential predictors of high severity fire included topographic and ecosystem characteristics. Topographic data (minimum, mean, maximum and standard deviation of elevation, slope, aspect) are 1/8th degree products derived from the GTOPO30 global 30 Arc Second (1 km) Elevation Data Set accessed online from the North American Land Data Assimilation System (LDAS) (http://ldas.gsfc.nasa.gov, Mitchell et al 2004).

We aggregated LANDFIRE existing vegetation type (EVT) and fire regime condition class (FRCC) variables to the modeling grid (Keane and Rollins 2007, www.landfire.gov). We used a reclassification of the EVT layer, extracting only the fractional area of forest types characterized by stand replacing fire regimes (Westerling et al 2011a). FRCC is a metric widely used to prioritize fuel treatments and characterizes ecosystem departure, in terms of vegetation structure and composition, from historical conditions (Hann 2004, Laverty and Williams 2000). We extracted the fractional area of each FRCC class: FRCC1 (departure <33%), FRCC2 (≥33% departure <66%), and FRCC3 (departure ≥67%) (Holsinger et al 2006, Keane and Rollins 2007). For this study, FRCC is fixed in time, reflecting year 2000 observed conditions.

We obtained hydrologic variables simulated with the Variable Infiltration Capacity (VIC) model and gridded climate data used to force VIC from the University of Washington National Hydrologic Prediction System (NHPS) (www.hydro.washington.edu/forecast/westwide/) (Liang et al 1994, Wood and Lettenmaier 2006). VIC calculates daily surface water and energy balances, estimating evaporation from vegetation canopy, bare soil surface, and transpiration by vegetation classes in each grid cell.

Climate data include temperature (Tmax, Tmin, Tave) and precipitation (PPT), while moisture deficit (MD), antecedent moisture deficit derivatives (e.g. 6 month prior moisture deficit), relative humidity (Rh), soil moisture, and snow water equivalent (SWQ) were derived from VIC and the Penman-Monteith equations (Penman 1948, Monteith 1965) on a monthly time step from 1915 present (Westerling et al 2011a, Westerling et al 2011b).

Stephenson (1998) showed that long term average MD and actual evapotranspiration (AET) are biologically meaningful drivers of the spatial distribution of vegetation types over multiple spatial scales. We used these and related variables as proxies for spatial variability in ecosystem and disturbance regime sensitivity to climate, including 30 years (1961–1990) means and standard deviations for Tave, PPT, cumulative MD, and AET. Standard deviations characterize inter-annual climate variability; for instance, higher precipitation standard deviation indicates locations with more dynamic precipitation regimes. We also created a thin plate spline interacting 1961–1990 average MD and AET to indicate biophysical site conditions for plant growth: different forest types fall along the gradient of MD and AET (Hastie et al 2001, Stephenson 1998).

We employed a multi-step process in modeling high severity fire occurrence. Previous research using parts of this dataset and similar questions used both logistic regression and classification and regression tree (CART) methods, individually and together (Collins et al 2007, Collins et al 2009, Dillon et al 2011, Preisler and Westerling 2007, Westerling and Bryant 2008, Westerling et al 2011a, Westerling et al 2011b). The number of potential predictor variables available for this analysis was 85. In order to limit the number of variables in building a predictive model, we first used a CART (Random Forest package in R; Liaw and Wiener 2002) model to identify the 20 most important variables in predicting fractional high severity. These provided the potential predictor set for developing conditional logistic regression models.

2.5.1. Logistic regression

We developed two conditional logistic regression models to predict high severity fire occurrence. Because MTBS contains only fires >400 ha, there exists an implicit a priori condition for our models: the occurrence of a >400 ha fire. Given this, we set a threshold for the presence of high severity equal to the median value of high severity fraction in the MTBS data, 0.042. Then, given presence of high severity fire, we set a threshold equal to the upper quartile cutoff, 0.1732. We define fires with high severity fraction above this threshold as high severity fires.

To model the probability of high severity fire presence, we use the logged odds, or logit:

$$\begin{equation} \begin{array}{*{20}c} {{\rm Model}\;{\rm Pa}:} \hfill & {{\rm Logit}\;(P|f > 0.042)} \hfill \\ {} \hfill & {\ln (P_{{\rm Pa}i} /(1 - P_{{\rm Pa}i} )) = \Sigma (b_0 + b_j X_{{\rm ij}} )} \hfill \\ \end{array} \end{equation}$$

where X_j is the set of predictor variables best fit to Model Pa, P_Pai is the probability of high severity fire presence, defined as the fraction f of high severity fire >0.042 for a given month and grid cell indexed by i.

Similarly, the model for occurrence of high severity fraction >0.1732 is:

$$\begin{equation} \begin{array}{*{20}c} {{\rm Model}\;{\rm Hi}:} \hfill & {{\rm Logit}\;(P|f > 0.1732|f > 0.042)} \hfill \\ {} \hfill & {\ln (P_{{\rm Hi}i} /(1 - P_{{\rm Hi}i} )) = \Sigma (b_0 + b_j X_{{\rm ik}} )} \hfill \\ \end{array} \end{equation}$$

where X_k denotes predictors best fit to Model Hi, where P_Hii is as P_Pai above. The probability of high severity fire occurrence (conditional on >400 total ha burned) for any month and location is the product of these two model probabilities: P_Pai^* P_Hii.

The thresholds for high severity fire presence and occurrence are necessarily arbitrary. Our goal was to be as objective as possible, while defining meaningful thresholds for relatively rare events, and to demonstrate the predictability of high severity fires. Our models specifically address the question:

Given that a fire burns >400 ha, and given that high severity fire is present, what is the probability that this is a high severity fire (i.e. in the upper quartile of high severity burned area)?

We use the Aikake Information Criterion (AIC) to evaluate model performance (Aikake 1974, Aikake 1981). The best model optimizes model fit (minimizes AIC) while penalizing excess predictive parameters. We performed leave-one-out cross-validation of the best model. We tested variable importance by calculating probabilities and testing model performance after removing variable groups.

2.5.2. Mapping probability of high severity fire occurrence

We tested trends for WUS, each state, and each month. We found no significant trend in WUS high severity fire occurrence over 1984–2014, except for Colorado (table S1 available at stacks.iop.org/ERL/12/065003/mmedia). While some studies have shown increasing fire season length, we saw no significant increase in high severity fire occurrence by month, May through October (figure S1). We found no correlation between fraction of high severity fire and total fire size, meaning increasing large fires does not necessarily increase fractional high severity fire area.

Seasonal occurrence of high severity fire coincided with WUS burned area (figure S1). However, the fine-grained distribution of high severity fire is quite variable. California and Idaho had the largest number of large fires and high severity fire occurrence, but many fires had no presence of high severity fire. Both Montana and Wyoming experienced fewer large fires than other states, often with no high severity fire present (figure S2).

The best predictive models for high severity fire presence included all variable groups—climate, topography, and vegetation (table 1). Pixel mean and maximum elevations are the only topographic variables in the best model. Climate normals include standard deviation of 1961–1990 cumulative annual water-year precipitation and moisture deficits and the thin plate spline interacting long term average MD and AET. Fire-year climate variables included average spring temperature (SPRT), average monthly temperature and normalized monthly moisture deficit at fire discovery (MD0), and MD the previous November (MD2).

Our best models included two vegetation variables. Fractional area of forest types with stand replacing fire regimes was important for Models Pa and Hi, while fractional area of FRCC3 was significant only in Model Pa.

Removing fire-year climate, our models do not predict extremes in inter-annual variability of high severity fire occurrence. Removing vegetation variables does not significantly alter predicted vs. observed high severity fire occurrences (high severity fire occurrence in vegetation-only models still varies with time because it is conditional on observed fires). Fixed site-specific variables describe a constant spatial distribution for probability of high severity fire given a large fire, identifying location and average recurrence rates for fire regimes where high severity fire occurs. While all model iterations perform well with regard to prediction vs. observation for locations with fires, the AIC is lowest for the full model (table 2, figure 1).

Zoom In Zoom Out Reset image size

To illustrate effects of variable groups, we mapped differences for two years with few (1991 [N = 16] and 1997 [N = 10]) and two years with many (1996 [N = 75] and 2000 [N = 88]) high severity fires (figure 2). Removing fire-year climate produces probability maps with a negligible shift in WUS average probabilities for high severity fire occurrence, but with distinct regional differences (figure 2). For the low year 1991, removing fire-year climate increases the maximum predicted probability from 79% to 90%. For 1996, a high occurrence year, removing fire-year climate increases the maximum from 76% to 90%. Probability of occurrence decreases over California; this pattern is opposite the Northern Rocky Mountain and North Cascade regions, where probabilities increase when we remove fire-year climate.

Zoom In Zoom Out Reset image size

Our results highlight the importance of inter-annually varying climate for high severity fire occurrence. Fire-year climate significantly modified conditional probabilities of high severity fire occurrence across the WUS, though with stronger effects in Sierra Nevada and Southwest forests with mixed-severity fire regimes than in Northwest and Northern Rockies forests with high severity fire regimes.

Others have shown a lengthening fire season and increasing large fire occurrence and burned area due to climate change (Westerling et al 2006, Dennison et al 2014, Jolly et al 2015, Westerling 2016). While conditions under which high severity fires in our record burned were warmer and drier than the long term average (figure S3), the lack of observed trends in fire severity occurrence (table S1) may be at least in part due to the short record, which begins right at the time of the largest temperature-driven increase in WUS forest burned area (Westerling et al 2006).

Similarly, the record for fire severity begins >70 years after the start of fire suppression. Ecosystems with short fire return intervals (7–10 years) may have missed up to 10 fire cycles by the time the MTBS severity record begins, while those with longer fire return intervals may have missed few or none, so we would only expect to see changes in severity due to fire suppression and fuel buildup on some of the landscape. Observed patterns of high severity fire occurrence may be influenced by changing climate and management, but without more observations, we can not quantify significant trends in high severity fire occurrence.

Our best models included both climate and land surface characteristics, e.g. topography and vegetation, to predict fire severity. The impact of removing vegetation or fire-year climate predictors is subtle. Considering first removing temporally fixed predictors such as fractional area in FRCC3 and in forest types with high severity fire regimes, inter-annual variability of model predictions is not significantly affected. Spatially, effects are most pronounced in mountainous areas, as one of the variables is forest fraction with high severity fire regimes (figure 2). We see a general decrease in probability of occurrence in the Sierra Nevada, western Nevada, and Northern Rocky Mountains. Both the Colorado Rocky Mountains and the Northern Cascades in Washington generally show a higher probability of high severity fire occurrence when we remove vegetation variables.

When we remove fire-year climate variables, conditional probabilities of high severity fire generally decrease in California and Nevada, but increase in the Pacific Northwest and Northern Rocky Mountain regions (figure 2). Many of the ecosystems in the Pacific Northwest and Northern Rocky Mountains are dominated by cool moist forests with historically infrequent stand replacing fire (Agee et al 1977, Agee 1993, Cansler and McKenzie 2014, Schoennagel et al 2004). These forests have abundant fuels, but are rarely hot and dry enough to burn. Increased probabilities in these regions after removing fire-year climate reflect climate's importance in determining whether these forests will burn and supports findings that fire-year climate controls severity occurrence in this forest type (Cansler and Mckenzie 2014, Schoennagel et al 2004). Previous studies found fire-year climate less important than topography and vegetation as a predictor of high severity fire occurrence in the North Cascades and Northern Rocky Mountains (Cansler and McKenzie 2014, Birch et al 2015, Dillon et al 2011). Our results indicate that in these systems with a propensity for high severity fire, fire-year climate still modulates high severity fire occurrence.

Probability maps show fire-year climate both amplifies and moderates probabilities of high severity fire based on biophysical site characteristics. Fire-year climate anomalies are highest for many high occurrence fire years: Tave, SPRT, MD0, MD2 (figure S3). November MD (MD2) is an indicator for moisture stress early in the water year; Van Mantgem et al (2013) found that high pre-fire climatic water deficit is related to increased post-fire tree mortality. These variables all indicate high severity occurrence years are drier and hotter than average, conditions that increase fuel flammability and enhance conditions for fire spread.

While climate variability significantly affected conditional probabilities of high severity fire, those probabilities are conditional on fires >400 ha. Though the predictive maps for years like 1991 and 1996 look very similar in terms of conditional probability of high severity fire occurrence, the actual fire record is quite different. The number of fires in the MTBS record for 1991 was 75 total (16 high severity) v. 272 total (75 high severity) fires for 1996. While ignitions resulting in large fires were few, modeled probabilities indicate that had more large fires occurred, high severity fire would have been likely in 1991.

The inclusion of inter-annually varying climate is critical for capturing high probability episodes in areas where fire severity is highly variable, especially California and the Southwest, where the coefficient of variation (thus influence of climate) is higher (figure 3). Our models imply substantial fraction high severity fire in mixed severity fire regimes requires more extreme climatic conditions. Conversely, in the Northern Rocky Mountain and Pacific Northwest, regions that are dominated by cool moist forests with high severity fire regimes, the coefficient of variation (and effect of climate) is lower.

Zoom In Zoom Out Reset image size

Our results are similar to Dillon et al (2011) regarding the importance of topography, but differ with respect to the importance of fire-year climate and the ability to predict extremes in high severity fire occurrence. This is likely due to differences in methodology and data. In our study each fire is classified as a high severity fire or not. Dillon et al (2011) selected a random subset of individual 30 m fire pixels classified to high severity or not. This means that they could select multiple pixels from a single fire with different classification, whereas their climate data would be identical for pixels from a given fire.

Dillon et al (2011) tested different independent climate variables, which could also lead to the differing conclusions. For each fire, we used hydroclimate variables for the voxel where the majority of area burned. Dillon et al (2011) interpolated monthly temperature and precipitation for the central latitude, longitude and mean elevation of each fire. Their soil moistures were simulated with VIC, as were our hydroclimate data, but at coarser scale. We also test a larger set of hydroclimate variables. Comparing approaches, we expect theirs would have more power to describe topographic controls on fire severity, whereas the present study may be better suited to demonstrating climatic controls.