Introduction 

Wildfire regimes across the globe are changing significantly, with fires becoming more frequent, severe, and widespread than in previous decades; current models predict continued increases due to climate change and forest management practices ([20]). This intensification represents a major shift in how fire behaves in natural systems, moving beyond historical patterns that ecosystems and species have adapted to over thousands of years ([46]). The combination of climate change effects, such as rising temperatures and altered precipitation patterns, with increased human pressures has created conditions that promote more destructive fire behavior, a trend particularly evident in Mediterranean-type ecosystems ([35], [8]). This global surge has been particularly dramatic since 2017, with unprecedented fire behavior occurring in the Mediterranean Basin ([49]), Siberia ([27]), Australia, the Amazon basin, and the Western US ([14]).

Anthropogenic activities, specifically land-use shifts, management practices, and human-caused climate change, serve as critical drivers of global wildfire intensification ([1]). These land-use shifts are often driven by socioeconomic factors such as rural abandonment and the loss of agricultural mosaics, which can lead to large accumulations of combustible biomass or the expansion of flammable plantations ([34], [21]). This results in vast accumulations of combustible material, making the landscape more prone to burning and leading to increasingly large, catastrophic wildfires ([21], [4]). Concurrently, the expansion of human development into the wildland-urban interface (WUI) has significantly increased both ignition sources and fire risks, as proximity to settlements facilitates higher anthropogenic ignition density ([31], [42]). In many regions, this is coupled with reduced cultural or prescribed burning, which, along with high-biomass invasive species, has changed spatial patterns of wildfire risk ([44], [4]). These changes depart from indigenous fire stewardship practices that historically maintained ecosystem balance ([21]).

Burn severity represents a fundamental measure of fire’s environmental impact, defined as the extent of fire-induced change to physical ecosystem components such as vegetation and soil ([22], [40]). This metric is most commonly operationalized through the Differenced Normalized Burn Ratio (dNBR), which allows for a standardized assessment of fire-induced ecological change across diverse forest types ([29]). This metric goes beyond simple fire occurrence to capture the degree of environmental change, making it a critical factor for understanding post-fire ecological dynamics ([40]). The need to better understand factors controlling fire severity stems from concerns about public safety, infrastructure, critical wildlife habitat, and watershed health ([39]). Research has consistently shown that fire severity is controlled by both top-down weather and climate factors and bottom-up spatial controls from fuels and topography ([41]). However, these relationships are complex, as burn severity is co-influenced by forest structural attributes that dictate flammability traits, such as fuel loads and density ([54]).

Modeling the environmental drivers of wildfire is a critical scientific field, as confirmed by global datasets such as the Global Fire Emissions Database (GFED) and MODIS burned area products, which have been developed specifically to link fire occurrence to key environmental drivers, including climate, vegetation, topography, and human factors ([51], [13]). To capture the non-linear and complex interactions among these diverse environmental drivers, recent studies have increasingly used machine learning models to predict fire outcomes with high accuracy. These models frequently utilize vegetation indices, such as the Normalized Difference Vegetation Index (NDVI), as robust proxies for fuel biomass and ecosystem state ([10]). While prior machine learning studies have successfully modeled burn severity using short-term pre-fire conditions ([12], [47]), they often overlook long-term vegetation decline. This is a significant limitation, as decadal NDVI trends reveal gradual forest degradation and mortality that fundamentally alter fuel characteristics and forest structure over time, making fire behavior in chronically stressed forests harder to predict ([23]). To address this, our study focuses on the Zagros Mountains, where forest decline is operationalized as decadal NDVI trends (2014-2023 Landsat data), serving as a proxy for chronic vegetation stress associated with Zagros oak decline syndrome ([28]). By integrating decadal forest degradation into wildfire modeling, this study advances the current state of the art from short-term predictive snapshots to a more holistic ecological framework. Therefore, the primary objectives of this study are: (i) to statistically test the hypothesis that long-term forest decline impacts recent wildfire burn severity; (ii) to assess the relative importance of this long-term trend against critical short-term variables; and (iii) to develop a comprehensive regression model that analyzes the complex interactions between these ecological factors and human factors (distance to roads and settlements) within the Zagros ecosystem.

This study offers practical value for Zagros forest management by developing a predictive framework to map vulnerability hotspots, integrating chronic decline with biophysical risks. It supports local stakeholders (forest services, communities) in preventive strategies like fuel management prioritization and early interventions, enhancing resilience against escalating fire threats ([28]).

  Material and methods 

Study area

The Kashkan Watershed is situated within the Central Zagros Mountains, which host Iran’s most extensive forest ecosystem, oak-dominated woodlands covering 5-6 million hectares across western Iran and representing 40-44% of the nation’s total forest area ([9]). These ancient semi-arid systems (origins ~5500 years ago - [11]) play a vital role in environmental balance ([28]), but face escalating wildfire frequency and human disturbances, including land use changes, unregulated grazing, deforestation, and tourism ([33], [43])

The study was conducted in the Kashkan Watershed (approximately 9540 km²) of Lorestan Province, Central Zagros Mountains, western Iran (Fig. 1). This region features rugged topography (elevation 500-2500 m a.s.l.) and a semi-arid Mediterranean climate, with cool-wet winters (mean annual precipitation ~500 mm) and hot-dry summers (mean temperature ~30 °C) that promote intense fire seasons. Vegetation comprises a mosaic of open oak woodlands, shrublands, grasslands, and rainfed croplands, dominated by Persian oak (Quercus brantii - [28], [45]).

Fig. 1 - Geographical location of the study area. The map illustrates the Kashkan Watershed within the Central Zagros Mountains, Lorestan Province, western Iran. The main panel shows the region’s rugged topography, while the inset map indicates its location in Iran. This area has a semi-arid Mediterranean climate and a vegetation mosaic dominated by Quercus brantii woodlands.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The area is an ideal site for fire-ecosystem analysis because of overlapping pressures: high wildfire susceptibility ([45]) and ongoing Zagros oak decline, which has simplified fuel structures from dense to open woodlands and shrublands ([28]). A major wildfire in June 2024 (~5000 ha burned - [18]) provides a documented case for evaluating burn severity.

Dependent variable: burn severity (dNBR)

Burn severity was mapped using dNBR from Sentinel-2 MSI Level-2A surface reflectance ([58]). Cloud-free composites were created for pre-fire (10-20 June 2024) and post-fire (23 June-12 July 2024) periods to capture the event while minimizing phenological bias. Clouds/ cirrus were masked via the QA60 bitmask. NBR was calculated as per eqn. 1 ([50]):

and dNBR as per eqn. 2:

where NBR_pre and NBR_post are the NBR values pre-fire and post-fire, respectively (higher values = greater severity - [50]). Rasters were reprojected to a 30-m UTM grid (EPSG: 32638, Zone 38N) for alignment. The nearest-neighbor resampling method was applied to all layers to preserve the original radiometric integrity of the spectral bands and avoid the artificial smoothing of fire-periphery pixels that results from bilinear interpolation. Land-cover classes such as water bodies and barren lands, which comprised a negligible fraction of the study area (~1.7%), were retained in the modeling phase to ensure a spatially continuous dataset. While these surfaces occasionally exhibited spectral artifacts (false-high dNBR values) due to moisture or turbidity changes, their limited frequency ensured they did not significantly influence the overall performance or generalization of the machine learning model. Additionally, we retained negative dNBR values to preserve the full statistical distribution of the dataset, allowing the model to learn from the complete range of post-fire spectral responses.

Independent variables (predictors)

To model burn severity drivers, we selected predictor variables encompassing long-term vegetation dynamics, pre-fire biophysical conditions, topography, fuel structures, and anthropogenic proximity; all features are detailed in Tab. 1. To address our primary research aim regarding forest decline, we operationalized long-term vegetation dynamics as annual May-September NDVI trends (2014-2023). We derived these trends on a pixel-wise basis from Landsat-8/9 Collection 2 imagery using ordinary least squares (OLS) slope calculation with the “linearFit” reducer in Google Earth Engine^® (GEE). To ensure data quality, we applied cloud and shadow masking via the QA PIXEL band. We further mitigated potential atmospheric artifacts by calculating the annual mean NDVI for each growing season before trend estimation. Negative slope values indicate chronic vegetation decline ([57]). To capture immediate pre-fire conditions, we utilized the Normalized Difference Moisture Index (NDMI) derived from multi-scene median composites of Sentinel-2 imagery (June 10-20, 2024) to represent canopy moisture (eqn. 3):

Additionally, Land Surface Temperature (LST) was calculated from multi-scene mean composites of Landsat-8 Collection 2 imagery (May 1 - June 20, 2024). These biophysical metrics were aggregated using temporal averaging (median for NDMI and mean for LST) to ensure a representative baseline that aligns with the phenology of the dNBR evaluation period ([50]).

We incorporated topographic controls on fire spread using a 30-m SRTM Digital Elevation Model to extract slope and aspect, with the latter encoded into sine and cosine components to account for its circular nature ([56]). To represent fuel models, we utilized the ESA WorldCover 2021 (10 m) product, reclassified into five fuel-oriented categories: Closed Forest, Open Forest, Shrub/Grassland, Cropland, and Bare/Water. The decision to merge shrublands and grasslands into a single class was based on their ecologically similar fire behavior in semi-arid ecosystems. Both vegetation types serve as fine surface fuels with high combustibility and rapid drying rates, acting as the primary carriers of fire. In the Zagros context, these classes often form a continuous surface fuel complex that facilitates fire spread, justifying their consolidation for improved modeling robustness. The categorical data were one-hot encoded (drop-first) to prevent multicollinearity ([19], [59]). Spatial structural heterogeneity was quantified using the Gray-Level Co-occurrence Matrix (GLCM) contrast. To ensure spatial consistency with the 30-m dependent variable (dNBR) and other environmental predictors, the Sentinel-2 NIR band was first aggregated to 30-m resolution before texture computation. We then applied a 3×3 kernel to capture stand-level horizontal fuel continuity and canopy complexity, providing the machine learning model with contextually relevant structural information ([61]). Finally, we mapped human access gradients, key drivers of ignition probability, using vector layers of road networks and settlements. Euclidean distances to the nearest road and village were calculated at a 30-m resolution across the entire study area to ensure a continuous and accurate distance gradient. We then masked these proximity layers to vegetated and burnable surfaces during final data stacking. To accommodate the non-linear decay effects of these proximity variables within our models, we retained both raw distances and their log1p-transformed versions ([12], [24]). We aligned and stacked all spatial layers into a unified raster mosaic for subsequent sampling.

Sampling strategy and statistical analysis

To construct a robust analytical dataset from the continuous raster mosaic, we employed a stratified random sampling approach based on burn severity classes. We stratified dNBR values into four categories to ensure a balanced, representative sample across the full spectrum of fire effects: Unburned/Regrowth (dNBR < 0.10), Low (0.10 ≤ dNBR < 0.27), Moderate (0.27 ≤ dNBR < 0.44), and High severity (dNBR ≥ 0.44). This stratification ensured balanced representation, extracting five points per class per square kilometer. Following multiple iterations and strict quality control to remove non-finite values, the final dataset comprised 89,314 valid points ([58] - see Appendix 1 in Supplementary material).

For statistical inference, we established a baseline using Ordinary Least Squares (OLS) regression. Multicollinearity was assessed by computing Variance Inflation Factors (VIF) after all transformations (log1p) and dummy encoding (drop-first) were applied. We used an iterative screening process in which we recursively removed the variable with the highest VIF and recalculated VIF values for the remaining predictors until all variables met the threshold (VIF < 8). Finally, we applied HC3 robust standard errors to account for potential heteroscedasticity ([25], [26]). To capture complex, non-linear relationships, we employed Random Forest and XGBoost algorithms. To ensure optimal performance and prevent overfitting, we tuned hyperparameters using a grid search across a predefined parameter space (e.g., tree depth, number of estimators, and learning rate). We then selected the best-performing configurations for the final models: Random Forest (500 estimators, max depth of 15) and XGBoost (600 estimators, learning rate of 0.08). Model performance was rigorously evaluated using five-fold cross-validation (4×4 grid bins) to account for spatial autocorrelation, relying on cross-validated R² and Root Mean Square Error (RMSE) to assess generalization ([3], [5], [47], [38]). Furthermore, feature importance was extracted to rank driver influence, utilizing “Mean Decrease in Impurity” (MDI) for the Random Forest model and the “Gain metric” for the XGBoost model. Finally, non-parametric tests (Mann-Whitney U and Kruskal-Wallis) were utilized to assess group-level differences, as preliminary Shapiro-Wilk tests indicated that the dNBR data violated the assumption of normality.

All geospatial data were processed in Google Earth Engine (GEE) for consistency and reproducibility ([32] - see also Appendix 1 in Supplementary material).

  Results 

Descriptive statistics and spatial patterns

The cleaned analytical dataset contains 89,314 valid records. Descriptive statistics for all continuous variables are summarized in Tab. 2 as the basis for subsequent inference and modeling. The dependent variable dNBR showed a mean of 0.25 with substantial spread (standard deviation, SD = 0.20), indicating wide heterogeneity in burn severity across the landscape. The long-term decadal NDVI trend prior to the fire was near-neutral on average (mean = 0.01; SD = 0.01), suggesting overall stability with localized decline or recovery signals rather than a strong basin-wide trend.

Pre-fire biophysical conditions depict moderately warm and variably moist canopies, with Land Surface Temperature (LST) averaging 33.19 °C (SD = 6.59) and NDMI averaging 0.08 (SD = 0.22), consistent with semi-arid summer conditions that modulate available fuel moisture. The landscape exhibits complex terrain, with a mean slope of 9.19° (SD = 8.72) and an aspect covering the full directional range (mean = 181.88°), creating diverse physical and microclimatic conditions that locally influence fire spread and resulting burn-intensity patterns.

Human-access gradients exhibit pronounced right-skewness typical of proximity variables: mean distance to the nearest road is 1976.79 m (SD = 2346.26), and mean distance to the nearest settlement is 1717.55 m (SD = 2163.50), motivating the use of log1p transformation in the learning models, along with raw distances, to capture non-linear decay effects. The transformed counterparts reflect this stabilization, with the log-transformed distance to roads averaging 6.84 (SD = 1.38) and the log-transformed distance to villages averaging 6.99 (SD = 0.93). This transformation improves numerical behavior without losing critical information at zero distances.

Fine-scale structural heterogeneity of the pre-fire landscape is evident in the spatial texture distributions. Specifically, texture contrast exhibits a high mean and substantial dispersion (mean = 169.02; SD = 346.05), texture homogeneity is concentrated near low values (mean = 0.18; SD = 0.11), and texture variance is broadly spread (mean = 83.45; SD = 168.94). These patterns collectively support the inclusion of a single screened texture metric (contrast) to reduce redundancy while retaining critical signals of fuel continuity.

The land cover composition over the 9540.20 km² study area is illustrated in Fig. 2. The landscape was predominantly characterized by Open Forest, which constituted 63.10% (6017.50 km²) of the total area. Shrub/Grassland was the second most prevalent category, covering 24.30% (2320.90 km²). Other significant land cover types included Cropland at 9.70% (927.3 km²) and Closed Forest at 2.70% (254 km²). The Bare/Impervious/Water class represented a negligible fraction of the total area at just 0.20% (20.5 km²).

Fig. 2 - Distribution of land cover and fuel types in the Kashkan Watershed. The donut chart illustrates the area (km²) and percentage of the five reclassified ESA WorldCover 2021 categories used in the modeling process. Open Forest and Shrub/Grassland represent the dominant fuel complexes in the study area, while non-burnable surfaces (Bare/Water) constitute a negligible fraction (0.2%).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Spatial distribution of key variables

The spatial distribution of burn severity (dNBR) across the basin was notably heterogeneous (Fig. 3a). While the landscape was predominantly affected by low and moderate-low severity burns, several distinct clusters of high severity were observed, particularly in the western and northern portions. Ecologically, this pattern corresponds spatially with the decadal NDVI trend (Fig. 3b), where many of the high-severity clusters coincided with areas exhibiting significant long-term vegetation decline. This visual correlation supports our hypothesis that prolonged ecosystem stress, which likely increases dead biomass and alters fuel continuity, exacerbates fire impacts. Conversely, areas characterized by vegetation stability or improvement were generally associated with lower burn severities. Pre-fire biophysical conditions also showed clear spatial structuring: Land Surface Temperature (LST) was highest in the lower-elevation, sparsely vegetated areas in the south, while NDMI indicated higher canopy moisture content in the densely forested northern highlands (Fig. 3c, Fig. 3d). Together, these spatial patterns highlight how moisture deficits and chronic forest decline converge to create pockets of high vulnerability, providing crucial context for interpreting the regression models.

Fig. 3 - Spatial distribution of environmental and fire-related variables. (a) dNBR: classified into Unburned (less than 0.10), Low (0.10-0.27), Moderate (0.27-0.44), and High severity (greater than or equal to 0.44) based on Key & Benson ([22]), adapted for Zagros canopy density. (b) NDVI Trend: 10-year Sen’s slope (2014-2023) categorized into Significant Decline (less than -0.04), Slight Decline (-0.04 to -0.01), Stable (-0.01 to +0.01), and Improvement (greater than +0.01). (c) LST: classified by 10 °C intervals to represent regional thermal gradients. (d) NDMI: thresholds represent moisture stress (less than 0) to high canopy moisture (greater than 0.4). (e) Land Cover: derived from ESA WorldCover 2021 standard classes.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Model performance and feature importance

Statistical comparisons of pixel groups further reinforced these spatial observations. Burn severity was significantly higher in areas with a negative long-term NDVI trend compared to areas with stable or positive trends (Z = 18.81, p < 0.001, r = 0.06). Group summaries showed higher severity for declining pixels (mean/median dNBR = 0.26/0.28, N = 38.619) than for healthy or stable pixels (mean/median = 0.24/0.23; N = 50.695). This quantitatively confirms the linkage between previous vegetation decline and elevated burn severity (Tab. 3). Across the primary land-cover classes, burn severity differed significantly (Kruskal-Wallis H = 13.616.97, p < 0.001). To account for the large sample size and ensure that our findings reflect meaningful ecological variation, we calculated the effect size using epsilon-squared (ε²), which was found to be 0.15, indicating a moderate and robust relationship between fuel type and fire severity. The shrub and grassland class exhibited the highest central tendency (mean = 0.30), while closed forest showed the lowest (mean = 0.07 - Tab. 4). This pronounced contrast, even after excluding non-burnable surfaces to mitigate spectral noise, highlights the role of fuel structure; the fine, contiguous surface fuels in shrublands facilitate more intense burning compared to the more moisture-buffered closed-canopy oak stands.

Fig. 4 illustrates the out-of-sample performance under 5-fold cross-validation, demonstrating that XGBoost achieved the best generalization (cross-validated R² ≈ 0.70, RMSE ≈ 0.11), followed by Random Forest (cross-validated R² ≈ 0.64, RMSE ≈ 0.12). Both models displayed high stability across data splits (R² standard deviation ≈ 0.00). The strong predictive capacity of these tree-based models indicates that burn severity in the Zagros ecosystem is driven by complex, non-linear interactions among predictors. The random k-fold cross-validation used here may yield optimistic performance metrics because of the strong spatial autocorrelation inherent in raster data. Consequently, the cross-validation metrics reported in Tab. 4 should be interpreted as upper-bound estimates of generalization.

Fig. 4 - Comparative performance analysis of Random Forest (RF) and XGBoost models. (a) Model goodness-of-fit based on training and 5-fold cross-validated R² values. (b) Predictive error based on mean cross-validated RMSE. The red values above the bars in panel (b) represent the standard deviation of R² across cross-validation folds, indicating the high stability and consistency of both ensemble learners.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Regarding variable importance (Tab. 5), the Random Forest model identified moisture (NDMI) as the most influential predictor, followed by thermal stress (LST) and slope. Secondary drivers included proximity to roads, the decadal NDVI trend, and spatial texture contrast. For the XGBoost model, categorical land cover emerged as the primary feature, followed by NDMI, slope, and distance to roads.

Multicollinearity among continuous predictors was generally low to moderate (Tab. S1 in Supplementary material). For instance, LST and NDMI showed Variance Inflation Factors (VIF) of 2.74 and 3.17, respectively, while topographic, spatial texture, and proximity variables remained near 1.00 to 1.60. As expected, land-cover indicators exhibited larger VIFs (up to 7.19 for the combined shrub and grassland class). These diagnostics justify our use of tree-based learners, which robustly handle remaining categorical collinearity while effectively capturing the underlying ecological drivers of fire severity.

Linear regression results

The ordinary least squares (OLS) baseline model offered an interpretable approximation of burn-severity drivers, yielding a moderate fit (R² = 0.34, adj-R² = 0.34, F = 5041, p < 0.001, N = 89.314). While these linear relationships explain roughly one-third of the variance in dNBR, the substantial performance gap between OLS and the ensemble algorithms highlights that the ecological mechanisms driving fire severity in the Zagros Mountains are fundamentally non-linear and rely heavily on the synergistic effects of moisture gradients, terrain, and long-term degradation.

  Discussion 

Interpretation of key findings

The results indicate that machine learning algorithms, particularly tree-based ensembles, provide a robust predictive framework for modeling burn severity (dNBR) within the heterogeneous Zagros landscape ([47]). However, compared to simpler parametric models, these enhanced predictive capacities should be interpreted contextually, bearing in mind that random cross-validation yields upper-bound estimates due to spatial autocorrelation.

These values are particularly noteworthy because the analyzed fire event was dominated by low-to-moderate severity, which compresses the dynamic range of the response variable and lowers the theoretical ceiling for explained variance. The ability of these models to recover spatially coherent patterns under such circumstances indicates non-linear structures and complex interactions among the drivers of burn severity, which linear models cannot capture.

Model explainability from both ensembles converges on a moisture-thermal mechanism as the primary driver, which is further modulated by fuel type and its spatial structure ([47]). In the Random Forest model, moisture (NDMI) and temperature (LST) emerged as the most influential predictors, followed by slope, proximity to roads, the decadal NDVI trend, and spatial texture contrast. Consistent with well-established literature, this hierarchy confirms that pre-fire fuel moisture and surface temperature impose first-order constraints on fire effects, while topography and anthropogenic access gradients shape where intense burning is most likely to occur ([12], [41]). The XGBoost model corroborates the primacy of NDMI but elevates the importance of categorical land cover classes, most notably bare/impervious land (class 5.0) and, to a lesser extent, cropland (class 4.0) and the merged shrubland/open-forest class (classes 2 and 3), above many continuous predictors. Collectively, these results suggest a multi-layered dynamic. At the broader watershed scale, fuel moisture (NDMI) and thermal stress (LST), the top continuous predictors in our models, establish the baseline environmental susceptibility. Conversely, at the local patch level, specific fuel types and spatial configurations drive the amplification or damping of burn severity. This local modulation is empirically supported by the XGBoost model, which prioritized categorical land cover as the most critical feature, and the Kruskal-Wallis test results (H = 13.616.97, p < 0.001). Specifically, the significant contrast between the high severity in contiguous fine fuels (Shrublands/Grasslands, mean = 0.30) and the lower severity in more resilient, moisture-buffered structures (Closed Forest, mean = 0.07) provides clear evidence of how local fuel characteristics modify the fire’s impact.

These model-derived patterns are reinforced by independent statistical tests. The Mann-Whitney U test revealed that pixels with a negative long-term NDVI trend experienced significantly higher dNBR values ([48]), directly supporting the hypothesis that chronic vegetation decline predisposes ecosystems to more severe fires by reducing live fuel moisture and increasing the load of receptive fine fuels. Furthermore, the Kruskal-Wallis test confirmed significant differences in burn severity among land cover classes. Among vegetated types, shrubland/grassland exhibited the highest median severity, followed by open forest and cropland, while closed-canopy forest was the least severely burned ([30]). This pattern aligns with expected differences in fuel continuity, loading, and moisture-buffering capacity across these vegetation types.

Synthesizing these findings, burn severity in the Zagros is best understood as an outcome of multi-scale interactions. At the macro-scale, pre-fire moisture and surface temperature establish the primary constraints on watershed-level fire potential. Concurrently, at the local scale, fuel type and neighborhood structure, captured by spatial texture, govern fire propagation and intensification. Finally, contextual factors such as terrain and anthropogenic access gradients further modify fire exposure and spread potential. This integrated understanding explains why flexible, non-linear models that can learn complex interactions generalize far better than their linear counterparts in an environment characterized by a limited dynamic range and complex spatial controls ([47]).

Comparison with existing research

Our cross-validated model performance (XGBoost R² ≈ 0.70 and Random Forest R² ≈ 0.65) is consistent with results from other studies in structurally complex forests where modeling relies solely on static variables. Previous research has consistently affirmed the power of Random Forest for predicting burn severity ([17], [6]). As noted in literature reviews, a key challenge is that much of the unexplained variance in such static models do not stem from algorithmic limitations, but from the absence of dynamic fire-weather and fine-scale fuel moisture data ([30]).

Our finding that moisture (NDMI) and temperature (LST) are the most influential predictors aligns with other research highlighting the critical role of these factors. For instance, studies have shown that Land Surface Temperature (LST) is a foundational variable in developing burn severity indices ([60]). Similarly, a large-scale analysis in Canada’s mountain national parks found that variables related to vegetation moisture and temperature (like NDMI and LST) were among the most important drivers of burn severity (dNBR - [55]). This convergence underscores a widely applicable moisture-thermal mechanism.

This distinction resolves apparent contradictions in previous NDVI-burn severity studies. Current research on Normalized Difference Vegetation Index (NDVI) and burn severity relationships reveals seemingly contradictory findings that highlight this complexity. Studies have found that higher NDVI values immediately before fire events are associated with increased burn severity ([16]). However, pre-fire forest health as measured by NDVI also plays a crucial role in post-fire recovery patterns. Areas with higher pre-fire NDVI values showed significant increases in post-fire forest recovery metrics compared to areas with lower pre-fire NDVI values, suggesting that healthier pre-fire forest conditions may result in faster forest recovery ([53], [2]).

Our novel finding that a decadal negative NDVI trend enhances prediction is consistent with the principles of time-series analysis for forest monitoring. In this approach, long-term trends in satellite data, such as from Landsat, are used to detect and quantify changes in forest health and cover over time ([37]). The observed association between the decadal NDVI trend and high burn severity suggests that long-term vegetation decline may act as a slow-acting predisposing factor. Rather than directly driving fire intensity, this trend could reflect an ecosystem’s increasing vulnerability to severe fire over time.

Our models revealed that landscape structure, particularly vegetative land-cover types, plays a significant role in predicting fire severity. However, interpretations of land-cover influence must account for spectral artifacts; non-burnable surfaces (e.g., Bare/Water) can exhibit artificially high dNBR values due to edge effects, potentially confounding machine-learning outputs if not masked. The high importance of the land cover variable is consistent with findings that emphasize the need for vegetation-specific information to improve burn severity prediction ([15]). Furthermore, the higher severity observed in shrublands than in dense forests aligns with studies showing that severity varies significantly across vegetation layers ([36]). The importance of spatial texture has also been confirmed elsewhere; for example, Hultquist et al. ([17]) identified texture as a key predictor in their Random Forest model. While high dNBR values over bare surfaces are often artifacts of spectral edge effects ([7]), our conclusions regarding land-cover controls remain robust. By explicitly excluding these non-burnable surfaces from our core statistical evaluations (Tab. 5), we confirmed that structural differences among actual vegetative fuel types (e.g., shrublands vs. forests) genuinely drive the observed fire severity patterns. The high importance of non-burnable surfaces in the XGBoost model reflects the algorithm’s ability to maximize predictive accuracy by identifying “zero-severity” boundaries, although we excluded these areas from ecological mean-comparison tests to avoid biasing the drivers of actual fuel consumption.

A substantial portion of the residual uncertainty in our model likely stems from the absence of dynamic fire-weather data. While we do not directly test this in our study, the literature suggests these variables are critical, and studies that integrate them often report significant improvements in model accuracy. For instance, Vanderhoof et al. ([52]) increased their Random Forest model’s R² from 0.46 to 0.78 by adding event-specific weather indices like wind speed. This confirms that the primary limitation in predictive precision is often limited access to top-down forcings, a challenge also identified as a barrier to comparing results across studies ([30]).

Novelty, limitations, and future research directions

The primary novelty of this work lies not in the individual methodological components, but in their regional integration; we combined long-term NDVI trends with short-term biophysical variables and texture metrics to create a comprehensive, multi-scale framework for burn-severity modeling specific to the Zagros ecosystem. This was achieved by jointly integrating (i) long-term ecosystem decline (the decadal NDVI trend), (ii) pre-fire biophysical status (NDMI, LST), and (iii) a collinearity-screened GLCM texture metric. While machine learning has been used to model fire in various contexts, integrating long-term ecological degradation with short-term biophysical drivers remains largely unexplored in semi-arid Mediterranean-type forest ecosystems. Addressing this gap is critical because chronic vegetation decline, such as the oak decline observed in the Zagros, represents a slow-acting predisposing factor that fundamentally alters fuel characteristics and fire vulnerability. This study therefore establishes a reproducible baseline for future research on fire effects in these chronically stressed environments. These contributions, however, should be considered in light of two principal limitations: (i) missing dynamic fire-weather forcings - hourly data fields for wind, humidity, and atmospheric stability were unavailable for the fire window; such variables are widely acknowledged as dominant “top-down” drivers of burn severity, and their exclusion is a recognized source of unexplained variance in remote sensing-based severity studies ([30], [52]); (ii) the proxy nature of fuel descriptors - satellite-derived indices and texture metrics capture fuel moisture and structure only indirectly. They cannot resolve fine-scale dead-fuel moisture content or live-fuel loading, both of which are identified in the literature as common blind spots that propagate uncertainty into model outputs ([47]).

These constraints define clear priorities for future research: (i) integrate high-resolution meteorology - future models should incorporate data from station networks or downscaled reanalysis products to explicitly represent event-scale fire-weather coupling; (ii) incorporate field-calibrated fuel data - collecting or accessing field-calibrated fuel measurements is essential to validate and refine the spectral and textural proxies used, thereby addressing a critical data gap highlighted in comparative severity assessments; (iii) conduct external validation - the framework’s generalizability must be tested across multiple fires in the Zagros region; this includes benchmarking its performance against standardized burn-severity products, whose classification thresholds and ecological relevance continue to face scientific scrutiny.

Strategic implications for fire management in the Zagros region

The modeling results provide a foundation for shifting Zagros fire management from a predominantly reactive posture to a proactive, data-driven strategy that anticipates, rather than merely responds to, severe fire events.

Proactive vulnerability mapping

The strong pre-fire signals associated with NDMI, LST, and the decadal NDVI trend suggest the potential to identify landscapes primed for high burn severity months in advance. However, as our current findings are derived from a single fire event, prospective validation across multiple independent fires is required to confirm this as an operational early-warning capability. Operationally, these variables can be composited in Google Earth Engine to generate dynamic “fire-vulnerability hotspot” maps at the start of each fire season. This allows agencies to strategically stage suppression resources and focus public-awareness campaigns where risk is greatest ([12], [45]).

Landscape-structure management

The significance of both spatial texture metrics and topographic variables indicates that fire severity is governed not only by terrain but also by the spatial arrangement and continuity of fuels. Management, therefore, should extend beyond simple fuel-type inventories to actively manipulate spatial patterns, for example, by constructing fuel breaks across high-continuity corridors and conserving highly fragmented patches as potential refugia.

Resilience-oriented stewardship

While the predominance of low-to-moderate severity in the studied event could suggest underlying ecosystem resilience, alternative factors like favorable weather conditions during the fire or effective local suppression efforts must also be considered. By running the model in reverse and identifying combinations of moisture status, vegetation trend, and texture associated with low dNBR, managers can prioritize maintaining landscape attributes that foster desirable, low-severity fire regimes. Such a shift aligns with emerging regional policy recommendations that call for resilience-focused planning rather than an exclusive reliance on full suppression.

Collectively, these implications argue for embedding the presented modeling framework within routine decision cycles, leveraging open-source platforms to deliver near-real-time risk products and guiding interventions that shape, not simply extinguish, future fire behavior.

  Conclusions 

This study aimed to develop a predictive model for burn severity within the complex, fire-prone ecosystems of the Zagros Mountains by integrating long-term ecosystem health trends, spatial texture metrics, and traditional biophysical variables. Our findings indicate that machine learning ensembles, particularly XGBoost, offer a robust predictive framework. The models revealed that pre-fire biophysical conditions (notably moisture and temperature) and fuel type (land cover) are strongly associated with fire effects. Furthermore, long-term vegetation decline (NDVI trend), topography, and landscape structural heterogeneity (texture) emerged as important secondary predictors of local burn severity. However, these findings must be interpreted in the context of several methodological limitations. Given the correlational design, the identified relationships represent plausible ecological associations rather than proven causal mechanisms. Additionally, the reported predictive performance likely reflects upper-bound estimates due to spatial autocorrelation, and accurate interpretation required the careful masking of spectral artifacts over non-burnable surfaces. Because this study is based on a single fire event, the current framework provides a foundational, data-driven understanding of regional vulnerability rather than an operationally ready management tool. Future research should prioritize prospective validation across multiple independent fires and integrate dynamic fire-weather data. Ultimately, this study advances the scientific understanding of spatial fire dynamics in the Zagros region, offering a vital first step toward developing proactive, spatially explicit landscape management strategies.

  List of abbreviations 

The following abbreviations have been used throughout the text:

• dNBR: differenced Normalized Burn Ratio
• GEE: Google Earth Engine
• GLCM: Gray-Level Co-occurrence Matrix
• HC3: Heteroscedasticity-Consistent Standard Errors (type 3)
• LST: Land Surface Temperature
• NDMI: Normalized Difference Moisture Index
• NDVI: Normalized Difference Vegetation Index
• NIR: Near-Infrared
• OLS: Ordinary Least Squares
• RF: Random Forest
• RMSECV: cross-validated Root Mean Square Error
• R²CV: cross-validated R-squared
• VIF: Variance Inflation Factor
• XGB: XGBoost

  Acknowledgements 

We thank the University of Tehran for providing access to research facilities. No external funding was received for this research.

  References

  Supplementary Material

Appendix 1 - Source code for preprocessing and burn-severity mapping.
Tab. S1 - Multicollinearity diagnostics for OLS (HC3): variance inflation factors (VIF) for all predictors.