Elevation Outlasts Fire Severity as the Dominant Driver of Post-Fire Vegetation Recovery in Kosciuszko National Park

1. Introduction

The 2019–2020 Australian bushfire season, popularly known as the Black Summer, burned an estimated 5.4 million hectares of New South Wales (NSW), including 2.7 million hectares of national park estate, and was the most extensive fire season in the state's recorded history [1]. The scale of the event has been attributed to the compounding effects of record-breaking drought, extreme fuel dryness, and anomalously high spring and summer temperatures across southeastern Australia [2]. Ecological impacts were correspondingly severe, with long-unburnt fire-sensitive vegetation communities, including Gondwanan rainforest relicts, burning for the first time in recorded history in some locations. Understanding how vegetation recovers from disturbance of this magnitude, and what governs the pace and completeness of that recovery, is therefore central not only to post-fire land management and biodiversity conservation in the immediate term, but to anticipating ecosystem response to a fire regime that is projected to intensify in frequency and severity under continued climate change.

Remote sensing time series of spectral vegetation indices, most commonly the Normalized Difference Vegetation Index (NDVI) and the Normalized Burn Ratio (NBR), have become the standard tool for monitoring post-fire recovery at landscape scale, offering a spatially continuous and temporally repeatable alternative to labour-intensive field survey [3]. The advent of cloud-based geospatial platforms such as Google Earth Engine [4] has made it computationally feasible to construct multi-year, cloud-free composite time series across large or remote study areas without local data storage or processing infrastructure and has substantially lowered the barrier to entry for this type of analysis. Combined with Random Forest regression [5] and, increasingly, SHapley Additive exPlanations (SHAP) [6] for interpreting the relative contribution of individual predictors, this workflow has become close to a standard template, and has been applied to fire events across boreal, Mediterranean, and temperate forest systems worldwide [7].

This template has recently been applied to Australian fires specifically. A 2026 study of the 2009 Black Saturday fires in Victoria used an equivalent NDVI/NBR/Random Forest framework to show that canopy greenness (NDVI) recovered within approximately 8 years while structural recovery (NBR) lagged, requiring around 11 years in high-severity zones [8]. A separate analysis of the Gospers Mountain fire, part of the same Black Summer event, used the TSS-RESTREND method to show recovery rate varied by both fire severity and land cover, and that above-average post-fire rainfall accelerated recovery by an estimated two years [9]. At a global scale, an analysis of post-fire productivity recovery found post-fire climatic conditions, particularly soil moisture, vapour pressure deficit and air temperature, dominated recovery time across roughly 89% of assessed burned area worldwide, more so than pre-fire productivity or fire-induced productivity loss [10], a pattern that, if it holds locally, predicts the severity-to-climate handover this study tests directly.

A second gap this study addresses is methodological. Random Forest models of post-fire recovery are almost universally validated using random train-test splitting, yet recovery values are spatially autocorrelated by construction, since neighbouring pixels experience correlated fire behaviour, terrain and climate. Random cross-validation under these conditions systematically underestimates true predictive error, because geographically close training and test points leak information across the split even with no spatial term in the model; block-structured cross-validation, holding out spatially contiguous groups together, is recommended wherever such dependence exists [11]. None of the studies above [7,8,9,10] report spatial cross-validation, leaving open how much of their reported skill reflects genuine, transferable relationships versus spatial autocorrelation exploited by random splitting.

This study therefore has four objectives: (1) to quantify vegetation recovery trajectories across fire severity classes in Kosciuszko National Park from 2021 to 2025; (2) to use SHAP analysis of Random Forest models to test whether the relative importance of fire severity, terrain, climate, soil, and vegetation-community predictors shifts over the first five years of recovery; (3) to test, using both random and spatial block cross-validation across a range of block resolutions, whether a comprehensive twelve-variable predictor set can meaningfully predict recovery at unseen locations, and to situate this result within the broader methodological literature on spatial cross-validation in ecology; and (4) to test explicitly whether elevation zone, as a proxy for vegetation type and growing-season length, produces divergent recovery trajectories distinct from, and more temporally persistent than, the effect of fire severity alone.

2. Materials and Methods

2.1. Study Area

During the 2019–2020 Black Summer fire season, the Dunns Road fire and the Adaminaby fire complex (originating from the Green Valley fire) burned through the northern and central portions of the park, consuming approximately 231,000 hectares, or 33.5% of the park's total area (Table 1) [12,13]. Burnt areas were concentrated in the montane and subalpine zones; the highest alpine areas near the park's southern boundary, including the Thredbo and Perisher resort areas and the summit region, were largely unaffected. This uneven distribution of burnt area across the elevation gradient is a feature of the fire event itself rather than a data-processing artefact and is reflected in the spatial coverage of the fire severity data used in this study (Section 2.2), shown in Figure 1.

2.2. Fire Severity and Fire History

Fire severity for the 2019–2020 fire season was obtained from the New South Wales Fire Extent and Severity Mapping (FESM) version 3 product, a statewide 10 m raster derived from Sentinel-2 imagery using a Random Forest classifier trained on aerial photo interpretation, with five severity classes (unburnt, low, moderate, high, extreme) plus an experimental reserved class not used operationally at the time of this study [1]. Independent aerial-photo-interpretation validation yielded balanced accuracies of 0.76–0.90 across the low-to-extreme severity classes, while post-fire field validation yielded 0.73–0.84 for the low-to-high classes; field validation for the extreme class was not reported [1]. The statewide raster was clipped to the park boundary, obtained from the NSW National Parks and Wildlife Service (NPWS) Estate dataset. Fire history back to the early twentieth century was obtained from the NSW NPWS Fire History dataset and clipped to the same boundary, and the number of previously recorded fires at each sample location was used as a predictor of fire-return-interval effects. The spatial distribution of fire severity across the park is shown in Figure 2.

2.3. Satellite Vegetation Indices and the Recovery Index

Sentinel-2 Level-2A surface reflectance imagery [14] was accessed via Google Earth Engine. For each year from 2019 to 2025, a cloud-masked median composite was constructed from the January–March period, chosen to hold seasonal phenology constant and avoid conflating seasonal variation with fire-driven change. From each composite, NDVI [15], NBR, the Normalized Difference Moisture Index (NDMI), and the Enhanced Vegetation Index (EVI) were calculated; NDVI was used as the primary index for the recovery calculation reported here, with NBR, NDMI, and EVI retained as ancillary indices for future cross-validation of the recovery signal.

A pixel-level recovery index was calculated for each year t from 2021 to 2025 as:

Recovery_t = (NDVI_t − NDVI_2020) / (NDVI_2019 − NDVI_2020)

where NDVI_2019 represents pre-fire vegetation condition and NDVI_2020 represents immediate post-fire condition. A value of 1.0 indicates full recovery to pre-fire NDVI; values above 1.0 indicate NDVI exceeding the pre-fire baseline. Because this ratio is unstable where the pre/post-fire NDVI difference is small, calculation was restricted to pixels classified as burnt (FESM severity class 2–5) and extreme outlier values were clipped to the range [−3, 3].

2.4. Environmental Predictors

Terrain variables (elevation, slope, aspect) were derived from the Copernicus GLO-30 digital elevation model (30 m). Climate variables (cumulative rainfall and mean temperature over 2020–2021, 2020–2023, and 2020–2025) were obtained from the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) [16] and ERA5-Land [17], which also provided volumetric soil moisture for the uppermost layer (0–7 cm) over the same windows. Static soil properties (organic carbon, clay, sand, pH, all at 0–5 cm depth) were obtained from SoilGrids 2.0 at 250 m resolution [18]. Generalised land cover came from the European Space Agency (ESA) WorldCover v200 (10 m) [19], and fine-scale vegetation community type from the NSW State Vegetation Type Map (SVTM, Extant edition; ~100 Plant Community Types (PCT) at 5 m resolution). All raster layers were accessed and pre-processed via Google Earth Engine [4] before extraction to point locations. Elevation and slope, the two topographic variables identified as most influential in subsequent analyses (Section 3.5), are mapped in Figure 3.

2.5. Sampling Design and Dataset Construction

5,000 sample points were generated within the burnt footprint, stratified to include 1,250 points from each of the four FESM severity classes (low, moderate, high, extreme). At each point, the recovery index for each year (2021–2025), pre-fire NDVI, and all environmental predictors listed above were extracted. After removing points with missing values arising from edge effects at raster boundaries, 4,864–4,907 points remained for analysis, depending on the specific variable set used.

2.6. Spatial Autocorrelation

Global spatial autocorrelation in the 2021 recovery values was tested using Moran's I [20] with a k-nearest-neighbour (k = 8) row-standardised spatial weights matrix and 999 permutations. Analyses used Python 3.8.20 with the libpysal (v4.8.1) and esda (v2.5.1) packages.

2.7. Random Forest Modelling and Cross-Validation

Three Random Forest regression models (300 trees, scikit-learn v1.3.2 implementation in Python 3.8.20 [21]), a widely used ensemble method for remote sensing and environmental prediction tasks [22], were built, predicting the recovery index for 2021, 2023, and 2025, using severity, pre-fire NDVI, elevation, slope, aspect, land cover, vegetation community (SVTM PCT), fire history count, soil organic carbon, clay, sand, pH, and the rainfall, temperature, and soil moisture variables for that year's recovery window. Each model was evaluated two ways, reporting both R² and root mean square error (RMSE): (i) a conventional 80/20 random train-test split, and (ii) five-fold spatial block cross-validation, in which sample points were assigned to spatial blocks using k-means clustering on point coordinates, and folds were built from held-out blocks rather than held-out points, reducing spatial dependence between training and test data relative to random splitting, following the block cross-validation approach recommended for spatially structured ecological data [11] and implemented in dedicated software such as the blockCV package [23]. To test sensitivity to the arbitrary choice of block count, the 2025 model was additionally evaluated using 5, 10, 15, 20, 30, and 50 spatial blocks, spanning coarse to fine partitions of the study area.

2.8. SHAP Analysis and Uncertainty Quantification

For each of the three fitted Random Forest models, SHAP (SHapley Additive exPlanations) values were used to quantify each predictor's contribution to model output [6], computed using the TreeExplainer implementation in the shap package (v0.44.1). The mean absolute SHAP value for each predictor was used as a measure of overall variable importance, allowing direct comparison of variable importance rankings across the three modelled years. To quantify uncertainty around these importance estimates, a bootstrap procedure was applied: for each modelled year, 20 bootstrap resamples of 2,000 points were drawn with replacement, a Random Forest (100 trees) was refitted on each resample, and mean absolute SHAP values were recomputed using shap package v0.52.0 (run in a separate computing environment from the primary analysis; numpy v2.5.2, scipy v1.18.1); 95% confidence intervals were taken as the 2.5th and 97.5th percentiles of the resulting bootstrap distribution for each variable.

2.9. Cross-Index Validation with NBR

To evaluate whether the NDVI-based recovery index used throughout this study is sensitive to the choice of spectral index, an equivalent recovery index was calculated using the Normalized Burn Ratio (NBR) [24], following the identical formula, masking, and clipping procedure described in Section 2.3, using NBR calculated from the same Sentinel-2 composites. Pearson and Spearman correlations between NDVI-based and NBR-based recovery were calculated at each of the 4,864 sample points for all five recovery years, and a paired t-test was used to test for systematic directional differences between the two indices.

2.10. Hotspot Analysis

Spatial clusters of significantly fast (“hotspot”) and slow (“coldspot”) recovery in 2025 were identified using the Getis-Ord Gi* local statistic [25] (esda v2.5.1 implementation), computed on the same k = 8 nearest-neighbour, row-standardised spatial weights matrix used for the Moran's I test, with 999 conditional permutations. Self-inclusion for Gi* used esda's default convention, with each point's self-weight set equal to its largest neighbour weight. Points with a permutation P-value below 0.05 were classified as significant hotspots (positive z-score) or coldspots (negative z-score); all other points were classified as not significant. Mean elevation was compared across the three classes as an independent test of the elevation effect identified through SHAP analysis, using both Welch's t-test and the Mann-Whitney U test to confirm the result was not dependent on distributional assumptions.

3. Results

3.1. Recovery Trajectories by Fire Severity Class

Mean recovery differed clearly among severity classes in the first-year post-fire, with low-severity areas closest to pre-fire NDVI (1.034) and high-severity areas furthest below it (0.722; Table 2, Figure 4). All classes rose through 2022–2023, with low- and moderate-severity areas briefly exceeding baseline. From 2023 onward, however, the ordering among severity classes broke down: by 2025, all four classes had converged to a narrow band between 0.88 and 0.92, indicating that the initial severity gradient in recovery had largely disappeared within five years.

3.2. Predictor Screening

Variance inflation factors (VIF) for the twelve-variable predictor set were all below 5, indicating no severe multicollinearity (Table 3, Figure 5). The highest VIF (4.66, clay content) reflects an expected soil-texture relationship with sand content but does not warrant exclusion of either variable.

3.3. Spatial Autocorrelation

The 2021 recovery values exhibited statistically significant positive spatial autocorrelation (Moran's I = 0.180; expected I under the null hypothesis = −0.0002; P = 0.001; z = 27.09), confirming that recovery values are not spatially independent and that random train-test splitting is likely to overstate model performance.

3.4. Random Forest Performance: Random Split Versus Spatial Cross-Validation

Random-split R² ranged from 0.100 to 0.161 across the three modelled years, indicating weak-to-moderate apparent predictive skill. Under spatial block cross-validation, however, R² fell to near zero or negative in every year (Table 4, Figure 6), indicating that the twelve-variable predictor set did not generalise to unseen locations despite showing weak positive skill under random splitting.

To test whether this gap reflected an arbitrary choice of spatial block resolution, the 2025 model was re-evaluated across six block counts (5–50). Spatial R² remained consistently near zero or negative at every resolution tested (range: −0.146 to +0.003), well below the random-split R² of 0.100 at every point (Figure 7), indicating the gap is a robust feature of the data rather than an artefact of one arbitrary parameter choice.

3.5. SHAP Variable Importance and Its Temporal Shift

SHAP analysis of the three fitted models revealed a marked temporal shift in variable importance (Table 5, Figure 8). Fire severity was the second most important predictor at year 1 (mean |SHAP| = 0.116) but declined to among the least important predictors by year 5 (0.019), a five-fold reduction that closely mirrors the convergence of severity-class recovery trajectories described in Section 3.1. Elevation, by contrast, was the most important predictor at year 1 (0.118) and year 5 (0.106), dipping only at year 3 (0.071), making it the most temporally stable predictor in the set. Pre-fire NDVI showed a similar pattern to elevation, peaking at year 3 (0.116). Rainfall, temperature, and soil moisture for the corresponding recovery year each showed modest importance (0.02–0.03) in the year they were most relevant, consistent with a growing role for climate/moisture variables as the influence of fire severity recedes.

The direction, not only the magnitude, of each variable's effect appears in the SHAP summary plots for each year (Figure 9, Figure 10 and Figure 11). In 2021 (Figure 9), high severity (red) is associated with negative SHAP values, confirming severely burned sites recovered more slowly, while low severity (blue) pushes predictions upward; this separation weakened by 2023 (Figure 10) and largely collapsed by 2025 (Figure 11), consistent with the declining severity importance reported above. Elevation shows a more persistent pattern: low elevation (blue) is associated with positive SHAP values in every year, while high elevation skews negative, consistent with the elevation-zone divergence in Section 3.7. Pre-fire NDVI shows the clearest direction: low pre-fire NDVI is consistently associated with negative SHAP values, indicating sparser pre-fire vegetation predicted less complete recovery.

To confirm that these year-to-year differences reflect genuine shifts rather than sampling noise, a bootstrap procedure (20 resamples of 2,000 points per year, 100 trees per resample) was used to compute 95% confidence intervals around each mean |SHAP| estimate (Table 6, Figure 12). Severity's confidence interval in 2021 (0.104–0.135) does not overlap at all with its 2025 interval (0.014–0.029), confirming that the decline in severity's importance is a robust feature of the data rather than an artefact of a single model fit. Elevation's confidence intervals show a similar pattern to Table 5's point estimates: elevation is significantly lower at 2023 (0.057–0.095) than at 2021 (0.108–0.144), but its 2021 and 2025 (0.072–0.130) intervals overlap, indicating that while elevation's dip at year 3 is statistically real, the claim that its importance is fully “unchanged” between year 1 and year 5 should be read as broadly stable rather than statistically identical.

3.6. Cross-Validation with NBR-Based Recovery

To test whether the NDVI-based recovery index is robust to the choice of spectral index, an equivalent recovery index was calculated using NBR (Section 2.9) and compared against the NDVI-based index at all sample points for all five recovery years (Table 7, Figure 13). NDVI-based and NBR-based recovery were moderately to strongly correlated in rank order (Spearman r = 0.70–0.77 across years) but weaklier correlated in absolute magnitude (Pearson r = 0.45–0.53), indicating that the two indices agree well on which locations recovered better than others but disagree more on the precise degree of recovery. A paired t-test at each year found that NDVI-based recovery was significantly higher than NBR-based recovery in every single year (mean difference 0.03–0.14, all P < 0.001; Table 7), consistent with the canopy-structure decoupling reported for the 2009 Black Saturday fires in Victoria, where canopy greenness (NDVI) recovered markedly faster than structural condition (NBR) [8]. The gap between the two indices was largest in 2021 and 2024 (mean difference 0.12–0.14) and smallest in 2025 (0.03), suggesting this canopy-structure decoupling is most pronounced early in recovery and narrows, though does not fully close, as structural regrowth gradually catches up to canopy greening.

3.7. Recovery Trajectories by Elevation Zone

Stratifying sample points into montane forest (<1,400 m) and subalpine woodland (1,400–1,850 m) revealed that the two zones followed broadly similar recovery trajectories through 2023 before diverging sharply (Table 8, Figure 14). No burnt sample points existed above 1,850 m, consistent with the fire's concentration in the park's montane and subalpine, rather than true alpine, zone. The 95% confidence intervals around each zone's yearly mean do not overlap at any point across 2021–2025, and the gap widens rather than narrows after 2023, indicating the two zones are statistically distinguishable throughout, not only in later stages. Subalpine woodland recovery declined from a 2023 peak of 0.957 to 0.769 by 2025, falling below its own 2021 starting value and well below full recovery. Montane forest, by contrast, declined more gently over the same period (1.020 to 0.927), remaining close to the pre-fire NDVI baseline throughout. The underlying spatial pattern of 2025 recovery, from which this elevation-zone divergence emerges, is mapped in Figure 15.

3.8. Hotspot Analysis

Getis-Ord Gi* analysis of 2025 recovery values identified 469 points (9.6%) as significant fast-recovery hotspots and 311 points (6.3%) as significant slow-recovery coldspots, with the remaining 4,084 points (83.2%) showing no significant local clustering. Mean elevation differed substantially and in the expected direction across these three classes (Table 9, Figure 16): hotspots occurred at a mean elevation of 877 m (s.d. 242), compared to 1,284 m (s.d. 275) for coldspots. This difference was confirmed as highly statistically significant by both a parametric and a non-parametric test (Welch's t = −21.20, P = 4.6×10⁻⁷⁵; Mann–Whitney U = 19,966, P = 3.2×10⁻⁶⁶), providing an independent confirmation, via a distinct spatial statistical method, of the elevation effect identified through SHAP analysis. The spatial distribution of these hotspots and coldspot clusters across the park is shown in Figure 17.

4. Discussion

4.1. Fire Severity is a Transient, not Persistent, Driver of Recovery

The clearest and most consistent finding of this study is that fire severity's influence on recovery in Kosciuszko National Park is transient rather than persistent. Severity strongly separated recovery trajectories at one-year post-fire, in both the descriptive comparison (Table 2) and SHAP-based importance (Table 5), with a roughly five-fold decline between 2021 (0.116) and 2025 (0.019). This is consistent with the Gospers Mountain TSS-RESTREND analysis, which similarly found severity-dependent differences that narrowed as post-fire rainfall accumulated [9], and with the global analysis of Xu et al. (2024), which found post-fire climatic conditions dominated recovery time across most global burned area, exceeding pre-fire productivity and fire-induced productivity loss once the acute phase has passed. Our results extend this by quantifying the shift within a single modelling framework, and showing the crossover from severity-dominated to climate/moisture-influenced regime occurs within a narrow window, between roughly the first and third year post-fire, for the communities sampled here.

4.2. Elevation as a Persistent, and Largely Independent, Driver

In contrast to severity, elevation was the most temporally stable predictor of recovery across the five-year study window, remaining the single most important SHAP predictor in both 2021 and 2025 and dipping only in the intervening 2023 model. This effect was corroborated independently through spatial hotspot analysis, a fundamentally different statistical approach based on local spatial clustering rather than variable-importance decomposition: fast-recovery clusters sat roughly 400 m lower in elevation, on average, than slow-recovery clusters (877 m versus 1,284 m; Table 9). The convergence of two independent analytical methods on the same conclusion substantially strengthens confidence that this is a genuine ecological signal rather than a modelling artefact specific to Random Forest or SHAP.

The persistence of pre-fire NDVI as an important predictor, particularly at year 3 (SHAP = 0.116, the highest value recorded for any variable at any year), is consistent with ecological memory, whereby the legacy of pre-disturbance ecosystem state shapes post-disturbance trajectories independently of the disturbance event itself [26]. Under this interpretation, sites with denser, more vigorous pre-fire vegetation retain greater capacity for recovery: at lower fire severity, eucalypts commonly resprout epicormically from surviving stems [27], whereas more severely burnt stems draw on root and lignotuber reserves and local seed bank as an information and material legacy, consistent with lignotuber-based vegetative recovery documented in Kosciuszko snow gum one year after the 2003 fires [28]. Sites with greater pre-fire vegetation density and vigour therefore recover more readily and rapidly once the initial post-fire growth phase is underway. This is consistent with the broader concept of ecological memory introduced above, in which antecedent ecosystem state measurably, and sometimes under appreciably, shapes disturbance recovery trajectories.

4.3. Limited Spatial Predictability, and Its Implications for the Field

The most methodologically important finding of this study, however, is the consistent failure of a comprehensive twelve-variable predictor set, spanning fire severity and history, terrain, climate, soil physicochemical properties, and fine-scale vegetation community type, to generalise under spatial cross-validation, regardless of the spatial block resolution used to define the cross-validation folds (Section 3.4, Figure 7). This has at least two distinct implications, one ecological and one methodological.

Ecologically, it suggests that recovery at any location retains a meaningful component not explained by any combination of the variables tested, even though several (severity, elevation, pre-fire NDVI) show clear relationships with recovery in aggregate: these variables shift the distribution of likely outcomes without narrowing it enough for confident pixel-level prediction. The unexplained variance may reflect fine-scale processes below the resolution of the data used here, including microtopographic variation in soil moisture and frost exposure, localised seed dispersal and recruitment limitation, patchy herbivory, or individual weather events during the critical early resprouting window. Feral horse grazing and trampling is a plausible source of this patchy variance: feral horses cause extensive, spatially uneven damage to vegetation, soil and waterways across Kosciuszko and the broader Australian Alps, concentrated in particular sub-catchments rather than distributed evenly [29]. A paired comparison of Sphagnum peatland sites within the park found that sites with feral horse presence were net carbon sources while similar sites without horses were net carbon sinks [30], a site-to-site difference driven by an unmapped biotic disturbance that none of the twelve landscape-scale predictors used here were designed to capture.

Methodologically, this result is a caution for the broader post-fire recovery literature. Block cross-validation, or an equivalent method respecting spatial dependence, should be used wherever such structure exists, even where no correlation is visible in fitted model residuals [11]. This pattern is not unique to post-fire recovery: an analogous demonstration in large-scale tropical forest biomass mapping found a model explaining over half of variance under random validation showed quasi-null predictive power once spatial autocorrelation was accounted for [31]. None of the directly comparable Australian or global post-fire recovery studies discussed here [7,8,9,10] report spatial cross-validation alongside random-split metrics, despite recovery values being spatially autocorrelated by construction in every case, a property we confirmed for our own study area via Moran's I (Section 3.3). The gap identified here, roughly 0.10–0.25 R² units, and its consistency across a six-fold range of block counts (Figure 7), suggests that predictive skill based solely on random splitting may be optimistic when strong spatial dependence is present, as it likely is in comparable published studies that do not report spatial cross-validation. We recommend that spatial cross-validation, or at minimum an explicit report of spatial autocorrelation, become standard practice when Random Forest methods are applied to spatially continuous ecological outcomes.

4.4. Implications for Post-Fire Land Management

From a management perspective, these results offer two practical points. First, the rapid convergence of severity-class recovery by year 3–5 suggests initial burn severity is a comparatively poor long-term indicator of which locations need sustained management attention; sites should instead be re-assessed periodically, since a site severely affected in year 1 may show no meaningful difference from a lightly burned site by year 3. Second, the persistent elevation effect and late-stage subalpine decline suggest post-fire monitoring and intervention (such as feral herbivore control or supplementary planting) may be more usefully prioritised toward higher-elevation burned areas from around year 3 onward, when the divergence from lower-elevation forest becomes apparent, rather than allocated primarily on initial severity mapping.

4.5. Limitations and Future Directions

This study has several limitations. First, the FESM severity product covers approximately half the park's area; the unmapped remainder was unburnt higher-elevation terrain, meaning recovery could not be assessed above ~1,850 m. Second, this is a single-park, remote-sensing-only study without direct field validation; related field-based studies exist for this park for other vegetation types [33,34], and we treat the present findings as hypothesis-generating rather than confirmatory, offering specific testable predictions such as field-measured resprouting vigour decoupling from canopy NDVI recovery in high-severity subalpine plots by year 3–5. Third, NDVI systematically overstated recovery relative to NBR in every year (Section 3.6), consistent with known NDVI saturation in dense canopy [8,32]; NDMI and EVI composites remain available for further cross-index validation. Finally, the five-year study window captures only the early phase of what is likely a multi-decadal recovery process for slower-growing subalpine communities, and continued monitoring will be needed to confirm whether the late-stage subalpine decline identified here persists or eventually re-converges with montane forest.

5. Conclusions

Vegetation recovery following the 2019–2020 Black Summer fires in Kosciuszko National Park shows a clear temporal shift in its dominant drivers: fire severity strongly shapes recovery in the first year but fades to near-negligible within five years, a decline confirmed as robust by bootstrap resampling, while elevation remains the most important predictor throughout, with subalpine snow gum woodland showing a statistically distinct decline across the entire study period not observed in lower-elevation montane forest. A cross-index comparison further showed that NDVI systematically overstates recovery relative to NBR every year, indicating canopy greening and structural recovery decouple over time. At the same time, a comprehensive twelve-variable predictor set failed to generalise to unseen locations under spatial cross-validation, robustly across a wide range of block resolutions, despite weak apparent skill under random splitting. Together these results suggest post-fire recovery in this alpine system is governed by a rapidly fading severity effect, a persistent elevation-mediated vegetation-type effect compounded by independently documented climate [35], dieback [36] and reburn stressors specific to this landscape, a canopy-structure decoupling detectable across spectral indices, and a substantial residual component that static covariates cannot resolve. We recommend that future Random Forest-based post-fire recovery studies routinely report spatial, in addition to random, cross-validation results and bootstrap-based uncertainty, and suggest alpine and subalpine systems, several of which are already classified as threatened by climate- and fire-regime-driven risk assessments [37], warrant further targeted investigation given the distinctive late-stage recovery divergence identified here.