Burn Severity and Post‐Fire Weather Are Key to Predicting Time‐To‐Recover From Australian Forest Fires

Climate change has accelerated the frequency of catastrophic wildfires; however, the drivers that control the time‐to‐recover of forests are poorly understood. We integrated remotely sensed data, climate records, and landscape features to identify the causes of variability in the time‐to‐recover of canopy leaf area in southeast Australian eucalypt forests. Approximately 97% of all observed burns between 2001 and 2014 recovered to a pre‐fire leaf area index (±0.25 sd) within six years. Time‐to‐recover was highly variable within individual wildfires (ranging between ≤1 and ≥5 years), across burn seasons (90% longer January to September), and year of fire (median time‐to‐recover varying four‐fold across fire years). We used the logistic growth function to estimate the leaf area recovery rate, burn severity, and the long‐term carrying capacity of leaf area. Time‐to‐recover was most correlated with the leaf area recovery rate. The leaf area recovery rate was largest in areas that experienced high burn severity, and smallest in areas of intermediate to low burn severity. The leaf area recovery rate was also strongly accelerated by anomalously high post‐fire precipitation, and delayed by post‐fire drought. Finally we developed a predictive machine‐learning model of time‐to‐recover (R2: 0.68). Despite the exceptionally high burn severity of the 2019–2020 Australian megafires, we forecast the time‐to‐recover to be only 15% longer than the average of previous fire years.


Introduction
The remarkable magnitude and scale of the Australian 2019-2020 Black Summer megafires have been well documented (Boer et al., 2020;Nolan et al., 2020) and have highlighted the importance of wildfire prediction and forest management.Even prior to the megafires, the Southeast Australian fire regime was marked by comparatively large and intense wildfires (Archibald et al., 2013).On the current path of carbon emissions, climate change will force more forests into a state of post-fire recovery as the frequency of fire-weather increases (Abatzoglou et al., 2018;Abram et al., 2021;Boer et al., 2017;Li et al., 2021).Increased frequency is particularly concerning for native forest regeneration (Gallagher et al., 2022).While research has established the influence of fire-weather (Clarke et al., 2011) and forest structure on forest fire occurrence (Cawson et al., 2020;Sullivan et al., 2012;Zylstra et al., 2016), less research has investigated the drivers and timescales of forest-level recovery from fire in Australia.This gap in our understanding of the drivers and timescales of post-fire recovery limits our capacity to build predictive ecosystem forecasts (Dietze et al., 2018) that could be operationalized to aid ecosystem management.
Natural resource management and disaster preparedness stand to benefit from readily updateable quantitative frameworks that can produce ecological forecasts of forest time-to-recover (TTR) from fire.This is a crucial knowledge gap for developing bushfire preparedness in fire adapted but flammable woody ecosystems, and for planning the management of individual species or ecological communities requiring conservation actions (Kelly et al., 2015).For example, we have limited capacity to predict the fate of broad suites of threatened tree species which are impacted by fires across their geographic range (Gallagher, Allen, et al., 2021).Changing fire regimes will place more species at risk of extinction, and effective management relies on species-specific knowledge of recovery across the geographic range (Bradstock et al., 1998).Observations of TTR from fire for plant species is typically gained from repeat field surveys, however remote sensing has emerged as a viable option to complement field surveys, particularly in inaccessible terrain.Further, it is only with the sample size of the remote sensing record that we can hope to identify generalizable patterns of TTR across the backdrop of southeast Australia's high year-to-year precipitation variability.Ultimately, without a robust understanding of TTR from disturbance, it is challenging to develop the theory to underpin any quantitative predictions that could aid the planning of management interventions under continued climate warming and climate extremes.
A framework to characterize forest recovery from disturbance is needed, where "recovery" can be defined as when a forest approaches a baseline (Ingrisch & Bahn, 2018) such as the long term mean annual value of a forest property (e.g., leaf area index; LAI).A forest's TTR can be defined as the interval between a fire and the point where the forest property of interest has achieved and stabilized near the baseline state.Standard field-measured forest ecology metrics such as standing biomass or species diversity can help characterize the multi-dimensional nature of forest recovery, but are difficult to gather in remote locations, and are seldom measured at broad spatial scales before and after fire.Satellite remote sensing data provides the only practical means to systematically observe both pre-and post-fire vegetation in order to examine TTR through time, and across the large and rugged range of southeast Australia's forests.Here we focus on the recovery of forest canopy leaf area, in part because it is amenable to satellite based remote sensing but also because it represents a critical dimension of forest recovery of relevance to regional water, energy, and carbon budgets.
TTR is an emergent property that is dependent on three sub-components: burn severity, the recovery rate, and the threshold defined as recovery.Burn severity in this context is the fraction of a forest property (e.g., leaf area, biomass, species richness) removed by the fire relative to its pre-fire condition, and is linked to the accumulation of understory fuels and anomalously dry and windy meteorological conditions (Bradstock et al., 2010), but is also mediated by species composition (Whitman et al., 2018) and forest structure (e.g., tree density, basal area; Cocke et al., 2005).Evidence is mixed on the impact of burn severity on eucalypt TTR, and likely dependent on both spatial scale of evaluation and forest type (Hislop et al., 2019).Crown fires in particular produce vastly different recovery trajectories in forest vegetation, where the microclimate buffering effects from the closed canopy are lost (Collins, Hunter, et al., 2021).Hot-burning fires may also slow recovery by incinerating the seedbank and surficial roots, the latter leading to erosion of topsoil if heavy rains follow fire (Auld et al., 2020).
The post-fire recovery rate is complicated to assess, frequently nonlinear, and dependent on many factors related to life-history strategy, some of which are still poorly understood (Fairman et al., 2015).Recovery rate is also dependent on intrinsic life-history strategies, and the disturbance severity caused by the fire.The ecological strategy governing fire regeneration sets the starting point for a burnt forest on its trajectory toward recovery.Unsurprisingly, obligate reseeding species whose mature individuals are killed by fire, have the longest TTR that can span to multiple decades (Bowman et al., 2016).Some species can withstand some understory fire but cannot resprout following a crown fire (e.g., E. regnans).However most eucalypts have the capacity to resprout following fire from a basal storage organ or from the epicormic tissue along trunks and stems (Burrows, 2013;Clarke et al., 2015).Resprouting species should be capable of recovering biomass faster than obligate reseeding species, and therefore maintain competitive ecological advantage in forests.However recovery via resprouting is contingent on sufficient supply of stored nonstructural carbohydrates to subsidize the reconstruction of the canopy, and the supporting root and sapwood structures (Pausas & Keeley, 2017).This is important to note because pre-fire drought conditions may reduce the pool of nonstructural carbohydrates (He et al., 2020;Mitchell et al., 2013), thereby reducing the stored resources available for post-fire growth.
In this study we used the Eucalyptus-dominated forests and woodlands of southeast Australia as a model system to examine the interactions between climate, burn severity, and species differences on forest TTR.Southeast Australia's woody ecosystems assemble along a moisture gradient spanning semi-arid woodlands in the interior to mesic temperate and subtropical forests along the coast.Southeast Australian vegetation also endures unusually high interannual climate variability (Jiang et al., 2017).Recent decades have witnessed the two most severe droughts in the observational record (2001-2009;van Dijk et al., 2013;and2017-2019;CSIRO & Bureau of Meteorology, 2020).However, both droughts were characteristically broken by La Niña events (King et al., 2020) which delivered exceptionally high rainfall in the years following the droughts.Australian forests are also highly distinct in species composition.Eucalypts (including the genera Angophora, Corymbia, and Eucalyptus) are the dominant tree species across 74% of Australia's forests (Australian Department of Agriculture, 2014), and are amongst the world's most fire adapted ecosystems with high rates of burning both historically (Crisp et al., 2011) and in the modern day (Miller & Murphy, 2017;van der Velde et al., 2021).The majority of eucalypt tree species resprout after fire (Burrows, 2013;Nicolle, 2006) owing in part to epicormic buds embedded underneath the bark or basal lignotubers, though obligate seeding species also exist.Eucalypts also vary in their flammability (Gill & Zylstra, 2005) whereby different species within wet and dry sclerophyll forests may influence how a major fire impacts the forest (Gill, 1997).
Here we ask: (a) How variable is the post-fire TTR across southeast Australian eucalypt forests?; (b) What spatial and temporal factors drive this variation?;(c) Can we use these factors to construct a predictive ecological TTR forecast of the forests impacted by the 2019-2020 Australian megafires.To do this, we use the 20-year MODIS remote sensing record to examine factors affecting TTR and its underlying components, burn severity and the post-fire recovery rate.

Study Region
We focused our analyses on Eucalyptus-dominated woody ecosystems of southeast Australia (hereafter called "forests," Figure 1a), specifically the vegetation structural classes denoted as "Eucalypt Open Forests", "Eucalypt Tall Open Forests", and "Eucalypt Woodlands" by the National Vegetation Information System 5.1 (https://www.environment.gov.au/land/native-vegetation/national-vegetation-information-system).The large majority of these mixed-species eucalypt dominated areas are within 200 km of the coast and encompass more than 237,000 km 2 of southeast Australia (Figure 1a).The Eucalypt (Tall/Low) Open Forests are wetter with a mean annual precipitation spanning 500-1,600 mm yr 1 , while the drier and more inland Eucalypt Woodlands receive 260-1,050 mm yr 1 .

Remote Sensing Data Processing
Several sensors could be used to track vegetation recovery from fire, but only the Moderate Resolution Imaging Spectroradiometer (MODIS) has the multidecadal duration and temporal frequency needed to construct gap-free monthly time series for even the most cloud prone locations in the study region.The finer spatial resolution of Landsat and Sentinel-2 enable more accurate fire mapping (Chuvieco et al., 2019;Massetti et al., 2019), but the size of Australian bushfires tends to be large enough to be well observed by MODIS's coarser pixel resolution.The duration of the Sentinel-2 is generally too short to track historical recovery, and the 16-day revisit schedule of Landsat would leave many data gaps, especially in regions with persistent cloud cover.Since our priority was to trace the temporal dimension of recovery, we opted to use the daily retrievals of 500 m MODIS to construct gapfree monthly composites.
We designed a multi-stage process to calculate TTR, examine how TTR and its underlying components are affected by landscape and meteorological drivers, and to use these drivers to forecast TTR (Figure 2).We tracked TTR with the MODIS LAI product (MOD15A2H; Myneni et al., 2015) because of its consistency after extensive evaluation of six other MODIS derived vegetation indices (Supporting Information S1 Results).We used the 500 m MODIS MCD64 burned-area product (Giglio et al., 2021) to identify individual bushfire pixel locations.In order to more closely relate fire to the dominant seasonal cycles of precipitation, we defined the "fire-year" as starting in March and ending in February of the following calendar year.Time series for burned pixel locations were constructed by aggregating daily MOD15A2H LAI retrievals to a monthly time step for the period covering 2001-2020.Daily retrievals that were not flagged as "ideal-quality" (Quality Assurance, bits 0-1) were discarded.All remote sensing data sets (Table S1 in Supporting Information S1) were initially processed using Google Earth Engine (Gorelick et al., 2017) and post-processed in R with the "stars" (Pebesma, 2020) R package.

Predictors of Time-To-Recover
We investigated a range of both static landscape variables and temporally dynamic variables that might affect TTR (Table 1).Specifically these variables broadly correspond to environmental modulators of growth rate, such as properties relating to soil drainage and fertility.Landscape variables included edaphic properties, topographic variables, mean climate and mean annual LAI, and twentieth century fire history.Topographic variables (elevation, slope, aspect) were calculated using the Hydrosheds hydrologically conditioned digital elevation model (Lehner et al., 2008).Edaphic variables were extracted from the modeled high resolution (0.09 km) Soil and Landscape Grid of Australia (Viscarra Rossel et al., 2015) and averaged across all soil depths.These variables    S1 in Supporting Information S1).From AWAP, we calculated the monthly sums of precipitation (P), mean monthly daily temperature maximum (T max ), minimum (T min ), and vapor pressure deficit (VPD) calculated from the 3:00 p.m. vapor pressure measurement and T max .Shortwave radiation and mean air temperature was used to calculate potential evapotranspiration (PET) using a variant of the Priestley-Taylor method (Davis et al., 2017).The monthly aggregate of AWAP data was then used to calculate moving window (3-12 months) estimates of pre-and post-fire meteorological anomalies corresponding to the location and timing of every burned pixel observation.Meteorological data for the first three months of 2021 was not available from the AWAP product, so we recalibrated monthly means from the ERA5-Land product (Muñoz-Sabater et al., 2021) to AWAP using generalized additive models (GAMs) (Wood, 2017).ERA5-Land precipitation and VPD were highly correlated with nominal difference after calibration (R 2 : 0.94, RMSE 83 mm yr 1 for 12-month precipitation; R 2 : 0.99 and RMSE: 0.072 kPa for 12-month VPD).For pre-fire LAI anomalies, we calculated the pre-fire 3-and 12-month moving window anomalies.Two metrics of remotely sensed burn severity were derived.As one proxy of burn severity, we calculated the anomaly of normalized burn ratio (NBR; Key & Benson, 2006) from the 500 m resolution surface reflectance of the near infrared and shortwave infrared bands from MCD43A4 Nadir Bidirectional Reflectance Distribution Function Adjusted Reflectance (MODIS-MCD43) (Table S1 in Supporting Information S1; Schaaf & Wang, 2015).The anomaly was calculated relative to the long-term mean of the full data set as opposed to immediate pre-fire NBR, as is done with dNBR.Our reason for doing so is that NBR was often already anomalously negative pre-fire, owing to antecedent drought conditions.As another proxy of burn severity, we used the MOD15A2H LAI product to calculate the fraction of LAI reduction relative to the long-term seasonal mean of LAI (i.e., 1 -LAI post-fire /LAI normal ).

Estimating Time-To-Recover
We used an ecological resilience framework (Ingrisch & Bahn, 2018) to define time-to-recover as the duration of time between the disturbance event and when the vegetation index reaches a baseline near its pre-disturbance state.To better isolate the signal of drivers on TTR, we only analyzed pixel locations in eucalypt dominated woody ecosystems that burned once between 2001 and 2014 and that did not burn again until at least the start of the Black Summer fires in September 2019.We did not extend the analysis on TTR beyond 2014 in order to allow a sufficient time window for effectively all burnt pixel locations to recover.Only pixel locations that burned once were included in the analysis, thereby removing a small proportion of burned pixel locations in the data set that burned twice (12.2%; or three times: 1.1%) prior to the 2019-2020 megafires.
We developed two approaches to provide a robust estimation of TTR with the criteria that the estimation of TTR must be robust to the variable timing of peak leaf area and not conflate sporadic and short-lived phenological pulses of understory vegetation with recovery.The first approach defines TTR (TTR MW ) as the minimum time required for a 12-month moving window of LAI to increase within 0.25 standard deviations of the mean annual LAI (using 2001-2019 as a reference period).By integrating over 12 months, the TTR MW approach reduces the contribution of seasonal or transient understory non-tree vegetation on post-fire vegetation dynamics.We specifically did not use the mean annual vegetation index (e.g., LAI) as the baseline state for estimating TTR MW because many locations did not, or may not reach the pre-fire baseline again, owing to changes in rainfall or decadal-scale rainfall variability (Wardell-Johnson et al., 2017).We examined a range of threshold values and found that setting the recovery baseline to less than 0.25 standard deviations beneath the mean annual value led to a high rate of transient false-recoveries (not shown).In comparison, the 0.25 standard deviation threshold reduced the rate of false-recoveries but was sufficiently low for more than 99% of the pixel locations to recover the forest canopy within 7 years, which is consistent with the range from both field observations (Fox et al., 1979;Gill, 1978;Volkova et al., 2019) and remote sensing (Heath et al., 2016).
Earth's Future 10.1029/2023EF003780 We used a second approach of estimating TTR to interrogate the causal nature of TTR by fitting logistic growth functions to the LAI time series for each burned pixel location.While the TTR MW approach can smooth over transient pulses in phenology, a limitation is that it cannot estimate rapid TTR ≤365 days, and it provides no insight on the post-fire LAI recovery rate, burn severity, or the magnitude of the vegetation index at time of recovery.By contrast, the parameters of the logistic continuous population growth function resolve these components, allowing us to probe deeper into what drives TTR.Each pixel location's post-fire leaf area trajectory was fit with: Here L t is the LAI at time t (days), r is the post-fire LAI recovery rate (LAI day-1), K is carrying capacity of LAI which is akin to the stabilized LAI at the time of recovery, and L 0 is the remaining LAI following the fire.The estimated quantity 1 L 0 /K is similar to our LAI based metric of burn severity.We computed the deseasonalized 3-month mean LAI to fit Equation 1 for each burned pixel location using nonlinear least squares (nls.multstart;Padfield & Matheson, 2020) in R v4.01).Parameters were constrained to plausible ranges, where both r and L 0 were specified to be positive, and K was specified to range within 0.25 standard deviations of the mean annual LAI.For less severe bushfires (e.g., not crown fires), the full deleterious effect on LAI sometimes lags for several months after the fire.Therefore we used the minimum LAI within 6 months following the date of the fire as the starting point for the time-series.The logistic function requires the time series to stabilize near a threshold of recovery, however other effects (e.g., drought) can cause unrelated declines of LAI which can prevent identification of the model parameters if the time series is much longer than needed to identify TTR.Therefore each specific pixel location's time series was truncated 6-months following the moving window estimate of the TTR.By setting L t to the recovery baseline (mean annual LAI minus 0.25 standard deviations of the interannual LAI), and solving for t, we used Equation 1 to calculate the logistic function time-to-recover (TTR LF ) as:

Attributing the Drivers of Time-To-Recover and Its Sub-Components
The logistic function parameters have clear interpretations and allowed us to examine the specific subcomponents of post-fire canopy recovery dynamics.We subsequently used GAMs with penalized thin plate regression spline smoothing functions ("mgcv" R package, Wood, 2017) to identify linkages between the LAI recovery rate, burn severity (i.e., 1 LAI post-fire :LAI normal ), and TTR with predictor variables (Table 1).Although mean annual LAI is a key component of TTR, we did not investigate it because its determinants within the Australian landscape have been documented elsewhere (Yang et al., 2018).Each GAM was specified with penalized univariate smoothing functions because the purpose of these models was to identify and show the individual marginal effect of each covariate on the response variables.Only the six largest contributors to model goodness of fit were retained for each model.Each model was fit with 100,000 randomly selected burned pixel observations and evaluated on the remaining burned pixel observations (approximately 24,000).
We also investigated the legacy effects of fire history  on burn severity and the annualized leaf recovery rate following a similar GAM fitting procedure; however, each model was only focused on identifying an effect from fire history covariates (i.e., time since last disturbance, number of wildfires and burns since 1950, history of burns in the 1990s).Therefore models only included a univariate smoothing function on a fire history covariate and a spatial covariate using a two-dimensional smoothing function on longitude and latitude.

Ecological Forecasting the Time-To-Recover of the 2019-2020 Mega Fires
We used a non-parametric gradient boosting machine learning approach to predict the TTR of forest recovery following the 2019-2020 megafires.Gradient boosting is a commonly used machine learning technique with high predictive performance (Grinsztajn et al., 2022) that can suggest an upper limit of empirical approaches to forecasting with limited information (e.g., only data derivable from until 1 year following fire).While gradient Earth's Future 10.1029/2023EF003780 RIFAI ET AL.
boosting does not provide statistical inference, it provides quantitative insight about the importance of predictor variables.Specifically variable importance scores can suggest which components (e.g., static, or dynamic) most structure TTR variability, while Shapley values can be used to calculate the individual additive contribution of each predictor variable to the predicted TTR of an observation.We specified TTR MW as the response variable in the gradient boosting model and the same set of predictor variables used in the GAMs (see above) to predict TTR.Specifying either a Gaussian or Gamma response distribution had effectively no impact on predictive performance, so we present the predictions using a Gaussian response in order to aid the interpretation of the Shapley values which are additive when calculated from a Gaussian response.The H 2 O machine learning framework (LeDell et al., 2021) was used to perform a grid search of 30 gradient boosting models with 149,204 observations of the data, and evaluated the model fit on a withheld 10,000 subset of the observations.The best fit model was selected by evaluation of cross-validation goodness of fit metrics (mean absolute error and root mean square error) on four partitions of the data.We report out-of-sample model performance on the withheld 10,000 observation subset of the data.

Overview
The eucalypt forests of study area have been subject to both high rates of annual burning and high interannual variability of burn area.More than 39.6% of the region's eucalypt forests (Figure 1a) burned at least once between the 2001-2019 fire years.Burning was highly concentrated toward the summer months with approximately 95% of burns having occurred between September through February.The Black Summer megafires of 2019-2020 produced more than double the annual burnt eucalypt forest area relative to any year since the beginning of the MODIS record in 2001, burning a total of 34% of coastal southeast Australia's eucalyptus forests.Apart from major fire years (2002,2006,2019), there were also years where fire was relatively scarce in eucalypt forests (e.g., 2004,2010,2011,2015).Both the spread and median of TTR for LAI were highly variable across fire-years (Figure 3c).Annual burn area was correlated with annual precipitation anomalies (ρ = 0.53; Figure 3a) and fire-season VPD (ρ = 0.52; Figure 3b).High VPD also correlated with greater LAI loss from fire (ρ = 0.31; Figure 3b).Post-fire LAI anomalies became progressively more negative as the 2017-2019 increased in severity (Figure 3c).

How Variable Is Time-To-Recover in Time and Space?
The range of TTR was dynamic across fire years.The interannual median of TTR MW and TTR LF varied by a factor of 3.75 and 5.4, with coefficients of variation 0.54 and 0.47, respectively (Figure 3c).The percentage of burned pixel locations that had recovered by one year post-fire was also highly variable across fire years (0.8% in 2006 to 66.3% in 2014; Figures S1 and S2 in Supporting Information S1).
Intrinsic spatial structure of TTR was evident within individual fires (Figure S3 in Supporting Information S1), but TTR was unrelated to the size of an individual burn (ρ = 0.11), and TTR for pixel locations within the same contiguous bushfire could span a difference of more than 5 years (Figures 1b-1d).Figure 1 shows the LAI recovery trajectories from three notable bushfires, the Moogem State Forest fire in 2002, the Canberra/ACT fires in 2003 (Bryant, 2008), and the catastrophic Black Saturday fires in Kilmore East that occurred in 2009 (Cruz et al., 2012).Within each of these fires there were pixel locations that recovered quickly (≤365 days), and others that required several years (Figures 1b-1d).On average it can be seen that the LAI anomaly trajectory to the baseline was hastened by positive precipitation anomalies, whereas recovery stagnated when negative precipitation anomalies occurred after the fire.

What Drives Differences in the Time to Recover?
Overall TTR MW was most correlated with the LAI recovery rate (ρ = 0.65; r in Equations 1 and 2), less correlated with burn severity (ρ = 0.31), and uncorrelated with the post-fire leaf area carrying capacity K (ρ = 0.06).The duration of TTR was relatively unaffected by the approach of estimation with TTR MW and TTR LF being highly correlated (ρ = 0.92), when poor logistic function model fits were (R 2 < 0.2, 15% of pixellocation fits).Notably, both TTR MW and TTR LF were longest when some leaf area (∼25%) remained following a burn (Figures 4a and 4b).This can be explained by the decaying response of the LAI recovery rate with decreasing burn severity (Figure 4c).Specifically, the LAI recovery rate was highest at the highest levels of burn severity and Earth's Future 10.1029/2023EF003780 RIFAI ET AL.
decreased to a minimum at intermediate to low levels of burn severity (Figure 4c).In Figure 4d we illustrate how TTR would vary if a forest experienced different levels of burn severity.Here, a forest with no remaining leaf area (burn severity, 1 L 0 /K = 1) would recover faster than a less burned forest that retained 25% of its leaf area (1 L 0 /K = 0.75) because of the steep decline in growth rate at intermediate burn severity values (Figure 4c).
Next we used GAMs to help characterize the drivers of the annualized post-fire leaf recovery growth rate (Figure 5), in addition to burn severity (Figure S4 in Supporting Information S1) and TTR (Figure S5 in Supporting Information S1).
The annualized LAI recovery rate r parameter had only moderate explainability (Figure 5b; R 2 :0.26,MAE: 0.73 (LAI yr 1 )).Regional mean annual climate PET, the moisture index P:PET, and post-fire precipitation had the largest effects on the LAI recovery rate r.Burn severity had the next largest effect where the LAI recovery rate was greatest at the highest levels of burn severity.The LAI recovery rate was also higher in forests that burned earlier in the fire season, or that had more negative pre-fire 3-month LAI anomalies.We also found no substantial effect of twentieth century fire frequency on the post-fire leaf growth rate (not shown).The effect of recent burns in the 1990s appear to contribute to an additional 0.2 (±0.01)LAI yr 1 of the LAI recovery rate, but compared to other drivers, this effect was small (Figure S6a in Supporting Information S1).No appreciable effect of twentieth century fire frequency on burn severity was found.The large majority of 2001-2014 burned pixel locations were not known to have experienced fire in more than 60 years (Figure S6b in Supporting Information S1).

How Long Will It Take Forest Canopies to Recover From the Black Summer Fires?
The gradient boosting model to predict TTR yielded a high overall predictive ability (R 2 : 0.68, mean absolute error: 177 days).The mesic forests in northern New South Wales are predicted to be the fastest to recover, most of which in less than 2 years.Forests with the longest predicted recovery times (≥5 years) are mostly higher in elevation, such as the forests in the Blue Mountains surrounding the Sydney metropolitan region and the Australian Alps southwest of Canberra.In Figure 6 we show a subset of the most important variables and Shapley additive contributions to predicted TTR for the forecast model.The rank order of these variables is approximate when variable importance score differences are small and can vary depending upon minor details of model fitting such as data stratification for training and evaluation, or the optimization of hyperparameters.Therefore we urge some caution in the interpretation of the relevance of these variables for actually structuring TTR.The NBR anomaly and elevation were the two most important predictors of TTR.Seven of the 10 most important variables were dynamic in nature, such as the pre-and post-fire 12-month accumulation of precipitation, the pre-fire 3- month mean VPD, and the NBR anomaly immediately following the fire.Notably, the pre-fire status of the canopy (pre-fire 3-month and 12-month LAI anomalies) was again indicated as being of high importance.The remaining three variables were static, corresponding to elevation, mean annual LAI, and sand content.Despite the high burn severity of these fires, we predict the recovery from the 2019-2020 fires to be only marginally longer than previous fire years.A raster of the predictions is presented in Zenodo repository (https://zenodo.org/records/10981885) with https://zenodo.org/doi/10.5281/zenodo.10981884.

Quantifying Time-To-Recover From Remotely Sensed LAI
We estimated TTR in this study with a 12-month moving-window approach that is simple, yet robust.To further interrogate the drivers of TTR, we also statistically estimated TTR from the logistic growth function, which while more complicated to estimate, it provides crucial information on the LAI recovery rate.With the logistic growth function (Equation 1), TTR is the outcome of the LAI recovery rate, the burn severity, and the baseline LAI (Equation 2).Quantifying the LAI recovery rate allowed us to show that it is a larger determinant of TTR than burn severity or the baseline LAI.Both approaches agreed on the duration of time required for a burned pixel location to reach the recovery threshold.The definition of the recovery threshold can be challenging, as an overly low threshold risks conflating recovery with transient flushes of understory vegetation, whereas an overly high threshold risks artificially extending TTR due to unrelated changes in hydroclimate.The pre-fire LAI may have already been reduced by drought given that extreme dryness is a common pre-condition for major wildfires.Even defining the threshold of recovery on long-term means, or an optimal LAI in equilibrium with local hydroclimate (Yang et al., 2018) is complicated because changes in climate and atmospheric CO 2 have also affected the longterm mean of leaf area in Australian forests (Donohue et al., 2009;Ukkola et al., 2016;Rifai et al., 2022).Ultimately, this definition of recovery is limited to being specific to a singular forest attribute.For example, while LAI might reasonably be expected to reach its long-term mean within 10 years, it may take many more decades to recover pre-fire levels of biomass or diversity -if ever.
We concede that quantifying TTR of canopy LAI will not address all aspects of forest management.Specifically, no remotely sensed metric can reasonably indicate the multiple dimensions of forest recovery from fire that may be most relevant to management (e.g., biomass, stand age, species diversity).Remaining questions regarding what vegetation functional groups are recovering (e.g., trees, shrubs, grass) may be addressed by remote sensing techniques, but will require extensive field validation.The recent satellite constellations of higher resolution satellites (e.g.Sentinel-2; Landsat 8 and 9) when coupled with field observations could facilitate a more expansive analysis of functional changes during forest recovery.The future is promising because research may be able to integrate additional different forms of remote sensing data (e.g., radar, LIDAR, solar induced fluorescence) with field-based measurements to map burn area (Chuvieco et al., 2019) and to characterize temporal changes in the recovery of biomass and canopy photosynthesis.Moreover, again we emphasize that recovery is multidimensional in nature and future research could advance the field by deriving field-validated remote sensing measures that can jointly characterize the multifaceted aspects of forest function, structure, and diversity.

Why Does Intermediate Burn Severity Result in the Longest Time to Recover?
Here we showed that the LAI recovery rate is generally highest with very high burn severity (Figure 4a).The balance between the low growth rate of intermediately burned forests (Figures 4a and 4b) and the burden of regrowing 50%-75% of the mean annual leaf area leads to longer TTR (Figure 4d).There are several possible reasons why the LAI recovery rate might increase following low levels of burn severity."Cool burns" of lower severity have a limited impact on the canopy, while providing "competitive release" by burning understory vegetation, incinerating part of the competitor seedbank and liberating soil nutrients from burnt vegetation (Carter & Darwin Foster, 2004;Tomkins et al., 1991).Soil phosphorus liberated from fire may also be particularly important for accelerating post-fire eucalypt seedling growth (Tng et al., 2014).Reduced competition may allow canopy trees increased access to newly accessible soil nutrients and water resources.We lack field measurements to verify this, however, this line of reasoning is in line with the prescribed burning forest management literature (Certini, 2005;Chambers & Attiwill, 1994).
It is unclear why greater burn severity should produce the fastest LAI recovery rates.As with "cool burns," there is a transient increase in nutrient availability, but soil organic matter can also be incinerated and a higher fraction of the standing stock of nitrogen and phosphorus in forest litter may be volatilized by the fire (Carter & Darwin Foster, 2004).Post-fire nutrient availability may also be complicated by successional species.For example, nitrogen fixing pioneer species may rapidly colonize burn scars, eventually facilitating succession by species previously killed by fire (Adams & Attiwill, 1984).
The large majority of the forest eucalypt species with ranges encompassing the study area have evolved fire regenerative strategies, such as basal resprouting or resprouting from epicormic buds along the stem (Clarke et al., 2015).We postulate that because most eucalypts also develop storage organs for nonstructural carbohydrates, trees recovering from fire are more limited by water availability (possibly to transport nonstructural carbohydrates) than by the nonstructural carbohydrates themselves.However, eucalypt species exposed to increased fire frequency may be unable to accumulate and store sufficient quantities of nonstructural carbohydrates ("resprouter exhaustion"; Nolan, Collins, et al., 2021) such that carbon costs of fire become greater than excess nonstructural carbohydrates sequestered in storage organs between fire intervals.Obligate seeders may also be imperiled because they may be burned too frequently to allow replenishment of the seed bank (Bowman et al., 2016;Fairman et al., 2015;Keith, 1996).The interactive effects of drought on post-fire recovery may also lead to excessive xylem cavitation (Nolan, Gauthey, et al., 2021), further impairing photosynthesis and diminishing the excess photosynthate supply available to sequester nonstructural carbohydrates.This presents a particular problem for narrow-ranged endemic species, where macro-scale disturbance regime shifts have already impacted historically less-burned refugia such as the catastrophic fires that threatened the last remaining grove of Wollemi pines and engulfed some temperate rainforests during the 2019-2020 megafires (Kelly et al., 2020;Mackenzie et al., 2021).It may be useful for conservation efforts to identify species that lack fire regenerative traits (Gallagher, Butt, et al., 2021), have limited ranges, or are distributed in regions known for exposure to current high fire frequency (Gallagher, Allen, et al., 2021).
The GAM analysis suggested that upland vegetation on ridges (higher height above nearest drainage) that was more sun exposed experienced greater burn severity (Figure S4 in Supporting Information S1).This is largely consistent with earlier remote sensing studies (Heath et al., 2016) where lower fuel moisture was attributed to drive greater fire intensity (Bradstock et al., 2010).Drought driven hydraulic canopy dieback as a result of the 2017-2020 drought (Nolan, Collins, et al., 2021) may also have produced an unusually high fuel load.However this speculation about fuel loads is controversial due to a lack of field survey data to confirm it (Bowman et al., 2021).Nevertheless, we found the VPD anomaly in the months leading up to fire to be a strong predictor of burn severity (Jain et al., 2022) not only because it dries out understory fuels making them more flammable, but also because of the strong effect it has on accelerating leaf senescence in Australia (Jolly et al., 2005) thereby potentially augmenting fuel loads.In recent decades, warming from climate change has caused increases in VPD (Rifai et al., 2019), particularly in Australia (Rifai et al., 2022), and correspondingly increased the Forest Fire Danger Index (Clarke et al., 2013), leading to an expectation of increases in fire frequency and severity (Clarke & Evans, 2019).Again this is increasingly of conservation concern for many species with life-cycles longer than intervals between fire, but most especially obligate seeder species that may not have enough time to re-seed between increasingly brief fire-intervals (Fairman et al., 2015;Gallagher, Allen, et al., 2021).

Functional Relevance of Forest Canopy Recovery Time
The recovery of forest canopy leaf area is directly important for re-establishing post-fire ecosystem functions such as evapotranspiration and photosynthesis, however there are also notable ecological interactions with biodiversity implications that are contingent on the recovery of forest canopies.For instance, understory forest microclimates Earth's Future 10.1029/2023EF003780 which support herbs and forbs are directly mediated by overstory leaf area.In the water limited regions (∼68% of the study area), herbaceous understory species (e.g., species in families such as Lamiaceae or Goodeniaceae, or ferns) and midstory vegetation (e.g., semi-woody shrubs) are unlikely to develop and persist through drought and heatwaves without microclimates created by overstory shade.However we caution that recovery to a pre-fire LAI baseline is far shorter than the required time to recover pre-fire forest structural elements (e.g., biomass).Many questions remain regarding how forests allocate carbon following fire, which may have implications for recovery timescales of biomass and biodiversity.For example, there are few reports of how long it takes forests to recover pre-fire levels of fruit and seed production (Denham & Whelan, 2000), an absence of which may also affect fauna which depend on these food resources.Further, the sequence and abundance of fauna returning to post-fire ecosystems will be mediated by the cover from predation afforded by leaf canopy (Gallagher et al., 2022;Nimmo et al., 2021).

Forecasting Time-To-Recover With Limited Information
The forecast of eucalypt forest TTR from the 2019-2020 megafires (Figure 6) are, on average, longer than previous years, although considerably shorter than would be expected considering both the extensive spatial scale and high burn severity as indicated by reduced LAI (Figure 3b; but see Collins, Bradstock, et al., 2021).
Recent analysis shows a high degree of spectral recovery over more than half of the burned areas in the first year following the fires (Gibson & Hislop, 2022).The TTR forecast has been notably reduced due to exceptionally high rainfall in 2020 (20% higher than normal) brought about by prevailing La Niña conditions.However the large effects of post-fire meteorological conditions likely produces an upper limit on the predictability of TTR, therefore some degree of caution will be required when forecasts are produced for future fire years.For example, these predictions would be invalidated if these locations burned again within the 5-year forecast window.While it is unlikely that an appreciable quantity of coarse woody eucalypt fuel could accumulate in this period, the anomalously high rainfall in 2020 extending through 2022 could facilitate a press and pulse cycle for competing herbaceous vegetation as happened in the years following the wet 2010/2011 La Niña event (Harris et al., 2018).Nevertheless, we hope that our approach to producing testable predictions of ecosystem recovery from fire can be used to facilitate data-driven planning of prescribed burns, the identification of slow-recovering ecosystems, and tracking of species ranges which are of conservation concern.

Figure 1 .
Figure 1.Geographic extent of the woody eucalypt ecosystems of coastal southeast Australia.(a) The number of times burned between 2001 and 2019, and insets showing Landsat 5/7 false color images of Moogem fire (2002), ACT Canberra fires (2003), and Kilmore East fire (2009).Regions in dark gray are classified as woody eucalypt ecosystems that did not burn between 2001 and 2020.(b)-(d) The median LAI anomaly (black line) and 95% quantile distribution (gray ribbon) are shown.Example pixel locations trajectories corresponding to the 1% (fast) and 99% (slow) time-to-recover quantiles are shown in purple and orange.The fast and slow trajectories terminate when they have met our definition of recovery (see methods).Background vertical bars show the moving 3-month precipitation anomaly over each focal region.

Figure 2 .
Figure 2.An overall schematic of workflows is presented for calculating the time-to-recover (TTR; Part 1) and forecasting TTR from the Black Summer fires (Part 2).Gray squares represent source data or outputs, and rounded blue squares represent data processing and modeling steps.

Figure 3 .
Figure 3. (a) Time series of burned eucalypt forest area in the southeast Australian study region (Figure 1a).The zonal rolling 12-month mean of the percent precipitation anomaly is represented as vertical color tiles.(b) Boxplot distributions of pixellevel minimum LAI anomalies following fire.The redline shows the median post-fire LAI anomaly across all fire years.The zonal rolling 3-month mean of the percent vapor pressure deficit anomaly is represented as vertical color tiles.(c) Boxplots of the time-to-recover (Logistic growth function estimate) by fire year.The gray region starting in 2015 shows the period where the observation window (truncated in late 2019) was too short to estimate time-to-recover.

Figure 4 .
Figure 4. (a) A contour density plot is shown for 184,065 observations of the 12-month moving-window estimate of time-to-recover (TTR MW ) across the range of remotely sensed burn severity as indicated by the MODIS LAI product.A generalized additive model (GAM) trend line fit with a penalized thin plate regression spline is shown in black overlaying a contour observation density.(b) A contour density plot is shown with a GAM trend line for the logistic function estimated time-to-recover across the range of logistic function estimated burn severity.(c) The LAI recovery rate (r is plotted against burn severity for 184,065 individual pixel locations.A nonlinear GAM trend line fit with a penalized thin plate regression spline is indicated in black.(d) A diagram illustrating LAI recovery trajectories with the same carrying capacity (K), but different LAI recovery rates (r) and remaining leaf area following fire (L 0 ).The LAI recovery rate (r) was simulated from its relationship with burn severity (1 L 0 /K) spanning 1%-95% remaining leaf area using the GAM fit shown in panel (c).

Figure 5 .
Figure 5.The partial effects from a generalized additive model are shown with the six largest univariate contributors to the annualized LAI recovery rate (365 × r)."Month of Fire" was fit as a parametric term, and all other covariates shown were fit with penalized thin plate regression spline smoothing terms.Rug plots along the xaxis indicate the distribution of covariate data.Note the y-axis range varies by subpanel.

Figure 6 .
Figure 6.A map of the predicted time-to-recover for southeast Australia forested regions affected by the 2019-2020 bushfires.The top left inset shows frequency distributions of time-to-recover from the 2019-2020 prediction, in comparison with the four fire years with the highest amount of burning between 2001 and 2014.The middle left inset shows a summary distribution of Shapley additive contributions to predicted TTR for the top eight variables with the largest variable importance scores from the gradient boosting predictive model.The contribution of these variables is colored by the normalized value of each variable that has been rescaled to range between zero and one.
S.W.R., M.G.D.K., A.J.P., L.A.C, and P. M., acknowledge support from the Australian Research Council Discovery Grant (DP190101823).S.W.R., M.G.D.K., and A.J.P. acknowledge support from the ARC Centre of Excellence for Climate Extremes (CE170100023).M.G.D.K. was also supported by the NSW Research Attraction and Acceleration Program.Open access publishing facilitated by University of New South Wales, as part of the Wiley -University of New South Wales agreement via the Council of Australian University Librarians.

Table 1
Jones et al., 2009)erive TTR and Variables Used to Predict Time-To-Recover Change, Energy, the Environment and Water, 2010) and Victoria (Victorian Government Data Directory, Department ofEnvironment, Land, Water & Planning, n.d.).We calculated "mean annual" climate (MA) variables from a 30-year climatology of meteorological data (see below) encompassing the period of 1982-2010 as was the recent standard(World Meteorological Organization, 2017).Mean annual LAI (LAI MA ) was calculated using the 2000/02-2019/09 time-series of the MOD15A2H data product, and specifically excluded the time periods following the Black Summer fires.Temporally dynamic variables included meteorological anomalies, remotely sensed proxies of burn severity, and pre-fire LAI anomalies.Meteorological data were obtained from the Australian Bureau of Meteorology's Australian Water Availability Project (AWAP;Jones et al., 2009)0.05°griddedclimate product (Table Note.Detailed data set references are listed in TableS1in Supporting Information S1.Climate