DOI: 10.5061/dryad.mw6m9063p
title: Climate is more influential to vegetation green-up than factors that contribute to erosion following high-severity wildfire

description: 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.
description: 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
 &gt; 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. 
 <em><strong>Methods for processing the
 data: </strong></em> 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.
description: This readme file was generated on 2023-12-04 by Joseph L. Crockett #
 GENERAL INFORMATION ## Title of Dataset: Climate is more influential to
 vegetation green-up than factors that contribute to erosion following
 high-severity wildfire ## Author/Principal Investigator Information Name:
 Joseph L. Crockett Institution: University of New Mexico Address:
 Department of Biology, University of New Mexico, MSC03 2020, Albuquerque,
 NM 87131-0001, United States Email:
 [jcrockett@unm.edu](mailto:jcrockett@unm.edu) ## Associate or
 Co-investigator Contact Information Name: Matthew D. Hurteau
 Institution:University of New Mexico Address: Department of Biology,
 University of New Mexico, MSC03 2020, Albuquerque, NM 87131-0001, United
 States Email: [mhurteau@unm.edu](mailto:mhurteau@unm.edu) ## Date of data
 collection: 2017-2023 ## Geographic location of data collection: New
 Mexico, Arizona ## Information about funding sources that supported the
 collection of the data: United States Department of Agriculture National
 Institute of Food and Agriculture Interagency Carbon Cycle Science program
 (Grant No. 2017-67004-26486/ project accession no. 1012226 ); United
 States Department of Agriculture, National Institute of Food and
 Agriculture, Agriculture and Food Research Initiative program (Grant No.
 2021-67034-35106/project accession no. 1026366; Joint Fire Science Program
 (Project JFSP 16-1-05-8, JFSP 20-1-01-9); New Mexico Space Grant
 Consortium. # SHARING/ACCESS INFORMATION # DATA & FILE OVERVIEW ##
 File List: **Data** `train_nbr15.rdata` **Code/Software**
 `Model_notebook_for_dryad.Rmd` `landtrendr_dryad_gee.txt` `ranger_func.R`
 # METHODOLOGICAL INFORMATION ## Instrument- or software-specific
 information needed to interpret the data: R version 4.1.0 (2021-05-18)
 'Camp Pontanezen' Packages used:
 rmarkdown,broom,caret,dplyr,fastshap,ggpubr,iml,knitr,lme4,pdp,ranger,reshape2,scales,shapviz,stringr,terra,tidyverse,ggthemes,grid,RSAGA In addition, RStudio is required to access `Model_notebook_for_dryad.Rmd` Google Earth Engine Code Editor is required to run the code in` landtrendr_dryad_gee.txt` # DATA-SPECIFIC INFORMATION FOR: train\_nbr15.rdata train_nbr15.rdata contains five datasets (train_nbr1, train_nbr2, train_nbr3, train_nbr4, train_nbr5) with the same structure, as follows: 1382557 obs. of 14 variables **Variable List:** nbr_X_year: Normalized Burn Ratio for year X (unitless) postfire_precipitation_total: Total daily precipitation (mm), fire date to end of fire year postfire_precipitation_coefvar: Daily Coefficient of Variation (mm), fire date to end of fire year ls_factor: Slope length factor (unitless) KFACTWS_DC: Soil Erodibility Factor (k-factor, unitless) nbr_0_year: Normalized Burn Ratio for fire year (unitless) vpdX: April- October mean vapor pressure deficit for year X (kPa) defX: April- October mean Climatic water deficit for year X (mm) pptX: April- October total precipitation for year X (mm) tmaxX: April- October mean monthly temperature maximum for year X (°C) month: Month of fire x: Longitude y: Latitude name: Fire Name/code from MTBS.gov Missing data codes: NA