Examining the Impacts of Pre-Fire Forest Conditions on Burn Severity Using Multiple Remote Sensing Platforms

: Pre-fire environmental conditions play a critical role in wildfire severity. This study investigated the impact of pre-fire forest conditions on burn severity as a result of the 2020 Bighorn Fire in the Santa Catalina Mountains in Arizona. Using a stepwise regression model and remotely sensed data from Landsat 8 and LiDAR, we analyzed the effects of structural and functional vegetation traits and environmental factors on burn severity. This analysis revealed that the difference normalized burn ratio (dNBR) was a more reliable indicator of burn severity compared to the relative dNBR (RdNBR). Stepwise regression identified pre-fire normalized difference vegetation index (NDVI), canopy cover, and tree density as significant variables across all land cover types that explained burn severity, suggesting that denser areas with higher vegetation greenness experienced more severe burns. Interestingly, residuals between the actual and estimated dNBR were lower in herbaceous zones compared to denser forested areas at similar elevations, suggesting potentially more predictable burn severity in open areas. Spatial analysis using Geary’s C statistics further revealed a strong negative autocorrelation: areas with high burn severity tended to be clustered, with lower severity areas interspersed. Overall, this study demonstrates the potential of readily available remote sensing data to predict potential burn severity values before a fire event, providing valuable information for forest managers to develop strategies for mitigating future wildfire damage.


Introduction
Wildfires are an essential forest process, but contemporary fires pose significant challenges to forests, resulting in substantial alterations to their structure and function [1][2][3].Fires can contribute to soil erosion, reduced capacity for carbon storage, and changes in ecosystem processes and forest health [4,5].The Southwestern U.S. has witnessed a substantial increase in wildfire occurrence, generally attributed to extreme droughts and the accumulation of uncharacteristically abundant fuels [6].Additionally, the legacy of past fires (or their suppression) can impact the severity of subsequent fires by influencing vegetation types and the time available for fuels to regenerate [7].
Burn severity, defined as the extent of environmental change caused by fire [8], carries multiple ecological implications.Severity is the direct consequence of fire behavior, exhibits a direct correlation with carbon emissions [9], and has significant impacts on vegetation and soils as well as the potential to profoundly influence the trajectory of ecosystem responses [10].Given the extensive damage caused by uncharacteristically high-severity contemporary fires, accurately revealing the relationship between pre-fire conditions and burn severity could provide a valuable method to anticipate wildfire impacts [11] that is crucial for prioritizing fuel management efforts [12], refining the quantification of firecarbon dynamics [13], and forecasting ecosystem impacts [10].The magnitude of burn severity is co-influenced by environmental factors that regulate fire behavior including wind, relative humidity, fuel loads, forest density, live and dead fuel moisture content, and terrain characteristics, making it challenging to attribute specific variables to variations in burn severity [14].
Prior studies have demonstrated that burn severity hinges strongly on the dominant vegetation cover types and forest structure because these factors are central to fire behavior [15,16].Research focusing on the spatial variability of burn severity commonly represents vegetation and fuels primarily through digital maps of vegetation types, land cover, and terrain factors [17] including slope, aspect, and the eastness and northness of the slope.In areas with rugged terrain, these maps are sometimes supplemented with a limited set of stand structural variables measured in the field [12].However, this is frequently accomplished using plot data, which only sample a small fraction of the landscape.The increasing availability of light detection and ranging (LiDAR) data presents an opportunity to obtain additional insights into the horizontal and vertical structure of pre-fire vegetation, which can significantly enhance our understanding of how vegetation and fuels contribute to fire severity [18].
The influence of topography on fire severity is intricate and variable, often yielding conflicting findings.These intricacies arise from the complex interactions between topography, fuels, and fire behavior [19,20].Some research has indicated that extreme fire weather conditions can override or alter the relationships between topography and fire severity [21].Assessing the robustness of these findings can be challenging, as only a limited number of studies have incorporated fire progression data to evaluate the alignment between topographical features and the advancing fire front [19].
To establish the relationship between pre-fire forest attributes and burn severity, it is essential to explore their connections to the actual expressed severity.Finer spatial resolution forest structure datasets, obtained through CubeSats, UAVs, or LiDAR, are becoming increasingly available as well as the traditional coarser multispectral remotely sensed datasets [22].Coupled with machine learning classification methods and various resolution datasets, forest structural classification can help explain the variations in burn severity by focusing on pre-fire structural diversity [23].
The primary aim of this research was to investigate how pre-fire forest conditions impacted the burn severity in a study landscape and to identify which pre-fire variables accounted for variations in fire severity across the landscape.To achieve this objective, we identified a suite of pre-fire structural and functional traits using aerial LiDAR and Landsat 8 data.These traits were used as predictor variables to model the burn severity of the 2020 Bighorn Fire in the Santa Catalina Mountains (SCMs), one of the Sky Island mountain ranges in southern Arizona, USA.Understanding the relationship between these variables and severity, by modifying fire behavior, is crucial for managing and mitigating the impact of wildfires on our forests.

Study Area and Fire Regime
For centuries up until the late-19th century, the southwest region of the United States experienced a fire regime very different from the present day.At that time, fires were frequent, occurring every 5-20 years, and typically of low intensity.These fires were predominantly ignited by lightning strikes [24] or by Native American practices [25].Because fires were frequent, often burning for weeks under moderate conditions, fuel accumulation was limited, and consequently, the fire severity was generally moderate [26].However, toward the end of the 19th century, changes in land management such as livestock grazing and fire suppression policies brought about a significant transformation in the fire regime [27,28].Urban expansion and agricultural activities further altered landscapes, disrupting natural fire cycles and fuel distribution [29].During the mid-20th century, aggressive fire suppression efforts became commonplace, with the aim of extinguishing all wildfires, regardless of their natural role or size [30].The combined impact of fire suppression, overgrazing, and changes in land use resulted in a reduced area being burned, leading to the accumulation of fuel and the uniformity of landscapes.This led to a decrease in vegetation diversity and an increase in the prevalence of fire-prone species like Bromus tectorum (cheatgrass) and other non-native grasses, contributing to more widespread and devastating fires [31].
In 2020, the southwest region of the United States, comprising Arizona and New Mexico, experienced a substantial wildfire season that led to the burning of 432,355 ha and the average area burned annually in the two states over the preceding ten-year period [32].Wildfires in Arizona covered a significantly larger area (376,164 ha) in 2020 compared to its ten-year annual average of 123,681 ha.New Mexico recorded a lower area (56,191 ha) burned by wildfires in 2020 compared to its ten-year average of 121,313 ha.
At the onset of 2020, southern Arizona including the SCM experienced a severe drought, reflecting below-average precipitation levels in the region during the preceding winter and spring seasons as well as the effect of rising temperatures.In June 2020, the Palmer drought severity index (PDSI) registered at −2.2 near Mt.Lemmon, indicative of severe drought conditions.Conversely, June 2019 exhibited a PDSI value of −0.88, reflecting a state of neutrality regarding drought conditions [33].Preceding the occurrence of the fire, the mean temperature exceeded the typical average by more than 5.6 • C, accompanied by gusty winds surpassing 40 miles per hour colliding with the mountain ridges [34].Insufficient rainfall in the second half of 2019 initially resulted in a reduction in the growth of fine fuels and below-average levels of warm-season fine fuel carryover, along with low soil moisture.However, above-average winter precipitation in Tucson, Arizona during the 2019-2020 season (11 mm monthly average, exceeding the historical average by 9 mm) led to increased fine fuel growth [35].
The study site (Figure 1) is located in the SCMs, one of the Arizona Sky Islands (Lat: 32.4 • N, Long: −110.8 • E) [36].The SCMs exhibit a rugged topography characterized by steep slopes and deep canyons [37].The highest point in the Catalinas is Mt.Lemmon, which sits at an elevation of 2791 m above sea level; the elevation change from the base is 1942 m.The mountains have a semi-arid climate, with an average annual precipitation varying from 300 mm at the mouth of Sabino Canyon to 750 mm on Mt.Lemmon [38].Precipitation is generally bi-seasonal, with two modes in winter (December to March) and summer (June to September) [39].Winter precipitation typically comes from Pacific storms, while summer precipitation comes from monsoonal moisture from the Gulf of Mexico.The vegetation in this region showcases remarkable diversity, supporting mixed coniferous forests on north facing slopes at high elevations comprised of mixed stands of Douglas fir (Pseudotsuga menziesii), Southwestern white pine (Pinus strobiformis), and white fir (Abies concolor).On many south facing slopes at higher elevations, mixed-conifer forests including stands of Ponderosa pine (Pinus ponderosa) form a landscape mosaic.The deep canyons at higher elevations, particularly those with flowing water, provide ideal conditions for thriving hardwood forests consisting of Bigtooth maple (Acer grandidentatum), Aspen (Populus tremuloides), New Mexico locust (Robinia neomexicana), Arizona walnut (Juglans major), Gambel oak (Quercus gambelii), and Velvet ash (Fraxinus velutina) [40].A species-rich zone of Madrean oaks dominates middle elevations, and lower elevations are covered mostly by desert scrub vegetation dotted with saguaros (Carnegiea gigantea) [41].
The Catalina Mountains have experienced multiple landscape fires in the past 20 years (Figure 1).In 2002, the Bullock Fire burned approximately 11,722 ha, affecting the eastern half of the Catalinas.The following year, the Aspen Fire perimeter included around 28,081 ha, primarily targeting the western half.These two events burned extensive portions of the mountain range with varying severity [42].Much of the vegetation on the mountain persisted through these fires, where fire severity was low to moderate; in areas that experienced higher severity, vegetation has recovered in some areas by seedlings or the root sprouts of some species such as aspens, oaks, and shrubs, which are now regenerating in the previously mixed-conifer forests [43,44].In areas of high tree mortality, conifers (which rely on seed regeneration) may recover slowly, or the vegetation may convert to an alternative state [45].In 2009, the Guthrie Fire perimeter covered 1966 ha, and in 2017, the Burro Fire burned through 11,059 ha, re-burning a section of the Bullock Fire's footprint and most of the Guthrie Fire area.The most recent event is the 2020 Bighorn Fire, which burned an area of 48,157 ha, encompassing the majority of mid-and upper-elevation forest areas over the period 5 June-23 July 2020 (https://en.wikipedia.org/wiki/Bighorn_Fire)[32].The Catalina Mountains have experienced multiple landscape fires in the past 20 years (Figure 1).In 2002, the Bullock Fire burned approximately 11,722 ha, affecting the eastern half of the Catalinas.The following year, the Aspen Fire perimeter included around 28,081 ha, primarily targeting the western half.These two events burned extensive portions of the mountain range with varying severity [42].Much of the vegetation on the mountain persisted through these fires, where fire severity was low to moderate; in areas that experienced higher severity, vegetation has recovered in some areas by seedlings or the root sprouts of some species such as aspens, oaks, and shrubs, which are now regenerating in the previously mixed-conifer forests [43,44].In areas of high tree mortality, conifers (which rely on seed regeneration) may recover slowly, or the vegetation may convert to an alternative state [45].In 2009, the Guthrie Fire perimeter covered 1966 ha, and in 2017, the Burro Fire burned through 11,059 ha, re-burning a section of the Bullock Fire's footprint and most of the Guthrie Fire area.The most recent event is the 2020 Bighorn Fire, which burned an area of 48,157 ha, encompassing the majority of mid-and upper-elevation forest areas over the period 5 June-23 July, 2020 (https://en.wikipedia.org/wiki/Big-The analysis area for this research was determined by the availability of geographically concurrent aerial LiDAR before and after the 2020 Bighorn Fire.It primarily covered the southern portion of the mountains and included land burned in the Bighorn, Aspen, Bullock, Guthrie, and Burro Fires.The total area of study was approximately 22,000 ha, which was covered by the pre-fire LiDAR flight in 2019.
The Bighorn Fire reburned nearly two-thirds of the 2017 Burro Fire [35].While the Burro Fire was the only fire to have burned a portion of the area affected by the Bighorn Fire in the previous ten fire seasons, it is probable that areas that burned in other events more than 10 years ago were still in a state of post-fire recovery.Therefore, it is likely that the condition of the areas affected by the Bighorn Fire had been significantly impacted by previous fires, affecting the vegetation composition and fuel loads, and thus the current wildfire severity [46].

Vegetation Structural and Topographical Variables Derived from Aerial LiDAR
Figure 2 shows the workflow and data sources, processing, and analysis steps that were performed to conduct this research.The study utilized a pre-fire LiDAR dataset obtained from the USGS 3DEP, which meets the quality standard of level 2 (0.7 nominal pulse spacing, NPS) and has an average point density of 2 points per square meter.The lidR package in R was employed to extract individual trees from the point cloud [47].Using the package, we executed a sequence of processing steps including ground classification [48], digital terrain and surface model (DEM) generation [49], height normalization, canopy height model (CHM) generation [50], and individual tree segmentation [51].Elevation, slope, aspect, eastness and northness were derived from the LiDAR data, as shown in the below (Equations ( 1)-( 4)): Aspect = (57.29578 Remote Sens. 2024, 16, 1803 6 of 21 estimate can be used to indicate canopy diversity by calculating the standard deviation of each tree height in the grid [58].

Spectral Indices to Derive Burn Severity and Vegetation Functional Variables Using Satellite Imagery
Satellite-derived burn severity data are typically quantified by spectral changes between the pre-and post-fire satellite images.We acquired Landsat 8 satellite images of the study area for 14 July 2019 and 1 August 2020, which represent the pre-and post-fire conditions, respectively.There is an active debate regarding which of the two main equations, the delta normalized burn ratio (dNBR) and its relativized form (RdNBR), or other newer formulations (RBR) is most suitable for quantifying burn severity [59].We used fire severity indices generated from multi-date LANDSAT data, published by the Burned Area DEM data allow for the direct calculation of slope and aspect by analyzing the infinitesimal changes in elevation (dz) along the x and y cell directions ( dx, dy) within the dataset.
Eastness is determined by computing the sine of the aspect, revealing the east-to-west orientation of the slope, while northness is derived from the cosine of the aspect, indicating the north-to-south orientation of the slope [52].
CHM generation in LiDAR deploys triangulation-based algorithms that differ from point-to-raster algorithms.While point-to-raster algorithms are straightforward and enable user-defined resolution, they may produce empty pixels when the grid resolution is too fine for the available point density [47].In contrast, triangulation-based algorithms utilize a Delaunay triangulation to interpolate first returns, which is more intricate but eliminates empty pixels irrespective of the input data resolution.
Tree density, tree height, canopy cover, and canopy diversity (standard deviation of canopy height) are estimated from a 1-m CHM, which is adopted by the individual tree segmentation (ITS) method [53].To calculate tree density, tree top points are identified directly from the straightforward method from the CHM and canopy coverage [54].These methods operate on the premise that the apex of a tree is the tallest point within its crown, and the crown's periphery is comparatively lower [55].A local maximum (LM) filter, typically with a designated size sliding window, is employed in a CHM, derived from a point cloud to pinpoint LMs that indicate treetops [56].Despite the perceived challenge of CHM-based techniques in identifying overshadowed trees, they can serve as initial points for algorithms that outline tree crowns in a CHM or point cloud space [57].Information for each treetop was recognized as a point feature and transformed into density features once placed on a grid with a 30 m × 30 m fishnet pattern.This was conducted to facilitate a comparison with the burn severity derived from the Landsat 8 data.Each tree height estimate can be used to indicate canopy diversity by calculating the standard deviation of each tree height in the grid [58].

Spectral Indices to Derive Burn Severity and Vegetation Functional Variables Using Satellite Imagery
Satellite-derived burn severity data are typically quantified by spectral changes between the pre-and post-fire satellite images.We acquired Landsat 8 satellite images of the study area for 14 July 2019 and 1 August 2020, which represent the pre-and post-fire conditions, respectively.There is an active debate regarding which of the two main equations, the delta normalized burn ratio (dNBR) and its relativized form (RdNBR), or other newer formulations (RBR) is most suitable for quantifying burn severity [59].We used fire severity indices generated from multi-date LANDSAT data, published by the Burned Area Emergency Response (BAER) program at the Geospatial Technology and Applications Center in the USDA Forest Service.The program released normalized burn ratio (NBR) and dNBR burn severity datasets based on Landsat 8 imagery.The RdNBR index is available through RAVG (Rapid Assessment of Vegetation Condition after Wildfire), which was published in Aug 2021, one year after fire containment [60][61][62].
Burn ratio calculations are based on the NBR index.The NBR uses near-infrared (ρNIR) and shortwave infrared (ρSWIR) spectral reflectance data to highlight the difference between healthy vegetation (which reflects NIR radiation strongly and absorbs SWIR radiation) and burned areas (which reflect both types of types of radiation more evenly) [62].The NBR is calculated using the following formula: The NBR values range from −1 to 1.A value of −1 indicates the absence of vegetation such as in areas severely affected by fire.On the other hand, a value of 1 indicates areas with dense vegetation with no reduction in reflectance.The differenced NBR (dNBR), which is adjusted with the amount of pre-fire vegetation, is calculated using the following formula and has a broader range (−2 to +2): Sites that have experienced severe burning are represented by high positive values within this range (Table 1).It would be biased in areas with low pre-fire biomass, where even a high-severity fire might result in a low dNBR simply because there was not much vegetation to burn in the first place.In this study area, the area already burned by previous fire would be biased when only using dNBR.To address this issue, the relativized differenced NBR (RdNBR) normalizes the dNBR by the square root of the pre-fire NBR.This normalization reduces the bias in the dNBR, making the RdNBR a more accurate measure of burn severity in areas with low pre-fire biomass [59].The relative dNBR is computed by applying an offset and a scale factor to the dNBR [62]: The RdNBR is a measure of the change in burn severity relative to the amount of vegetation that was present before the fire.This generally provides a more representative measure of the actual impact of the fire on the landscape [61].
To understand the functional traits of the mountain ecosystem, we calculated the prefire normalized difference vegetation index (NDVI) and normalized difference moisture index (NDMI) using red, NIR, and SWIR reflectance values.These vegetation indices are commonly used to characterize plant health, biomass, and water content [63].These indices quantify the density and greenness of vegetation based on the difference in reflectance between the near-infrared (NIR) and red bands of the electromagnetic spectrum.It has been suggested that shortwave infrared is well-suited for the remote sensing of the canopy water content [64].Consequently, it might also reflect the pre-fire canopy conditions [65].
By analyzing the NDVI and NDMI values of an area before a fire, we can understand the baseline vegetation conditions when comparing them with structural traits from a LiDAR dataset and quantify how these vegetation and structural variables interact with each other to affect the burn severity.The indices were derived from Landsat 8 satellite imagery from the pre-fire scene in July 2019, prior to the 2020 Bighorn Fire, and were calculated using the following formula: The highest NDVI of 2019, a year prior to the fire, was determined by comparing the mean NDVI values over the year using Landsat 8 images of the boundary of the burned area.This analysis was performed through Google Earth Engine.The maximum index was included to represent the fuel loads and vegetation conditions before the fire.As a result, the Landsat 8 image from 2 October 2019, displays the peak mean NDVI value for the region.All spectral indices used in the analysis are shown in Figure 2.

Pre-Fire Land Cover Classification
Burn severity is significantly influenced by pre-fire environmental conditions, with a particular emphasis on land cover.As illustrated in Figure 1, a substantial portion of the study area had experienced multiple previous fires, leading to varying rates of recovery across different locations.Consequently, we used thematic land cover classification results to represent the pre-fire conditions of the forest.
To conduct thematic imagery classification and estimate basic land cover categories throughout the entire SCMs over a pre-fire time series, we employed a combination of passive and active remote sensing data sources (Figure 2).These sources included Planet Scope satellite images (4-band, 3-m resolution) and the USGS 3DEP LiDAR canopy height model.The Planet Scope satellite imagery was chosen for its frequent temporal coverage before as well as its relatively high spatial resolution of 3 m, surpassing that of Landsat 8.
We utilized pre-fire aerial LiDAR data to generate canopy height models for use as input in our classification process.Our approach involved employing these remote sensing data sources in a supervised thematic classification task, wherein individual pixels were categorized as belonging to either a vegetation or abiotic feature type.To accomplish this, we utilized the C5.0 decision tree classifier algorithm, implemented within RStudio, to assign pixels into one of four distinct classes (Table 2).

Herbaceous
Green vegetation spectral signature but have low height (<1 m).This includes green grasses, forbs, and ferns.
Decision trees represent a class of supervised learning algorithms employed extensively for the purposes of both classification and regression [66].The C5.0 algorithm employs an iterative, top-down, and recursive approach in the construction of decision trees [67].This method initiates the process by designating the entire dataset as the root node of the tree and subsequently divides the data into subsets based on the attributes or features, seeking to optimize the maximization of information gain or minimization of impurity as its guiding principle.In the context of our study, this functionality of C5.0 was harnessed to effectuate the classification of distinct land cover classes with the overarching goal of maximizing the information gain across these varying categories.
Given the acknowledged risk of overfitting stemming from an imbalanced distribution of samples [68], we adopted a training-validation sample ratio of 75% to 25% for each class.Additionally, we employed 58 reference points distributed across the study area to assess the overall accuracy of the final classification outcome.The land cover data were a crucial factor in quantifying the regression results, enabling us to discern variations in these results based on the attributes of the land cover classification outcomes.

Correlation Analysis and Stepwise Regression Using Mallow's C(p)
We employed Pearson correlation analysis to investigate the interrelations among various environmental factors (Figure 2).We used hierarchical clustering based on correlation coefficients among these variables based on the coefficient of correlation among the environmental variables.This analytical approach facilitated the discernment of linear relationships among paired variables.Moreover, it assisted in selecting the burn severity index (be it RdNBR or dNBR) that most accurately encapsulated the dominant trends and relationships within the dataset.
We used stepwise regression as a statistical method to select a subset of relevant predictor variables from a larger set of variables in a linear regression model.The objective was to identify a subset of predictors that would optimally explain the variation in the response variable while maintaining model simplicity [69].The leaps package in R provided tools to perform stepwise regression with different selection criteria such as the Akaike information criterion (AIC) or Bayesian information criterion (BIC) [70].We used this routine to perform stepwise regression, exhaustively exploring all possible combinations of predictor variables to identify the best subset based on the selected criterion (AIC or BIC).The best subset was the one that provided the best trade-off between model fit and complexity.We evaluated the best combination of regression model using Mallows' Cp, also known as "Mallow's C(p)" statistic [71,72], which measures the predictive accuracy used in the context of linear regression and model selection.It is used as part of the process for selecting the best subset of predictor variables in stepwise regression.The purpose of Mallows' Cp is to assess the goodness-of-fit of a model while considering its complexity.It is particularly useful when comparing different subsets of predictor variables in stepwise regression, where subsets with different numbers of predictors are evaluated.
The formula for Mallows' Cp statistic is as follows: where: SEE p is the sum of squared errors (residual sum of squares) for the model with p predictor variables; σ 2 is the estimated error variance for the full model (model with all predictor variables); n is the number of observations in the data; p is the number of predictor variables in the model.
Mallow's Cp value is contingent upon the variance between the anticipated error of a model with p predictors and the error variance of the complete model.A lower Mallows' Cp value signifies an optimal balance between the model's fit and complexity.When the Cp value is proximate to, or slightly exceeds, the number of predictors (p), it indicates a model that strikes a balance between complexity and simplicity.In our research, we examined many environmental variable combinations to identify the most suitable regression model that could account for variation in fire severity using this criterion.

Analyzing Residual Values with Geary's C
We used the regression model results to compute the residual values for dNBR, which represent the difference between the estimated burn severity value obtained from the stepwise regression model and the actual burn severity value derived from the Landsat satellite data (Figure 2).These residual values exhibit natural autocorrelation, as described by Legendre [73], due to the inherent spatial autocorrelation in ecological data.Spatial autocorrelation in burn severity arises because species at one location are strongly influenced by both biotic and abiotic factors, which tend to be spatially coherent over a landscape area.
Since wildfires cannot be replicated, it is crucial to maximize the number of data points (N) and implement blocking within burn severity categories to ensure that our statistical analysis yields valuable insights without being affected by pseudoreplication or spatial autocorrelation.This challenge is further complicated by landscape heterogeneity, which varies depending on the scale of analysis [74].Therefore, our approach in this study involved analyzing all residuals to identify potential autocorrelation patterns, rather than relying on small plots that may not adequately capture such autocorrelation.Global Moran's I and Geary's C are spatial autocorrelation coefficients used to assess the degree of spatial dependency among observations in geographical space [75].Global Moran's I quantifies the overall spatial autocorrelation within a dataset, indicating whether adjacent observations tend to exhibit similar (positive autocorrelation) or dissimilar (negative autocorrelation) values [76].While both metrics evaluate spatial autocorrelation, they exhibit differences in sensitivity and focus.Moran's I is generally more suitable for capturing global spatial trends, while Geary's C tends to be more responsive to local spatial patterns and variations because it utilizes the sum of squared distances, in contrast to Moran's I, which relies on standardized spatial covariance [77].Geary's C, due to its use of squared distances, is less influenced by linear associations and may detect autocorrelation that Moran's I might underweight [77].These two coefficients offer complementary insights into the underlying spatial structures and dependencies inherent in geographical data.In this study, we employed Geary's C to quantify the spatial pattern of residual values in relation to the final coefficient of determination obtained from the stepwise regression analysis.

Correlation between Environmental Variables and Burn Severity Indices
The relationship between RdNBR and dNBR from Landsat 8 images in 2020 was generally linear (Figure 3), with some outliers exceeding index values of 2000.Positive dNBR values signified a rise in near-infrared radiation reflectance post-fire compared to its pre-fire state.For both dNBR and RdNBR, most areas exhibited positive values.
Since wildfires cannot be replicated, it is crucial to maximize the number of data points (N) and implement blocking within burn severity categories to ensure that our statistical analysis yields valuable insights without being affected by pseudoreplication or spatial autocorrelation.This challenge is further complicated by landscape heterogeneity, which varies depending on the scale of analysis [74].Therefore, our approach in this study involved analyzing all residuals to identify potential autocorrelation patterns, rather than relying on small plots that may not adequately capture such autocorrelation.Global Moran's I and Geary's C are spatial autocorrelation coefficients used to assess the degree of spatial dependency among observations in geographical space [75].Global Moran's I quantifies the overall spatial autocorrelation within a dataset, indicating whether adjacent observations tend to exhibit similar (positive autocorrelation) or dissimilar (negative autocorrelation) values [76].While both metrics evaluate spatial autocorrelation, they exhibit differences in sensitivity and focus.Moran's I is generally more suitable for capturing global spatial trends, while Geary's C tends to be more responsive to local spatial patterns and variations because it utilizes the sum of squared distances, in contrast to Moran's I, which relies on standardized spatial covariance [77].Geary's C, due to its use of squared distances, is less influenced by linear associations and may detect autocorrelation that Moran's I might underweight [77].These two coefficients offer complementary insights into the underlying spatial structures and dependencies inherent in geographical data.In this study, we employed Geary's C to quantify the spatial pattern of residual values in relation to the final coefficient of determination obtained from the stepwise regression analysis.

Correlation between Environmental Variables and Burn Severity Indices
The relationship between RdNBR and dNBR from Landsat 8 images in 2020 was generally linear (Figure 3), with some outliers exceeding index values of 2000.Positive dNBR values signified a rise in near-infrared radiation reflectance post-fire compared to its prefire state.For both dNBR and RdNBR, most areas exhibited positive values.We found a robust correlation (r = 0.83) between the two burn severity indices, RdNBR and dNBR.This implies that employing either index for regression analysis might produce analogous outcomes.Vegetation indices manifest higher correlation coefficients, in contrast to the lower coefficients among terrain attributes such as eastness or northness (Figure 4).Notably, elevation emerged as a pivotal variable correlating with alterations in the vegetation indices (r = 0.85 with NDVI and r = 0.78 with NDMI).This may be attributed to the properties of the study locale, often referred to as a "Sky Island", which possesses clines in abundance of vegetation and fuel along the elevation gradient [36,78].In contrast, when comparing the burn severities with all other environmental variables, dNBR consistently exhibited higher correlation coefficients compared to the RdNBR across nearly all variables.Thus, dNBR appears to better represent the entire set of variables.
We found a robust correlation (r = 0.83) between the two burn severity ind RdNBR and dNBR.This implies that employing either index for regression analysis m produce analogous outcomes.Vegetation indices manifest higher correlation coefficie in contrast to the lower coefficients among terrain attributes such as eastness or north (Figure 4).Notably, elevation emerged as a pivotal variable correlating with alteration the vegetation indices (r = 0.85 with NDVI and r = 0.78 with NDMI).This may be attribu to the properties of the study locale, often referred to as a "Sky Island", which posse clines in abundance of vegetation and fuel along the elevation gradient [36,78].In cont when comparing the burn severities with all other environmental variables, dNBR sistently exhibited higher correlation coefficients compared to the RdNBR across ne all variables.Thus, dNBR appears to better represent the entire set of variables.

Stepwise Regression Analysis Using Land Cover Classification Results
The transitional zone (Figure 5), which lies between elevations of 1500 to 1800 m and is characterized by yellow and light blue colors, was previously affected by the Aspen Fire in 2003.The vegetation of this area is a mix of desert, woodland, and montane plants, many of which recovered during the post-fire era after the Aspen Fire.However, the northeast side of the mountain, which is highlighted by the red and orange colors in Figure 5, was previously damaged in 2002 by the Bullock Fire and showed more difference compared to the mid-elevation transition area.This area is mostly classified as a mixed-conifer and subalpine forest range, covered by tall trees.
in 2003.The vegetation of this area is a mix of desert, woodland, and montane plants, many of which recovered during the post-fire era after the Aspen Fire.However, the northeast side of the mountain, which is highlighted by the red and orange colors in Figure 5, was previously damaged in 2002 by the Bullock Fire and showed more difference compared to the mid-elevation transition area.This area is mostly classified as a mixedconifer and subalpine forest range, covered by tall trees.Considering the heterogeneous vegetation coverage within the study area, we evaluated areas with diminished residual values in a regression model.To gauge the sensitivity of dNBR to various environmental factors, we developed a stepwise regression model with dNBR values as the dependent variable and 13 other environmental factors as independent variables to elucidate disparities in the residual values amongst the groups.
The results show that different types of land cover contained different combinations of important variables in the stepwise regression model (Table 3).
Vegetation-associated indices and variables including the NDVI prior to the fire, canopy cover, tree density, and the standard deviation of tree height were significantly Considering the heterogeneous vegetation coverage within the study area, we evaluated areas with diminished residual values in a regression model.To gauge the sensitivity of dNBR to various environmental factors, we developed a stepwise regression model with dNBR values as the dependent variable and 13 other environmental factors as independent variables to elucidate disparities in the residual values amongst the groups.
The results show that different types of land cover contained different combinations of important variables in the stepwise regression model (Table 3).
Vegetation-associated indices and variables including the NDVI prior to the fire, canopy cover, tree density, and the standard deviation of tree height were significantly associated with burn severity across all land cover types.However, variables such as the NDMI and the maximum NDVI (NDVI_Max) did not exhibit a significant relationship with burn severity within areas covered by healthy vegetation.In unclassified land cover areas, aspect was incorporated within the amalgamation of variables within the stepwise regression model.Because most of the study area has a south-facing slope, the variable northness was not significantly correlated with burn severity of the Bighorn Fire across any classifications.Given the considerable discrepancy in sample sizes amongst the herbaceous (n = 4249), ground (n = 68,754), and alive tree (n = 65,611) areas, a direct comparison was not feasible.Correlation analysis indicates that NDVI exhibited the strongest correlation with changes in dNBR.Additional comparisons emerged when contrasting the outcomes of the regression model solely utilizing NDVI with those of the stepwise model.A simple linear regression analysis encompassing all land cover types using only NDVI as the predictor variable produced a lower coefficient of determination (r 2 ) compared to the coefficient derived from the stepwise model.Specifically, the coefficient obtained from the model employing solely NDVI was 0.33, compared to 0.16 for herbaceous cover, 0.26 for ground, and 0.25 for the live vegetation area.In general, these coefficients were marginally low; however, we noted more pronounced disparities in the case of the herbaceous area compared to the other two land cover categories.
The outcomes derived from the stepwise regression analysis indicated a notable association between all environmental variables and the dNBR burn severity index (Table 3).Each land cover type exhibited a comparable coefficient of determination, yielding a value of 0.30, while the significant variables within the stepwise regression model varied among the cover types.Notably, NDVI, canopy cover, and tree density emerged as significant variables across all land cover types, underscoring their substantial influence on burn severity outcomes.Consistent with findings from the correlation analysis, vegetationrelated variables such as NDVI prominently featured in the stepwise regression, whereas terrain variables such as slope and aspect demonstrated a negligible impact on burn severity, at least within the context of this study.
The frequency and spatial distribution of the residual values derived from the regression models were specific to three distinct land cover types (Figures 6 and 7).These residual values elucidate the disparities between the predicted and observed burn severity indices.In Figure 6, the distribution of residual values across all classes demonstrated adherence to a normal distribution, as evidenced by their alignment with the theoretical quantiles.These quantiles were calculated based on the assumed parameters of the distribution such as mean and standard deviation.However, deviations from this expected pattern were observed among the three different land cover classes.The residual values associated with the alive trees and shrubs class exhibited a departure from the expected straight line, indicating a discrepancy between the predicted and observed values.Conversely, the herbaceous class displayed a downward trend, suggesting lighter tails, while the bare ground/rock/sparse vegetation class exhibited an upward trend, indicative of heavier tails.
straight line, indicating a discrepancy between the predicted and observed values.Conversely, the herbaceous class displayed a downward trend, suggesting lighter tails, while the bare ground/rock/sparse vegetation class exhibited an upward trend, indicative of heavier tails.
Figure 7 provides a spatial representation of these residual value patterns.Elevated regions tended to display higher residual values, correlating with more severe burn severity indices.Lower elevation areas showed minimal differences between the predicted and observed residual values.These spatial trends aligned with the terrain characteristics of the sky island; areas of higher elevation typically featured larger biomass coverage, whereas lower elevation areas were characterized by shrubs and bare ground.Figure 7 provides a spatial representation of these residual value patterns.Elevated regions tended to display higher residual values, correlating with more severe burn severity indices.Lower elevation areas showed minimal differences between the predicted and observed residual values.These spatial trends aligned with the terrain characteristics of the sky island; areas of higher elevation typically featured larger biomass coverage, whereas lower elevation areas were characterized by shrubs and bare ground.

Calculating Spatial Distribution Using Geary's C
To ascertain the relationship between the residual values and land cover types, we evaluated the spatial autocorrelation within the dataset.The absence of spatial autocorrelation implies that the residual values were randomly distributed without interdependencies among the neighboring data points.Given that the land cover data indicated areas of homogeneous or similar cover types, we anticipated the presence of spatial autocorrelation within the dataset.
In the present study, Geary's C was computed to evaluate the spatial characteristics of the residual values, specifically investigating the presence or absence of spatial autocorrelation.Given that each set of residual values was derived from distinct regression models, each corresponding to a specific land cover type, we anticipated that spatial autocorrelation would manifest across all classes (Table 4).All Geary's C values were notably larger than 1, suggesting a strong negative spatial autocorrelation within each category.

Calculating Spatial Distribution Using Geary's C
To ascertain the relationship between the residual values and land cover types, we evaluated the spatial autocorrelation within the dataset.The absence of spatial autocorrelation implies that the residual values were randomly distributed without interdependencies among the neighboring data points.Given that the land cover data indicated areas of homogeneous or similar cover types, we anticipated the presence of spatial autocorrelation within the dataset.
In the present study, Geary's C was computed to evaluate the spatial characteristics of the residual values, specifically investigating the presence or absence of spatial autocorrelation.Given that each set of residual values was derived from distinct regression models, each corresponding to a specific land cover type, we anticipated that spatial autocorrelation would manifest across all classes (Table 4).All Geary's C values were notably larger than 1, suggesting a strong negative spatial autocorrelation within each category.

Discussion
When comparing dNBR and RdNBR through the distribution plots, we found that dNBR exhibited fewer outlier values in the examined study area.The results indicate that dNBR generally aligned more closely with the actual burn severity than the RdNBR for this area.However, as previously mentioned, the accuracy can fluctuate based on various environmental factors.While this study primarily examined the correlation between burn severity and pre-fire conditions, future research could benefit from comparing the burn severity to the actual changes in biomass.Such an investigation could shed light on the effectiveness of different burn severity indices in environments like sky islands, characterized by a diverse gradient of pre-fire vegetation types and structures with variable fire return histories.
Through our correlation and regression assessments, we found that certain vegetation and biomass attributes exhibited a stronger relationship with burn severity compared to other factors such as terrain.Prior research primarily concentrated on evaluating the fuel structure and loading solely through the utilization of dNBR indices [79], in contrast to the present investigation, which integrated an expanded array of structural and functional variables encompassing terrain and canopy attributes.Despite the relatively diminished overall predictive precision concerning the potential burn severity discerned within this study, it offers a means of identifying variables with greater likelihoods of predicting the burn severity.Additionally, the predictive accuracy exhibited variability, contingent upon the distinct ecosystem types that typify diverse landscapes [80,81].This study delineated the distinctive attributes of sky islands in the southwestern United States, characterized by a blend of vegetation amidst rugged mountainous terrain.
It is axiomatic that regions with more vegetation inherently have more fuel, leading to heightened burn indices [82].Previous studies have shown that densely forested areas with conifer trees have a higher regression value, while Mediterranean shrublands have a lower one.Generally, the correlation coefficient (r 2 ) for NDVI and dNBR values ranges from 0.5 to 0.8 [62,83,84].Given that our study zone encompassed a range from shrublands to thick forests situated on steep terrain, the derived R 2 value of 0.358 might appear counterintuitive at first glance.In a corresponding investigation, specifically in the context of a comparative analysis involving a field-based fire severity index (GeoCB)I and dNDVI derived from Landsat TM data, the R 2 value fell within the range of 0.3 to 0.5, aligning closely with the outcomes obtained within the present study [85].Nevertheless, the precise regression value can fluctuate based on various determinants such as the nature of pre-fire vegetation, burn severity, and the specific remote sensing equipment utilized.
The study indicates that the residuals of the stepwise regression model exhibit spatial variability across mountainous terrains.As the elevation increased, there was a marked increase in the residual values.This trend aligns closely with the inherent attributes of sky island mountains.Specifically, regions at lower elevations consist predominantly of sparse shrublands, while higher elevations are characterized by a dense canopy of coniferous trees.Notably, an intriguing pattern discerned from the vegetation classification map indicates a prevalence of herbaceous species at higher elevation areas.These herbaceous regions overlapped predominantly with zones that experienced previous fire events, notably the Bullock and Aspen Fires, indicating a possible post-fire community.Residual errors tended to be lower in proximity to these herbaceous zones in comparison to densely forested areas, even at comparable elevations.
Accurate discrimination between rocks and exposed terrain presented inherent challenges, necessitating their amalgamation into a unified classification category.Notably, the discernment of sparse vegetative features, particularly arid shrubs lacking pronounced green pigmentation, proved to be a formidable task, with heightened complexity observed at lower elevations within the cartographic representation.Consequently, the category encompassing sparse vegetation was amalgamated with the class designated for bare ground and rock.We acknowledge the potential for the occasional misclassification of isolated boulders as deceased arboreal or shrub entities.

Conclusions
This study presented a novel approach for investigating the predictability of post-fire burn severity using pre-fire structural and functional vegetation traits.We leveraged a well-established regression model retrospectively to identify a viable methodology for such predictions.Among the various factors analyzed, vegetation attributes emerged as the most significant predictors of burn severity.Stepwise regression analysis revealed that vegetation-related indices, NDVI, canopy cover, and tree density were significant variables across all land cover types.These variables were highly correlated with burn severity values throughout the mountain range.In the scope of fire behavior, NDVI acts as a proxy for essential functional traits influenced by drought conditions, and tree density and canopyrelated variables reflect the available fuel loads in the mountain ranges.It also implies that the areas with denser vegetation and higher levels of greenness in this mountain range had experienced more severe burns.It is worth noting that the study sites have a history of multiple fires, resulting in more diverse environmental conditions compared to other potential areas.Despite this, we successfully extracted variables indicative of functional traits from the LiDAR and satellite datasets, demonstrating valuable methods for future research in this direction.
Not all study locales consistently furnish adequate geospatial data, thus imposing limitations on the evaluation of post-fire severity.Nevertheless, the judicious adoption of the stepwise regression methodology, as demonstrated in this investigation, offers a viable means to retroactively estimate the impact of pre-fire conditions on post-fire burn severity by utilizing widely accessible satellite datasets such as Landsat or CubeSat.In this framework, these findings underscore that a comprehensive array of pre-fire vegetation structural, functional, and environmental variables can effectively facilitate the estimation of burn severity.However, caution must be exercised in extrapolating these methodologies to other regions and ecosystem types, as variations in tree species and terrain characteristics necessitate nuanced adaptations.
While this study examined numerous variables, the advent of newly developed sensors and satellites promises a wealth of additional information.It is imperative, however, to scrutinize these advances to learn more about what they may tell us regarding the factors that govern burn severity and ascertain the extent of their influence, particularly amidst the backdrop of recent extreme climate patterns in the southwestern United States.These insights not only augment our comprehension of burn severity dynamics, but also catalyze further inquiry in this field, thereby fostering advancements in both knowledge and methodology.

Figure 1 .
Figure 1.Boundaries of significant fires since 2002 in the Santa Catalina Mountains north of Tucson, AZ, USA.The area of the 2020 Bighorn fire (red polygon) encompassed most areas that had already experienced burning within the last two decades.The study area (black outline) is located on the south side of the mountain.

Figure 1 .
Figure 1.Boundaries of significant fires since 2002 in the Santa Catalina Mountains north of Tucson, AZ, USA.The area of the 2020 Bighorn fire (red polygon) encompassed most areas that had already experienced burning within the last two decades.The study area (black outline) is located on the south side of the mountain.

Figure 2 .
Figure 2. Structural and functional traits of the forest and topographic variables were extracted from multiple data sources (marked with blue circles) including aerial LiDAR point cloud and multispectral Landsat 8 satellite spectral surface reflectance data.

Figure 2 .
Figure 2. Structural and functional traits of the forest and topographic variables were extracted from multiple data sources (marked with blue circles) including aerial LiDAR point cloud and multispectral Landsat 8 satellite spectral surface reflectance data.

Figure 3 .
Figure 3. dNBR and RdNBR for the Bighorn Fire for the Santa Catalina Mountains (2020).The values of dNBR and RdNBR were multiplied by 1000 for visualization purposes.A total of 138,312 points were extracted from a 30-m grid derived from Landsat 8 images.Negative values indicate regrowth following the fire, while positive values indicate a higher degree of burn severity.

Figure 4 .
Figure 4. Correlation of all structural and functional variables.Correlations among the t graphic and pre-fire vegetation explanatory variables and the two burn severity indices (Rd and dNBR), where the dNBR represents a higher correlation coefficient overall.

Figure 4 .
Figure 4. Correlation of all structural and functional variables.Correlations among the topographic and pre-fire vegetation explanatory variables and the two burn severity indices (RdNBR and dNBR), where the dNBR represents a higher correlation coefficient overall.

Figure 5 .
Figure 5. Vegetation classification map, classified by Planet Scope and USGS LiDAR using classification and regression tree (CART), prior to the Bighorn Fire (overall accuracy 0.88).This indicates a variety of vegetation types across the study area.Lower elevations are characterized by sparse vegetation, whereas as the elevation increases, a richer diversity of vegetation species can be observed.

Figure 5 .
Figure 5. Vegetation classification map, classified by Planet Scope and USGS LiDAR using classification and regression tree (CART), prior to the Bighorn Fire (overall accuracy 0.88).This indicates a variety of vegetation types across the study area.Lower elevations are characterized by sparse vegetation, whereas as the elevation increases, a richer diversity of vegetation species can be observed.

Figure 6 .
Figure 6.Distributions of residuals of all classes.It is pertinent to acknowledge that the three land cover types under consideration exhibited unequal population sizes.The live trees/shrubs class demonstrated the smallest residual errors in contrast to the remaining two vegetation types.

Figure 6 .
Figure 6.Distributions of residuals of all classes.It is pertinent to acknowledge that the three land cover types under consideration exhibited unequal population sizes.The live trees/shrubs class demonstrated the smallest residual errors in contrast to the remaining two vegetation types.

Figure 7 .
Figure 7. Residual values representing the difference between the predicted and observed burn severity index values.Larger values were observed in higher elevation areas, where more severe burn severity was observed.Smaller differences were observed between the predicted and observed dNBR values in the lower elevation area.

Figure 7 .
Figure 7. Residual values representing the difference between the predicted and observed burn severity index values.Larger values were observed in higher elevation areas, where more severe burn severity was observed.Smaller differences were observed between the predicted and observed dNBR values in the lower elevation area.

Table 2 .
Thematic classes for vegetation cover classification.

Table 3 .
Stepwise regression model result based on land cover classification.The significance of variables depends on differences in cover types; not all variables are always significant.(* indicates significance, p < 0.05).

Table 4 .
Comparison of spatial autocorrelation between land cover types.