Does increased forest protection correspond to higher fire severity in frequent-fire forests of the western United States ?

There is a widespread view among land managers and others that the protected status of many forestlands in the western United States corresponds with higher fire severity levels due to historical restrictions on logging that contribute to greater amounts of biomass and fuel loading in less intensively managed areas, particularly after decades of fire suppression. This view has led to recent proposals—both administrative and legislative—to reduce or eliminate forest protections and increase some forms of logging based on the belief that restrictions on active management have increased fire severity. We investigated the relationship between protected status and fire severity using the Random Forests algorithm applied to 1500 fires affecting 9.5 million hectares between 1984 and 2014 in pine (Pinus ponderosa, Pinus jeffreyi) and mixedconifer forests of western United States, accounting for key topographic and climate variables. We found forests with higher levels of protection had lower severity values even though they are generally identified as having the highest overall levels of biomass and fuel loading. Our results suggest a need to reconsider current overly simplistic assumptions about the relationship between forest protection and fire severity in fire management and policy.


IntroductIon
It is a widely held assumption among federal land management agencies and others that a lack of active forest management of some federal forestlands-especially within relatively frequent-fire forest types such as ponderosa pine (Pinus ponderosa) and mixed conifers-is associated with higher levels of fire severity when wildland fires occur (USDA Forest Service 2004, 2016. This prevailing forest/fire management hypothesis assumes that forests with higher levels of protection, and therefore less logging, will burn more intensely due to higher fuel loads and forest density. Recommendations have been made to increase logging as fuel reduction and decrease forest protections before wildland fire can be more extensively reintroduced on the landscape after decades of fire suppression (USDA Forest Service 2004, 2016. The concern follows that, in the absence of such a shift in forest management, fires are burning too severely and may adversely affect forest resilience (North et al. 2009, Stephens et al. 2013, Hessburg 2016. Nearly every fire season, the United States Congress introduces forest management legislation based on this view and aimed at increasing mechanical fuel treatments via intensive logging and weakened forest protections.
However, the fundamental premise for this fire management strategy has not been rigorously

Methods
We used Gap Analysis Program (GAP) protection classes (USGS 2012), as described below, to determine whether areas with the most protection (i.e., GAP1 and GAP2) had a tendency to burn more severely than areas where intensive management is allowed (i.e., GAP3 and GAP4). We compared satellite-derived burn severity data for 1500 fires affecting 9.5 million hectares from years for which there were available data  among four different forest protection levels ( Fig. 1), accounting for variation in topography and climate. We analyzed fires within relatively frequent-fire forest types comprised of pine and mixed-conifer forests mainly because these are the predominant forest types at low to midelevations in the western United States, there is a large data set on fire occurrence, and they have been a major concern of land managers for some time due to decades of fire suppression. We defined geographic extent of forest types from the Biophysical Settings data set (BpS) (Rollins 2009;public communication, http://www.landfire.gov) that derived forest maps from satellite imagery and represents plant communities based on NatureServe's Ecological Systems classification. Baker (2015) noted that some previous work found ~65% classification accuracy of this system with regard to specific forest types and, accordingly, he analyzed groups of related forest types in order to improve accuracy. We followed his approach (see Appendix S1: Table S1). The categories selected from the Biophysical Settings map were ponderosa/Jeffrey pine and mixed-conifer forest types with relatively frequent-fire regimes (e.g., Swetnam and Baisan 1996, Taylor and Skinner 1998, Schoennagel et al. 2004, Stephens and Collins 2004, Sherriff et al. 2014), compared to other forest types with different fire regimes such as high-elevation forests and many coastal forests not studied herein. Forest types in our study totaled 29.2 million hectares ( Fig. 1; Appendix S1: Table S1). We used the BpS data to capture areas that were classified as forests before fire, because postfire vegetation maps can potentially show these same areas as temporarily changed to other vegetation types. We sampled our response and predictor variables on an evenly spaced 90 × 90 m grid within these forest types using ArcMap 10.3 (ESRI 2014). This created a data set of 5,580,435 independent observations from which we drew our random samples to create our models. The 90-m spacing was chosen because it was the smallest spacing of points that was computationally practical with which to operate.

Fires
The Monitoring Trends in Burn Severity project (MTBS, public communication, http://www. mtbs.gov) is a U.S. Department of Interior and Department of Agriculture-sponsored program that has compiled burn severity data from satellite imagery, which became available in 1984, for fires >405 ha, and was current up to 2014 (Eidenshink et al. 2007). The MTBS Web site allows bulk download of spatial products that include two closely related indices of burn severity: differenced normalized burn ratio (dNBR) (Key and Benson 2006) and relative differenced normalized burn ratio (RdNBR) (Miller and Thode 2007). Both indices are calculated from Landsat TM and ETM satellite imagery of reflected light from the earth's surface at infrared wavelengths from before and after fire to measure associated changes in vegetation cover and soil characteristics. We defined burn severity with the RdNBR index because it adjusts for prefire conditions at each pixel and provides a more consistent measure of burn severity than dNBR when studying broad geographic regions with many different vegetation types (Miller et al. 2009a, Norton et al. 2009). RdNBR values typically range from negative 500 to 1500 with values further away from zero representing greater change from prefire conditions. Negative values represent vegetation growth and positive values increasing levels of overstory vegetation mortality. The RdNBR values could be used to classify fires into discrete burn severity classes of low, medium, and high but this was not performed in our study, as we desired to have a continuous response variable in our models.
We intersected forest sampling points with fire perimeters downloaded from MTBS to determine fires that occurred in our analysis area, and censored fires with <100 sampling points (81 ha). The remaining points represented sampling locations from 2069 fires (Fig. 1). We extracted RdNBR values at each sampling point as our response variable as well as predictor variables that included topography, geography, climate, and GAP status. These sampling points were used to investigate the relationship between forest protection levels and burn severity (Appendix S1: Tables S2 and S3). We chose topographic and climatic variables based on previous studies that quantified the relationship between burn severity, topography, and climate (Dillon et al. 2011, Kane et al. 2015.

Topographic and climatic data
To account for the effects of topographic and climatic variability, we derived several topographic indices (Appendix S1: Table S2) from seamless elevation data (public communication, http://www. landfire.gov/topographic.php) downscaled to 90m 2 spatial resolution due to computational limits when intersecting sampling points. These indices capture categories of topography, including percentage slope, surface complexity, slope position, and several temperature and moisture metrics derived from aspect and slope position. We used the Geomorphometry and Gradient Metrics Toolbox version 2.0 (public communication, http:// evansmurphy.wix.com/evansspatial) to compute these metrics. We also computed several temperature and precipitation variables (Appendix S1: Table S3) by downloading climatic conditions for each month from 1984 to 2014 from the PRISM climate group (public communication, http://prism. oregonstate.edu). Climate grids record precipitation and minimum, mean, and maximum temperature at a 4-km grid scale created by interpolating data from over 10,000 weather stations. To determine the departure from average conditions, we subtracted each climate grid by its 30-yr mean monthly value. These "30-yr Normals" data sets were also downloaded from the PRISM Web site and reflected the mean values from the most recent full decades . We determined mean seasonal values with summer defined as the mean of July, August, and September of the year before a given fire; fall being the mean of October, November, and December of the previous year; winter the mean of January, February, and March of the current year of a given fire; and spring the mean of April, May, and June of the current year.

Protected area status and ecoregion classification
We used the Protected Areas Database of the United States (PAD-US; USGS 2012) to determine forest protection status, which is the U.S. official inventory of protected open space. The PAD-US includes all federal and most State conservation lands and classifies these areas with a GAP ranking code (see map at: http://gis1.usgs.gov/csas/ gap/viewer/padus/Map.aspx). The GAP status code (herein referred to interchangeably as GAP class or protection status) is a metric of management to conserve biodiversity with four relative categories. GAP1 is protected lands managed for biodiversity where disturbance events (e.g., fires) are generally allowed to proceed naturally. These lands include national parks, wilderness areas, and national wildlife refuges. GAP2 is protected lands managed for biodiversity where disturbance events are often suppressed. They include state parks and national monuments, as well as a small number of wilderness areas and national parks with different management from GAP1. GAP3 is lands managed for multiple uses and are subjected to logging. Most of these areas consist of non-wilderness USDA Forest Service and U.S. Department of Interior Bureau of Land Management lands as well as state trust lands. GAP4 is lands with no mandate for protection such as tribal, military, and private lands. GAP status is relevant to the intensity of both current and past managements.
We made one modification to GAP levels by converting Inventoried Roadless Areas (IRAs) from the 2001 Roadless Area Conservation Rule (S_USA.RoadlessArea_2001, public communication, http://data.fs.usda.gov/geodata/edw/datase ts.php) to GAP2 unless these areas already were defined as GAP1. We considered most IRAs as GAP2 given they are prone to policy changes and because they allow for certain limited types of logging (e.g., removal of predominately small trees for fuel reduction in some circumstances). However, we note that very little logging has occurred within IRAs since the Roadless Rule, although there occasionally have been proposals to log portions of some IRAs pre-and postfire, and fire suppression often occurs.
We modified level III ecoregions (U.S. Environmental Protection Agency (EPA) 2013) to create areas of similar climate and geography (Fig. 1). We did this by extracting ecoregions and combining adjacent provinces in our study region.

Random Forests analysis
We investigated the relationship between protection status and burn severity using the datamining algorithm Random Forests (RF) (Breiman 2001) with the "randomForestSRC" add-in package (Ishwaran and Kogalur 2016) in R (R Core Team 2013). This algorithm is an extension of classification and regression trees (CART) (Breiman et al. 1984) that recursively partitions observations into groups based on binary rule splits of the predictor variables. The main advantage of using RF in our study is that it can work with spatially autocorrelated data (Cutler et al. 2007). It can also model complex, nonlinear relationships among variables, makes no assumption of variable distributions (Kane et al. 2015), and produces accurate predictions without overfitting the available data (Breiman 2001).
Our independent observations were a random subset of our 5.5 million points, from which we drew three random samples of 25,000 points each. Each sample consisted of 500 fires randomly selected without replacement from the pool of 2069 fires. Fifty points were then randomly selected within each of the 500 fires. Our dependent variables were all continuous (Appendix S1: Tables S2 and S3) except for the main variable of interest, protected area status, which included the four GAP levels. The three observation samples were used to create three RF model runs, each consisting of 1000 regression trees. We conducted three RF model runs to assess whether our random samples of 25,000 points produced fairly consistent results.
The RF algorithm samples approximately 66% of the data to build the regression trees, and the remaining data are used for validation and to assess variable importance. We used this validation sample to determine the amount of variance explained and variable importance.
The algorithm also produces individual variable importance measures by calculating differences in prediction mean-square-error before and after randomly permuting each dependent variable's values. Variable importance is a measure of how much each variable contributes to the model's overall predicative accuracy.
Unlike linear models, RF does not produce regression coefficients to examine how a change in a predictor variable affects the response variable. The analogy to this in RF is the partial dependence plot which is a graphical depiction of how the response will change with a single predictor while averaging out the effects of the other predictors, such as the climatic and topographic variables (Cutler et al. 2007). We used this approach, in addition to using RF to determine overall variable importance as described above, in order to determine the effect of GAP status, in particular, on fire severity, while averaging out effects of climate and topography.

Mixed-effects analysis
We performed a linear mixed-effects analysis using the "nlme" add-on package in R (Pinheiro et al. 2015). We used a random intercept model and identified year of fire (n = 31) and ecoregion (n = 10) as random effects. Similar to our RF models, our independent observations were a random subset of our 5.5 million points but for these models we drew three random samples of 50,000 points each. Each sample consisted of 500 fires randomly selected without replacement, and within each of those fires, 100 points were randomly selected. Our dependent variables were the same used in our RF models, and we logtransformed the non-normal variables of slope, surface roughness, and topographic radiation aspect index. We removed dependent variables that were correlated with each other (Pearson's r > 0.5), retaining 21 of 45 candidate dependent variables, and centered these on their means. Model reduction was performed in a stepwise process using bidirectional elimination with Bayesian information criterion selection criterion.

Spatial autocorrelation analysis
Spatial autocorrelation (SA) is the measure of similarity between pairs of observations in relationship to the distance between them. Ecological variables are inherently autocorrelated because v www.esajournals.org BRADLEy ET AL. landscape attributes that are closer together are often more similar than those that are far apart.
We assessed the SA in the Pearson residuals with inspection of Moran's I autocorrelation index using the "APE" package add-in in R (Paradis et al. 2004) after removing points that shared the same x and y coordinates. Moran's I is an index that ranges from −1 to 1 with the sign of the values indicating strength and direction of SA. Values close to zero are considered to have a random spatial pattern. Our mixed-effects models all had a Moran's I values statistically different from 0 at the 95% confidence level (P < 0.001) so we included a spatial correlation structure in our model using the "nlme" package in R. Of Gaussian, exponential, linear, and spherical spatial correlation structures, we determined that the exponential structure produced the lowest Akaike's information criterion (AIC). Despite these additions, our second measurements still found relatively small, but significant, autocorrelation (Moran's I for model runs 1, 2, 3 = 0.10, 0.08, 0.10, all P < 0.001).

results
With regard to ranking of variables in the model runs, variable importance plots from the three RF model runs show that protection status was consistently ranked as one of the 10 most important of the 45 variables in explaining burn severity (Appendix S1: Table S4). The most important variable explaining burn severity was ecoregion for models 1 and 2 and maximum temperature from the previous fall for model 3.
With regard to the GAP status variable in particular, after averaging out the effects of climatic and topographic variables, the RF partial dependence plots show an increasing trend of fire severity with decreasing protection status (Fig. 2). Fires in GAP4 had mean RdNBR values greater than two standard errors higher than all other GAP levels. Fires in GAP3 had mean RdNBR values two standard errors higher than GAP1 in all model runs. GAP3 differences with GAP2 were less pronounced with only one model showing differences greater than two standard errors. Fires in GAP1 were consistently the least severe, being two standard errors less than GAP3 in all model runs and two standard errors less than GAP2 in two of three model runs.
Our mixed-effects models validated these findings with similar results (Fig. 3, Appendix S1: status levels 2 and 3 were not significantly different in the mixed-effects models. Although the level of autocorrelation was significant, it was small in our model (Moran's I ~0.1) and not enough to account for such a substantial difference in burn severity among protection classes.

Protected forests burn at lower severities
We found no evidence to support the prevailing forest/fire management hypothesis that higher levels of forest protections are associated with more severe fires based on the RF and linear mixed-effects modeling approaches. On the contrary, using over three decades of fire severity data from relatively frequent-fire pine and mixed-conifer forests throughout the western United States, we found support for the opposite conclusion-burn severity tended to be higher in areas with lower levels of protection status (more intense management), after accounting for topographic and climatic conditions in all three model runs. Thus, we rejected the prevailing forest management view that areas with higher protection levels burn most severely during wildfires.
Protection classes are relevant not only to recent or current forest management practices but also to past management. Millions of hectares of land have been protected from logging since the 1964 Wilderness Act and the 2001 Roadless Rule, but these areas are typically categorized as such due to a lack of historical road building and associated logging across patches >2000 ha, while GAP3 lands, for instance, such as National Forests lands under "multiple use management," have generally experienced some form of logging activity over the last 80 yr.
We expect that the effects of historic logging from nearly a century ago to gradually lessen over time, as succession and natural disturbance processes reestablish structural and compositional complexity, but it was beyond the scope of this study to attempt to assess the relative role of recent vs. historical logging. Similarly, industrial fire suppression programs that intensified in the 1940s influenced fire extent across forest protection classes. While more recent let-burn policies have been applied in GAP1 and GAP2 forests in some circumstances, evidence indicates that protected forests nevertheless remain in a substantial fire deficit, relative to the prefire suppression era , 2016, Parks et al. 2015. Thus, we believe it is unlikely that recent decisions to allow some backcountry fires to burn, largely unimpeded, account for much of the differences in fire severity among protection classes that we found, simply because such letburn policies have not been extensive enough to remedy the ongoing fire deficit.
While forests in different protection classes can vary in elevation, with protected forests often occupying higher elevations, our results indicate that protection class itself produced notable v www.esajournals.org BRADLEy ET AL. differences in fire severity after averaging out the effects of elevation and climate (see Fig. 2 and Results above). In our study, GAP1 forests were 284 m on average higher in elevation than GAP4 forests, while GAP1 forests experienced lower fire severity. This is the opposite of expectations if elevation was a key influence because higher elevation forests are associated with higher fire severity (see, e.g., Schoennagel et al. 2004, Sherriff et al. 2014. We note that we are not the first to determine that increased fire severity often occurs in forests with an active logging history (Countryman 1956, Odion et al. 2004).

Prevailing forest-fire management perspectives vs. alternative views
An extension of the prevailing forest/fire management hypothesis is that biomass and fuels increase with increasing time after fire (due to suppression), leading to such intense fires that the most long-unburned forests will experience predominantly severe fire behavior (e.g., see USDA Forest Service 2004, Agee and Skinner 2005, Spies et al. 2006, Miller et al. 2009b, Stephens et al. 2013, Lydersen et al. 2014, Dennison et al. 2014, Hessburg 2016). However, this was not the case for the most longunburned forests in two ecoregions in which this question has been previously investigated-the Sierra Nevada of California and the Klamath-Siskiyou of northern California and southwest Oregon. In these ecoregions, the most longunburned forests experienced mostly low/ moderate-severity fire (Odion et al. 2004, Odion and Hanson 2006, van Wagtendonk et al. 2012. Some of these researchers have hypothesized that as forests mature, the overstory canopy results in cooling shade that allows surface fuels to stay moister longer into fire season (Odion andHanson 2006, 2008). This effect may also lead to a reduction in pyrogenic native shrubs and other understory vegetation that can carry fire, due to insufficient sunlight reaching the understory (Odion et al. 2004(Odion et al. , 2010. Another fundamental assumption is that current fires are becoming too large and severe compared to recent historical time lines (Agee and Skinner 2005, Spies et al. 2006, Miller et al. 2009b, Stephens et al. 2013, Lydersen et al. 2014, Dennison et al. 2014, Hessburg 2016). However, others have shown that this is not the case for most western forest types. For instance, using the MTBS (www. mtbs.gov) data set, Picotte et al. (2016) found that most vegetation groups in the conterminous United States exhibited no detectable change in area burned or fire severity from 1984 to 2010. Similarly, Hanson et al. (2009) found no increase in rates of high-severity fire from 1984 to 2005 in dry forests within the range of the northern spotted owl (Strix occidentalis caurina) based on the MTBS data set. Using reference data and records of high-severity fire, Baker (2015) found no significant upward trends in fire severity from 1984 to 2012 across all dry western forest regions (25.5 million ha), nearly all of which instead were too low or were within the range of historical rates. Parks et al. (2015) modeled area burned as a function of climatic variables in western forests and non-forest types, documenting most forested areas had experienced a fire deficit (observed vs. expected) during 1984 to 2012 that was likely due to fire suppression.
Whether fires are increasing or not depends to a large extent on the baseline chosen for comparisons (i.e., shifting baseline perspective, Whitlock et al. 2015). For instance, using time lines predating the fire suppression era, researchers have documented no significant increases in high-severity fire for dry forests across the West Baker 2012a, Odion et al. 2014) or for specific regions (Williams and Baker 2012b, Sherriff et al. 2014, Tepley and Veblen 2015. Future trends, with climate change and increasing temperatures, may be less simple than previously believed, due to shifts in pyrogenic understory vegetation (Parks et al. 2016). This is more than just a matter of academic debate, as most forest management policies assume that fire, particularly high-severity fire, is increasing, is in excess of recent historical baselines, and needs to be reduced in size, intensity, and occurrence over large landscapes to prevent widespread ecosystem damages (policy examples include USDA Forest Service 2002, Healthy Forests Restoration Act 2003, USDA Forest Service 2009, HR 167: Wildfire Disaster Funding Act 2015. However, large fires (landscape scale or the so-called megafires) produce myriad ecosystem benefits underappreciated by most land managers and decision-makers (DellaSala and Hanson 2015a. High-severity fire patches, in particular, provide a pulse of "biological legacies" (e.g., snags, down logs, and native shrub patches) essential for complex early seral associates (e.g., many bird species) that link seral stages from new forest to old growth (Swanson et al. 2011, Donato et al. 2012, DellaSala et al. 2014, Hanson 2014, DellaSala and Hanson 2015a. Complex early seral forests are most often logged after fire, which, along with aggressive fire suppression, exacerbates their rarity and heightens their conservation importance (Swanson et al. 2011, DellaSala et al. 2014, Hanson 2014).

Limitations
One limitation of our study is that, due to the coarseness of the management intensity variables that we used (i.e., GAP status), we cannot rule out whether low intensities of management decreased the occurrence of high-severity fire in some circumstances. However, the relationship between forest density/fuel, mechanical fuel treatment, and fire severity is complex. For instance, thinning without subsequent prescribed fire has little effect on fire severity (see Kalies and yocum Kent 2016) and, in some cases, can increase fire severity (Raymond and Peterson 2005, Ager et al. 2007, Wimberly et al. 2009) and tree mortality (see, e.g., Stephens andMoghaddas 2005, Stephens 2009: Figure 6)-the effects depend on the improbable co-occurrence of reduced fuels (generally a short time line, within a decade or so) and wildfire activity (Rhodes and Baker 2008) and can be over-ridden by extreme fire weather (Bessie and Johnson 1995, Hély et al. 2001, Schoennagel et al. 2004, Lydersen et al. 2014. Empirical data from actual fires also indicate that postfire logging can increase fire severity in reburns (Thompson et al. 2007), despite removal of woody biomass (tree trunks) described by land managers as forest fuels (Peterson et al. 2015). While our study did not specifically test for these effects, such active forest management practices are common on GAP3 and GAP4 lands. Recognizing these limitations, researchers have stressed the need for managers to strive for coexistence with fire by prioritizing fuel reduction nearest homes and allowing more fires to occur unimpeded in the backcountry (Moritz 2014, Dunn and Bailey 2016, Moritz and Knowles 2016.
Follow-up research at finer scales is needed to determine management emphasis and history in relation to fire severity. However, we believe our findings are robust at the subcontinental and ecoregional scales.

conclusIons
In general, our findings-that forests with the highest levels of protection from logging tend to burn least severely-suggest a need for managers and policymakers to rethink current forest and fire management direction, particularly proposals that seek to weaken forest protections or suspend environmental laws ostensibly to facilitate a more extensive and industrial forest-fire management regime. Such approaches would likely achieve the opposite of their intended consequences and would degrade complex early seral forests . We suggest that the results of our study counsel in favor of increased protection for federal forestlands without the concern that this may lead to more severe fires.
Allowing wildfires to burn under safe conditions is an effective restoration tool for achieving landscape heterogeneity and biodiversity conservation objectives in regions where high levels of biodiversity are associated with mixed-intensity fires (i.e., "pyrodiversity begets biodiversity," see DellaSala and Hanson 2015b). Managers concerned about fires can close and decommission roads that contribute to human-caused fire ignitions and treat fire-prone tree plantations where fires have been shown to burn uncharacteristically severe (Odion et al. 2004). Prioritizing fuel treatments to flammable vegetation adjacent to homes along with specific measures that reduce fire risks to home structures are precautionary steps for allowing more fires to proceed safely in the backcountry (Moritz 2014, Moritz and Knowles 2016. Managing for wildfire benefits as we suggest is also consistent with recent national forest policies such as 2012 National Forest Management Act planning rule that emphasizes maintaining and restoring ecological integrity across the national forest system and because complex early forests can only be produced by natural disturbance events not mimicked by mechanical fuel reduction or clear-cut logging (Swanson et al. 2011, DellaSala et al. 2014 wishing to maintain biodiversity in fire-adapted forests should appropriately weigh the benefits of wildfires against the ecological costs of mechanical fuel reduction and fire suppression (Ingalsbee and Raja 2015) and should consider expansion of protected forest areas as a means of maintaining natural ecosystem processes like wildland fire.