Exploring the Inﬂuence of Forest Tenure and Protection Status on Post-Fire Recovery in Southeast Australia

: Research Highlights : We used Landsat time series data to investigate the role forest tenure and protection status play in the recovery of a forest after a ﬁre. Background and Objectives : Changing ﬁre regimes put forests in southeast Australia under increasing pressure. Our investigation aimed to explore the impact of different forest management structures on a forest’s resilience to ﬁre by looking at the post-ﬁre recovery duration. Materials and Methods : The analysis included a total of 60.6 Mha of land containing 25.4 Mha of forest in southeast Australia. Multispectral time series data from Landsat satellites and a local reference dataset were used to model attributes of disturbance and recovery over a period of 33 years. Results: Protected public forest spectrally recovered 0.4 years faster than protected private forest. No other signiﬁcant effects in relation to different tenure and protection status were found. Climatic and topographic variables were found to have much greater inﬂuence on post-ﬁre spectral recovery. Conclusions: Protected area status in public forests resulted in slightly faster recovery, compared with the private protected forest estate. However, factors outside the control of land managers and policy makers, i.e., climatic and topographic variables, appear to have a much greater impact on post-ﬁre recovery.


Introduction
In Australian forests, disturbances caused by wildfire have occurred regularly for millions of years. However, although many local species have adapted to regular fire through mechanisms such as epicormic resprouting and fire-resistant canopy-stored seeds [1,2], it has been suggested that changing fire regimes with higher burn frequencies and severities might pose a threat to some of Australia's unique forest ecosystems [3]. Recent studies have indicated that current forest management strategies in combination with climate change might lead to irreversible changes in some forested landscapes [4][5][6]. The fire season of summer 2019/20 can be regarded as an indication of changing fire regimes, as it was unprecedented in scale [7] and affected highly vulnerable and diverse forest ecosystems that are not typically exposed to regular fire disturbance [8]. Despite the scale and severity of the 2019/20 fire events, a relatively fast recovery has been documented, which has been attributed to above-average rainfall in the months after the fires [9]. However, other studies that examined native Australian forests have indicated that post-fire recovery usually takes several years [10,11].
Forest management policies in Australia vary with tenure and protection status. The management of protected areas under the National Reserve System (NRS) requires adherence to strict guidelines, usually not permitting timber extraction, and fire mitigation measures such as prescribed burning only under strict conditions [12]. Most protected areas can be found on public land; however, some are managed privately. On the other hand, most of the Australian forest is not under NRS protection, which permits private owners and public land stakeholders to have more invasive management strategies, including timber extraction and fire prevention measures. To ensure Australian forests retain their ecological key functions in the future it is necessary to evaluate how different tenure and protection status might impact forest resilience to fire. Of particular interest is recovery duration, as some forested ecosystems rely on a certain period of recovery before re-developing the ability to reproduce [13]. This study aims to explore the relationship between different forest management regimes and the length of post-fire recovery in southeast Australia.
Forest recovery is a complex process that includes a broad range of ecological and structural variables [10,14]. However, collecting detailed terrestrial information on these variables is cost intensive and often not feasible across large areas. Satellite-based remote sensing, on the other hand, offers wall-to-wall coverage at regular temporal intervals. For example, the Landsat archive provides multispectral satellite imagery reaching back to the 1970s, and has proven particularly valuable for forest change detection, as its spatial and temporal resolution is well suited to map forest-related phenomena [15][16][17]. Although multispectral remote sensing is not capable of capturing the full complexity of structural forest change, many studies have demonstrated that the spectral signals of disturbance and recovery are sufficiently correlated with the structural changes to address a broad range of monitoring questions [14,18]. The Normalised Burn Ratio (NBR) spectral index, in particular, has been found to be capable of capturing forest change [19,20], as it exploits the sensor's near infrared and short-wave infrared bands, which are more sensitive to structural complexity within the forest canopy [21].
To leverage data from the Landsat archive, a range of tools have been developed. A frequently used algorithm for processing Landsat time series data is LandTrendr [20,22], which fits linear trend segments to the spectral time series of each pixel, filtering noise and deriving generalised trends [21]. For this study, LandTrendr segmentation based on the NBR index was used to extract spectral metrics of disturbance and recovery for the study area in southeast Australia. The disturbance information obtained from LandTrendr was used in conjunction with further predictor variables to model higher-level disturbance attributes with the Random Forests algorithm, including fire severity and causal agent. To describe forest recovery, we adopted the Years to Recovery (Y2R) approach introduced by Pickell et al. [23], defined as the years it takes a pixel to return to 80 percent of its predisturbance spectral value (Y2R80) [14]. Subsequently, the Y2R80 metric was used to assess whether forest tenure and protection status had a significant correlation with recovery length. The key objective was to investigate whether land management and protection status influence a forest's recovery from wildfire, and how these factors compare with the impact of other environmental variables (e.g., climate and topography) in post-fir recovery.

Materials and Methods
To investigate how land tenure and protection status relate to post-fire recovery, a number of processing steps were undertaken. Firstly, disturbance and recovery maps were created for the study area using Landsat time series data from 1988 to 2021 and the segmentation algorithm LandTrendr [21,24]. Based on these maps, a random forest (RF) classifier was used to predict disturbance severity on the pixel level [25]. Disturbance patches were identified and attributed a disturbance agent using a second RF model [17,26], retaining only disturbances caused by fire. The remaining disturbance events were then used to examine how forest recovery behaviour varies across areas of different protection status and tenure in southeast Australia. The methods are discussed in detail below.

Study Area
The study area covered 25.4 Mha of forested land in southeast Australia. Its boundaries were defined as the extent of all 11 biogeographic regions (Interim Biogeographic Regionalisation for Australia, or IBRA) that can be found in the state of Victoria, some of which stretch across the border into the neighbouring states New South Wales and South Australia, as shown in Figure 1. Each IBRA region features distinct climatic, geological and ecological properties [27]. gionalisation for Australia, or IBRA) that can be found in the state of Victoria, some of which stretch across the border into the neighbouring states New South Wales and South Australia, as shown in Figure 1. Each IBRA region features distinct climatic, geological and ecological properties [27].
Sclerophyll forests account for a large percentage of all woody vegetation in the study area [28]. Many eucalypt species found in the study area have adapted to regular fire disturbances through mechanisms such as epicormic resprouting [1,10].

Restricting the Area of Interest
We restricted our analysis to forested areas within the study area, defined in Australia as areas with greater than 20 percent canopy cover with the potential to reach more than 2 m of height, over an area of more than 0.2 ha [29]. To create the corresponding forest mask, we used the National Forest and Sparse Woody Vegetation Data product of the Department of the Environment and Energy. The dataset provides a total of 23 nearannual layers of forest extent for Australia between 1988 and 2018, based on a classification of Landsat imagery [29]. Using these maps, a composite forest mask was produced to restrict the analysis to areas that were continuously forested between 1988 and 2021, whilst allowing for periods of disturbance and regrowth. Sparse woody vegetation, defined as areas with a canopy cover between 5 and 20 percent, was considered non-forest for the purpose of our analysis.
To investigate the correlation of forest tenure and protection status with recovery duration, we used tenure information from the Department of Agriculture and Water Resources Australia for the year 2013 [30] and data on protection status from the Commonwealth of Australia for the year 2016 [31]. Sclerophyll forests account for a large percentage of all woody vegetation in the study area [28]. Many eucalypt species found in the study area have adapted to regular fire disturbances through mechanisms such as epicormic resprouting [1,10].

Restricting the Area of Interest
We restricted our analysis to forested areas within the study area, defined in Australia as areas with greater than 20 percent canopy cover with the potential to reach more than 2 m of height, over an area of more than 0.2 ha [29]. To create the corresponding forest mask, we used the National Forest and Sparse Woody Vegetation Data product of the Department of the Environment and Energy. The dataset provides a total of 23 near-annual layers of forest extent for Australia between 1988 and 2018, based on a classification of Landsat imagery [29]. Using these maps, a composite forest mask was produced to restrict the analysis to areas that were continuously forested between 1988 and 2021, whilst allowing for periods of disturbance and regrowth. Sparse woody vegetation, defined as areas with a canopy cover between 5 and 20 percent, was considered non-forest for the purpose of our analysis.
To investigate the correlation of forest tenure and protection status with recovery duration, we used tenure information from the Department of Agriculture and Water Resources Australia for the year 2013 [30] and data on protection status from the Commonwealth of Australia for the year 2016 [31].

Forest Change Maps
We used Landsat L1T Surface Reflectance imagery [32] from the sensors ETM, ETM+ and OLI to create annual image stacks for the years 1988 to 2021. Only imagery from between February 1 and May 1 was considered, as this time window allowed most disturbance events happening during the Australian fire season to be detected in the same year as they occurred. For each year, a medoid composite was created [33][34][35]. Clouds and cloud shadow were filtered out using the FMASK algorithm [36]. All processing was carried out on the Google Earth Engine platform [37]. The medoid composites were then used to create forest change maps with LandTrendr [21,38], an algorithm widely acknowledged for its suitability for spectral segmentation of Landsat time series [15,16]. The LandTrendr algorithm uses the values of a given spectral index to fit trend segments to the time series of each pixel. We chose the Normalised Burned Ratio (NBR) index, as it has been shown to be well suited to the monitoring of fire disturbance [19,39] and is sensitive to structural forest change [14].
We chose parameters for the LandTrendr algorithm that were similar to the approach of Kennedy et al. [40]. Areas that were affected by a secondary disturbance were excluded from further analysis, given the objective of this study was to analyse uninterrupted recovery from fire.
In addition to the change maps derived from the LandTrendr analysis of the primary index NBR, the algorithm allows for the creation of change maps from additional spectral indices. The spectral values of the chosen indices are fitted to the trend segments obtained from the segmentation of the primary index [41]. For our analysis, we exported additional change map products created using the Tasseled Cap indices Brightness (TCB), Greenness (TCG) and Wetness (TCW). These indices have been found to respond to disturbance and recovery in different ways [20,22,42], and thus were expected to broaden the spectrum of forest changes that we could distinguish, such as different disturbance severities and causal disturbance agents.
The final outputs from LandTrendr included maps of the pre-disturbance value as well as the magnitude and relative magnitude of disturbance. Furthermore, a map of the disturbance onset year was generated. To exclude very small disturbance events, we applied a minimum mapping unit to the output, restricting the size of a change area to a minimum of 11 pixels (1 ha). Lastly, pixels that experienced an NBR change of less than 50 were also filtered out.

Predicting Disturbance Severity at a Pixel Level
The year of detection (YOD) of a pixel's primary disturbance event was directly obtained from LandTrendr; however, the disturbance severity was assigned through a RF classifier at the pixel level. RF is an ensemble classification technique that uses a large number of individual decision trees, each trained with a random subset of the reference data [43]. For our model, we set the number of trees to 200. Widely used in remote sensing applications, RF has been found to be a reliable classification tool even when trained on non-Gaussian distributed data [44]. In our RF model, we used a range of predictor variables that included spectral, climatic and topographic variables.
To train the RF classifier we relied on a reference dataset. The Victorian Forest Monitoring Program (VFMP) features 786 plots, with locations systematically stratified according to biogeographic region and public land tenure, namely multiple-use forests and nature conservation reserves [45]. Based on the VFMP framework, Soto-Berelov et al. [46] created an extensive reference dataset for forest disturbance in Victoria. In a 4 km 2 area surrounding each of the 786 VFMP plots, ten Landsat pixels were randomly chosen and attributed by human interpreters to determine disturbance events and associated severity.
To develop the RF model, we used disturbance severity information from the VFMPbased reference dataset. It contains 7823 reference plots with associated disturbance history, including severity as indicated by a number between 0 and 10 [46]. For the analysis, the values were reclassified into 3 categories: low (1-3), medium (4-7) and high severity (8)(9)(10). To ensure that the reference dataset and the LandTrendr-derived raster datasets were describing the same disturbance event for each reference pixel, only pixels where the onset years coincided (+/−1 year tolerance) were considered for training purposes. Undisturbed pixels or pixels that experienced a secondary disturbance were excluded, leaving a total of 613 disturbed training pixels. Predictor variables were extracted from LandTrendr disturbance segments and included NBR, TCB, TCG and TCW information on the pre-disturbance value (start value), as well as the spectral magnitude and relative spectral magnitude. Furthermore, climatic and topographic variables were included (mean temperature, annual rainfall, elevation, slope and aspect).
A total of 75 percent of the reference pixels were used to train the RF model, while 25 percent were held back for the assessment of model performance using the accuracy score metric of sklearn in Python [47]. Furthermore, we evaluated the relative importance of predictor variables using permutation importance [47].

Predicting Disturbance Agent at a Patch Level
Since our study was focused on recovery after fire disturbances, we proceeded to remove other types of disturbances (e.g., logging and drought) by classifying disturbance agents at the patch level. We defined a disturbance patch as a collection of adjacent pixels that were disturbed in the same year and with the same severity [20,41] by intersecting the YOD layer obtained from LandTrendr with the disturbance severity layer obtained from the first RF model described above.
To infer the causal disturbance agent for individual patches, we assigned a range of predictor variables to each patch, namely the mean and standard deviation of all contained pixels, computed from the same raster datasets as used in the pixel-level classification described in Section 2.4. Additionally, the area size and the circumference of each patch were determined and added to the range of predictor variables.
We used the same reference dataset as described in Section 2.4 to train the RF classifier for the patch-level classification of disturbance agents [46]. The original description of a disturbance agent in the reference data comprised a total of 10 different types of disturbance; however, we simplified these classes for our analysis and regrouped the values into three categories that covered the two most relevant disturbance agents within the study area: fire, logging and other. Every patch that contained one or more reference pixels was considered a candidate for the training dataset.
To ensure the integrity of the training data, a range of consistency checks were performed. If prospective training patches contained two or more reference pixels, which suggested different disturbance agents or onset years, the patch was excluded from training purposes. Furthermore, for each training patch, the YOD value as indicated by the reference dataset was checked against the value suggested by LandTrendr. If the two differed by more than one year, we assumed that different disturbance events were depicted, and removed the patch from the training dataset. A total of 216 patches passed the consistency checks, of which 75 percent were used to develop the RF model for the prediction of the disturbance agent on the patch level. The resulting RF model was assessed using a held-back portion of the patches (25 percent), and the relative importance of predictor variables was investigated using permutation importance [47].

Assessing Forest Recovery
For the recovery analysis, we used the annual NBR values of the time series fitted to the segments provided by the LandTrendr algorithm. The advantage of this approach is that, for each pixel, it yields a stabilised time series that is less affected by year-to-year noise resulting from sun angle, atmospheric effects, phenology or other factors [41].
To assess the impact of land tenure and protection status on the recovery process of disturbed forest, we defined recovery duration as the time it takes a pixel to return to 80 percent of its pre-disturbance NBR value [23]. The Years to Recovery 80 percent metric (Y2R80) has been previously tested by White et al. [14] for describing forest regrowth after clearfell logging, and found to be a good indicator of the state of actual structural recovery. Other studies also found the Y2R80 metric to be a suitable descriptor of recovery after fire disturbance [23,48].
To be eligible for the final analysis of recovery duration, a range of conditions were placed on Y2R80 pixels. The mean NBR value of the six years prior to the onset of a disturbance event was used to establish a pre-disturbance reference value. Hence, pixels disturbed before the year 1994 were excluded from the analysis due to a lack of pre-disturbance years. Furthermore, no Y2R80 pixels burnt after the year 2011 were considered, leaving pixels disturbed in more recent years a maximum of 10 years to recover. Consequently, we imposed this maximum as a rule on all disturbed pixels for consistency. Thus, every pixel that needed longer than 10 years to recover to 80 percent of its pre-disturbance NBR value was considered non-recovered.

Examining the Effects of Tenure and Protection on Recovery
As a means of evaluating the relative importance of forest tenure and protection status for recovery, we used a random subset (n = 100,000) of all Y2R80 pixels in a RF regression model to predict recovery length. A range of other predictor variables were also included in the model, including disturbance severity, elevation, annual rainfall and mean temperature.
Additionally, violin plots were used to assess post-fire recovery across various forest tenure and protection status categories. To examine the findings further, a series of one-way Analyses of Variance (ANOVA) tests were conducted.

Forest Change Maps
The change mapping using the LandTrendr algorithm resulted in a total area of 7.3 Mha of disturbed forest, including pixels that experienced a secondary disturbance. Figure 2 shows an example of resulting disturbance information (namely, year of disturbance) for a selected part of the study area. Investigation of the temporal agreement (year of disturbance) between the VFMP-based reference data and the LandTrendr output yielded an overall agreement of 63.5 percent. This comparison allowed for a tolerance of +/− one year, which accounted for disturbance events happening outside the annual Landsat observation window from 1 February to 1 May. Year of detection (YOD) of disturbance mapped for a selected part of the study area. The YOD information was derived from a Landsat time series using the LandTrendr segmentation algorithm. Non-forested areas and pixels that experienced a secondary disturbance were masked out.

Classification of Disturbance Severity
For pixel-level modelling of disturbance severity, the most important predictor variables were the start value and relative change magnitude of NBR, the relative magnitude of TCW, the year of disturbance and temperature ( Table 1). The accuracy score of the RFdriven severity classification was 0.65. Table 1. Feature importance (permutation importance) of predictor variables in a random forest model predicting pixel-level disturbance severity. Year of detection (YOD) of disturbance mapped for a selected part of the study area. The YOD information was derived from a Landsat time series using the LandTrendr segmentation algorithm. Non-forested areas and pixels that experienced a secondary disturbance were masked out.

Classification of Disturbance Severity
For pixel-level modelling of disturbance severity, the most important predictor variables were the start value and relative change magnitude of NBR, the relative magnitude  Table 1). The accuracy score of the RF-driven severity classification was 0.65.

Classification of Causal Agent of Disturbance
The RF classifier used to predict the disturbance causal agent was based on a total of 38 predictor variables. The most important predictor variables were the patch-based mean and standard deviation of the TCW relative magnitude of disturbance, as well as the mean TCW pre-disturbance value. Furthermore, the year of disturbance, mean slope and standard deviation of elevation also demonstrated a high predictive importance. The model had an accuracy score of 0.82. Figure 3 shows examples of classified disturbance patches with the associated disturbance agent.

Impact of Tenure and Protection on Recovery
In total, 14.5 percent of all Y2R80 pixels did not spectrally recover within 10 years, and thus were considered unrecovered and excluded from the subsequent analysis ( Table 2). The remaining eligible Y2R80 pixels (Figure 4) in our analysis together covered an area of 942,629 ha. The RF classifier used to predict the disturbance causal agent was based on a total of 38 predictor variables. The most important predictor variables were the patch-based mean and standard deviation of the TCW relative magnitude of disturbance, as well as the mean TCW pre-disturbance value. Furthermore, the year of disturbance, mean slope and standard deviation of elevation also demonstrated a high predictive importance. The model had an accuracy score of 0.82. Figure 3 shows examples of classified disturbance patches with the associated disturbance agent. Figure 3. Random forest-based predictions of causal disturbance agent for disturbance patches in a selected part of the study area. Non-forested areas and pixels that experienced a secondary disturbance were masked out.

Impact of Tenure and Protection on Recovery
In total, 14.5 percent of all Y2R80 pixels did not spectrally recover within 10 years, and thus were considered unrecovered and excluded from the subsequent analysis ( Table  2). The remaining eligible Y2R80 pixels (Figure 4) in our analysis together covered an area of 942,629 ha.     Central to the analysis was to investigate the correlation of forest tenure and protection status with the recovery duration. Firstly, a random subsample (n = 100,000) of the recovery dataset was used to predict Y2R80 using RF regression [47]. Table 3 shows the relative importance of all examined predictor variables. The RF model ranked protection status as the least important predictor variable, while the forest tenure had a slightly higher importance. Both factors, however, were of minor importance when compared with environmental predictor variables and fire severity. Subsequently, we created violin plots of the Y2R80 distribution across forest tenure and protection status categories and conducted ANOVA tests using randomly subsampled groups (n = 1000 in each). Figure 5 demonstrates that the mean values of all categories are similar. ANOVA tests and subsequent pairwise comparisons yielded no significant differences between most groups (p > 0.05), except for those between the protected public and protected private group (p < 0.005) which on average took 4.5 and 4.9 years to recover to pre-disturbance spectral levels, respectively. Subsequently, we created violin plots of the Y2R80 distribution across forest tenure and protection status categories and conducted ANOVA tests using randomly subsampled groups (n = 1000 in each). Figure 5 demonstrates that the mean values of all categories are similar. ANOVA tests and subsequent pairwise comparisons yielded no significant differences between most groups (p > 0.05), except for those between the protected public and protected private group (p < 0.005) which on average took 4.5 and 4.9 years to recover to pre-disturbance spectral levels, respectively. However, the described differences between tenure and protection categories were eclipsed by the impact of other factors that we identified as strong predictor variables of Y2R80, such as fire severity, as demonstrated by the violin plots of Figure 6. According to ANOVA tests, the fire severity categories low, medium and high all differed significantly. Each of the severity categories showed a unique distribution of Y2R80 values. Most pixels disturbed by low-severity fires recovered within less than 2 years, while pixels disturbed by medium-severity fires displayed a wider range of recovery responses. High-severity fires caused a bi-modal distribution of Y2R80 values, with a statistical mode of 7 years. However, the described differences between tenure and protection categories were eclipsed by the impact of other factors that we identified as strong predictor variables of Y2R80, such as fire severity, as demonstrated by the violin plots of Figure 6. According to ANOVA tests, the fire severity categories low, medium and high all differed significantly. Each of the severity categories showed a unique distribution of Y2R80 values. Most pixels disturbed by low-severity fires recovered within less than 2 years, while pixels disturbed by medium-severity fires displayed a wider range of recovery responses. High-severity fires caused a bi-modal distribution of Y2R80 values, with a statistical mode of 7 years.

Discussion
In this study, we examined the impact of forest tenure, protection status and environmental variables on spectral recovery duration (Y2R80) of forest after fire. The results showed that forest tenure and protection status were mostly unrelated to Y2R80. Almost no significant differences were found between the recovery of private and public tenure forests, and between protected and unprotected forests. Our analysis furthermore confirmed that, compared with other factors such as climatic and topographic variables, forest tenure and protection status were very weak predictors of recovery length after fire events. Considering these two factors are the only variables investigated in this study that land managers and policy makers can influence, the results of this study suggest that forest resilience to fire might not depend on forest-management-related factors as much as other studies have suggested [4]. On the other hand, mean annual temperatures in southeast Australia have exceeded the average temperature of the reference period (1961 to 1990) in every year between 2000 and 2022 [49]. The results of this study suggest that this development might have far greater implications for post-fire forest recovery than forest tenure and protection status. However, our analysis found one significant difference, showing that protected forests under public tenure spectrally recovered 0.4 years faster than privately held forest under protection status. This finding was backed by the low number of unrecovered pixels (Y2R80 > 10 years) in public protected forests, where only 6.2 percent of the disturbed forest had not spectrally recovered within 10 years of the fire, compared with 21.3 percent in private protected forests. Although reviewing different protection policies across Australia was out of the scope of this study, this result might be of interest to land managers and should be further explored in future studies.
One limitation of our analysis was its focus on pixels that experienced only one disturbance during the whole study period (1988 to 2021). While this ensured a spectrally stable pre-disturbance period for each pixel and prevented recovery periods from being interrupted by a secondary disturbance, it also might have concealed some effects of Australian fire management policies. Prescribed burns are a common feature of fire mitigation strategies, aiming at reducing the fuel load, but can vary in their utilisation depending on

Discussion
In this study, we examined the impact of forest tenure, protection status and environmental variables on spectral recovery duration (Y2R80) of forest after fire. The results showed that forest tenure and protection status were mostly unrelated to Y2R80. Almost no significant differences were found between the recovery of private and public tenure forests, and between protected and unprotected forests. Our analysis furthermore confirmed that, compared with other factors such as climatic and topographic variables, forest tenure and protection status were very weak predictors of recovery length after fire events. Considering these two factors are the only variables investigated in this study that land managers and policy makers can influence, the results of this study suggest that forest resilience to fire might not depend on forest-management-related factors as much as other studies have suggested [4]. On the other hand, mean annual temperatures in southeast Australia have exceeded the average temperature of the reference period (1961 to 1990) in every year between 2000 and 2022 [49]. The results of this study suggest that this development might have far greater implications for post-fire forest recovery than forest tenure and protection status. However, our analysis found one significant difference, showing that protected forests under public tenure spectrally recovered 0.4 years faster than privately held forest under protection status. This finding was backed by the low number of unrecovered pixels (Y2R80 > 10 years) in public protected forests, where only 6.2 percent of the disturbed forest had not spectrally recovered within 10 years of the fire, compared with 21.3 percent in private protected forests. Although reviewing different protection policies across Australia was out of the scope of this study, this result might be of interest to land managers and should be further explored in future studies.
One limitation of our analysis was its focus on pixels that experienced only one disturbance during the whole study period (1988 to 2021). While this ensured a spectrally stable pre-disturbance period for each pixel and prevented recovery periods from being interrupted by a secondary disturbance, it also might have concealed some effects of Australian fire management policies. Prescribed burns are a common feature of fire mitigation strategies, aiming at reducing the fuel load, but can vary in their utilisation depending on forest tenure and protection status [50]. However, as prescribed burns can be registered by LandTrendr as a secondary disturbance signal, it can be assumed that some pixels experiencing such provisions were masked out from our analysis. Thus, we expected our recovery analysis to be partly insensitive to possible effects of these fire management strategies on forest and its recovery.
The use of Landsat imagery in conjunction with an algorithm for time series segmentation such as LandTrendr is widely used for forest disturbance mapping [15,35,41]. In our study, we combined LandTrendr output with a local reference dataset of forest disturbance to model attributes of forest disturbance and recovery. The temporal agreement between disturbance maps derived from LandTrendr and the VFMP-based reference dataset was 63.5 percent. It can be assumed that the NBR magnitude threshold of LandTrendr affected this level of agreement, which allowed only disturbances with an NBR loss of greater than 50 to be detected. This threshold value might not have been sensitive enough to detect all disturbance events registered using the multiple sources of evidence training data [46]. However, it can be assumed that the effect of this discrepancy on the RF models was limited since only reference pixels in agreement with the LandTrendr output were used for training purposes.
Our first RF model predicted disturbance severity at a pixel level and achieved an OOB score of 0.61. In comparison, the authors of the VFMP-based training dataset reported an overall accuracy of 68 percent for a similar RF-based modelling approach [46]. One source of confusion named by the authors was the high rate of disagreement between interpreters for pixels of low disturbance severity, which might have affected the results of our study as well. Our second RF classifier for modelling the causal disturbance agent on patch level achieved an OOB score of 0.90. Our findings on the importance of predictor variables align with earlier studies conducted in the area of southeast Australia, confirming spectral magnitude and relative magnitude as valuable predictors of disturbance agent [20]. One notable difference, however, was the relatively low importance that disturbance patch shape metrics had in our model.
In the second part of our study, we examined recovery after fire disturbance for forested pixels. The Y2R80 metric relies on a stable pre-disturbance period and an uninterrupted recovery period, so we decided to examine only the recovery of disturbance events that occurred between 1994 and 2011. This left at least six years of pre-disturbance history for each pixel from the beginning of our Landsat time series in 1988, and a maximum of ten years for pixels disturbed in 2011 to recover before the end of the time series. Consequently, we decided to label all pixels that did not recover within 10 years as unrecovered, to avoid bias relating to events occurring late in the time series [23,51] and make the comparison of Y2R80 pixels more coherent. This resulted in the discarding of 14.5 percent of otherwise eligible pixels from the analysis. Lowering the recovery threshold (e.g., using a 60 percent instead of an 80 percent threshold) would have yielded a higher percentage of pixels being eligible for the final analysis. However, White et al. [14] observed that the Y2R60 metric overestimated recovery and showed a narrower distribution of values, indicating a decrease in sensitivity.
Lastly, it must be emphasised that spectral recovery, as measured in this study, does not adequately reflect all elements of structural recovery [10,18,52]. Vertical forest structure is a key factor in biomass monitoring [53,54] and fire modelling [55], as well as habitat mapping [10,56]. Despite the sensitivity of the NBR index to some elements of vegetation structure [14,39], passive remote sensing data generally have limits in capturing vertical forest structure in detail. Thus, it is increasingly used in conjunction with active remote sensing technologies such as spaceborne lidar [18,[57][58][59] or synthetic aperture radar [60], as well as airborne lidar [55,[61][62][63][64]. Future investigations into the role of forest tenure and protection policies in post-fire recovery would benefit from taking advantage of such approaches.

Conclusions
We used Landsat imagery from 1988 to 2021 and the LandTrendr algorithm to create disturbance and recovery maps for southeast Australia across 11 bioregions. Applying random forest models and a local reference dataset, disturbance severity and causal agent were modelled. We then conducted an analysis of recovery (Y2R80) for a total of 942,629 ha of forested pixels affected by fire events between 1994 and 2011, excluding areas that experienced a secondary disturbance. The findings suggest that the variables of interest, forest tenure and protection status, were mostly unrelated to recovery duration. The only significant difference was found for protected forests, where public forest recovered 0.4 years faster than its private equivalent. However, other factors such as climatic and topographic variables were much more strongly associated with post-fire recovery duration.