Background

In the southwestern United States, post-fire vegetation recovery is increasingly variable in forest burned at high-severity. Many factors, including temperature, drought, and erosion, can reduce post-fire vegetation recovery rates. Here, we examined how post-fire precipitation variability, topography, and soils influenced post-fire vegetation recovery in the southwestern United States as measured by greenness. We modeled relationships between post-fire vegetation and these predictors using Random Forest and examined changes in post-fire normalized burn ratio across fires in Arizona and New Mexico. We incorporated growing season climate to determine if year-of-fire effects were persistent during the subsequent five years or if temperature, water deficit, and precipitation in the years following fire were more influential for vegetation greenness.

Results

We found reductions in post-fire greenness in areas burned at high-severity when heavy and intense precipitation fell on more erodible soils immediately post-fire. In highly erodible scenarios, when accounting for growing season climate, coefficient of variation for year-of-fire precipitation, total precipitation, and soil erodibility decreased greenness in the fifth year. While the effects of year-of-fire factors related to erosion were significant, they were small, and the variability explained by growing season vapor pressure deficit and growing season precipitation were significantly greater.

Conclusions

Our results suggest that while the factors that contribute to post-fire erosion and its effects on vegetation recovery are important, at a regional scale, the majority of the variability in post-fire greenness in high-severity burned areas in southwestern forests is due to climatic drivers such as growing season precipitation and vapor pressure deficit. Given the scale of area burned at high-severity, the likelihood that high-severity burned area will continue to increase, and the potential for more post-fire erosion that can result in different vegetation trajectories, quantifying how these factors alter the trajectory of greenness and what that means in terms of ecosystem development is central to understanding how different ecosystem types will be distributed across these landscapes with additional climate change.

We used the LANDFIRE Existing Vegetation Type layer (2016) to select ponderosa pine, mixed-conifer, and sub-alpine forests, then selected burned areas using Monitoring Trends in Burn Severity (MTBS) fire perimeters for the period 1985-2017. The MTBS program maps fire severity for fires with a minimum burned area of 404ha (Eidenshink et al. 2007). We only included fires through 2017 to ensure that we had at least five years of post-fire data for all fires in the dataset.

Within each fire perimeter, we aggregated gridded climate data, remotely sensed vegetation metrics, topography, and soil metrics. Unless otherwise stated, the following descriptions refer to pixelwise calculations. We calculated Normalized Burn Ratio (NBR) for each fire year -1 to fire year +5 (Eq. 1), which describes both the amount and greenness of vegetation in a pixel. Compared to other metrics such as Normalized Differenced Vegetation Index, NBR is more sensitive to post-fire recovery than other indices (Pickell et al. 2015). The study period encompassed multiple Landsat missions (5 TM, 7 ETM+, and 8 OLI), with 30m resolution and 16 day return intervals. We imported perimeters accessed from MTBS into Google Earth Engine (GEE), then for each perimeter we excluded pixels identified by the USGS quality band as clouded, water, etc. and applied coefficients and offsets from Roy et al. (2016) to Landsat 5 TM and 7 ETM+ products to allow for continuous data between Landsat 5 TM, 7 ETM+, and 8 OLI imagery (Gorelick et al. 2017). For the non-masked pixels, we calculated the normalized burn ratio (NBR, eq. 1) as:

Eq. 1 NBR = SWIR – RED / SWIR + RED

where SWIR is the short-wave infrared band, and red is the red band. We selected scenes between June 20 and September 20 to minimize snow-cover and capture peak summer greenness in southwestern United States (Notaro et al. 2010). We then used a GEE implementation of the Landtrendr algorithm to segment and fit pixelwise curves to time series of imagery to summarize each year’s peak greenness of vegetation with a single image of each pixel’s maximum NBR (Kennedy et al. 2018). Finally, we excluded pixels in which the pre-fire NBR was less than 50, a threshold that corresponds to hard surfaces and rock. 

We selected year-of-fire effects based on the potential for substrate loss on post-fire slopes and used variables related to those in the Universal Soil Loss Equation (USLE), which predicts long term average soil loss based on rainfall, topography, soil characteristics, and vegetation characteristics in agricultural settings (Wischmeier and Smith, 1978). Here, we used slope length factor, soil erodibility, total precipitation, and coefficient of variation of precipitation between the fire date and the end of the fire year, which are similar to the topographic, soil conditions, and rainfall used in the USLE to predict long term soil loss. These variables are strongly predictive of first year hillslope erosion in mountainous environments, though we substituted soil erodibility factor for percent rock to accommodate a wider array of soil textures and makeups (Miller et al. 2011). 

We accessed gridded, daily precipitation data (resolution: 1000m) from the Daymet V4 repository hosted on GEE (Thornton et al. 2020). We selected precipitation between the fire date and December 31st of the year of fire, then calculated the coefficient of variation and total precipitation that occurred in that period. We used the RSAGA tool ls_factor to calculate landscape slope factor from the USGS 3DEP 10m Digital Elevation Model, in which lower values correspond to shallower, shorter slopes and higher values correspond to steeper, longer slopes (U.S. Geological Survey). We accessed soil erosion likelihood from SSURGO and used ArcMap 10.8.1 to calculate the k-factor erosion potential of soils, in which values > 0.4 correspond to more erodible soils (Soil Survey Staff).

We accessed monthly Terraclimate data to use for the growing season climate variables (Table 1, 4000m) (Abatzoglou et al. 2018). We calculated the growing season (April through October) mean and sum of the climate variables for each post-fire year. We thinned the full dataset to 20% to reduce the likelihood of spatial autocorrelation and retained the remaining 80% for final model fit assessment. 

Methods for processing the data: 

To determine whether the year-of-fire effects persisted through the five years post-fire, we modeled relationships between post-fire vegetation and predictors with Random Forest using the Ranger package (Wright & Ziegler 2017). The Random Forest algorithm can capture complex, nonlinear relationships and interactions between response and predictor variables, without making a priori assumptions about these relationships. 

For each of the five post-fire years, we built five-fold cross-validated models with all predictors using the Ranger implementation of Random Forest in R with 500 trees and the square root of the number of variables per node split (R Core Team 2023). To determine the optimal number of variables, we used recursive variable elimination in a ten-fold cross-validation, repeated five times using the ‘rfe’ function from the ‘caret’ package and used Root Mean Squared Error (RMSE) as our performance metric, with sets of 1, 3, 6, and 9 variables. We refit a final model based on the optimal number of variables from the algorithm and calculated the permutation feature importance. We assessed final model fit by calculating the root mean squared error (RMSE) of model predictions and the coefficient of determination (R2) of a linear regression between observations and predictions of the 80% of data withheld from model development.

To examine the underlying relationships between variables and outcomes, we used tools from the ‘iml’ package to make global and local inferences (Molnar et al. 2018). and the strength of interactions between all variables to determine 1) whether variables did not meaningfully contribute to the models and 2) whether year-of-fire variables were influential when interacting with growing season variables, indicating their effects over longer periods. The interaction strength is the excess variance of the 2-dimensional partial dependence function greater than the sum of the two 1-dimensional partial dependence functions, measured by Friedman's H-statistic. 

To examine how a particular variable affected outcomes over its entire range, we calculated Accumulated Local Effects (ALE) for all individual variables and for pairwise combinations of year-of-fire variables to show the relative change in the prediction caused by a change in the variable, setting grid intervals to 20 units to extract clear trends. Finally, we calculated Shapley Additive exPlanations (SHAP) values for a likely highly erodible scenario using the year-five models to decompose how changes in year-of-fire variables can affect the baseline NBR value while still accounting for any interacting effects with growing season climate. (Table 2, Strumbeli & Kononenko 2014). The SHAP values decompose predictions into the contributions from each variable, showing how the difference from a baseline value is achieved. We used a scenario from a dataset consisting of mean climate, fires occurring in June, and 100 samples from the top 10th percentile of year-of-fire variables to decompose how year-of-fire variables can affect the baseline NBR value while still accounting for any interacting effects with growing season climate in a highly erodible post-fire scenario.