Identifying Post-Fire Recovery Trajectories and Driving Factors Using Landsat Time Series in Fire-Prone Mediterranean Pine Forests

Abstract

Wildfires constitute the most important natural disturbance of Mediterranean forests, driving vegetation dynamics. Although Mediterranean species have developed ecological post-fire recovery strategies, the impacts of climate change and changes in fire regimes may endanger their resilience capacity. This study aims at assessing post-fire recovery dynamics at different stages in two large fires that occurred in Mediterranean pine forests (Spain) using temporal segmentation of the Landsat time series (1994–2018). Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) was used to derive trajectory metrics from Tasseled Cap Wetness (TCW), sensitive to canopy moisture and structure, and Tasseled Cap Angle (TCA), related to vegetation cover gradients. Different groups of post-fire trajectories were identified through K-means clustering of the Recovery Ratios (RR) from fitted trajectories: continuous recovery, continuous recovery with slope changes, continuous recovery stabilized and non-continuous recovery. The influence of pre-fire conditions, fire severity, topographic variables and post-fire climate on recovery rates for each recovery category at successional stages was analyzed through Geographically Weighted Regression (GWR). The modeling results indicated that pine forest recovery rates were highly sensitive to post-fire climate in the mid and long-term and to fire severity in the short-term, but less influenced by topographic conditions (adjusted R-squared ranged from 0.58 to 0.88 and from 0.54 to 0.93 for TCA and TCW, respectively). Recovery estimation was assessed through orthophotos, showing a high accuracy (Dice Coefficient ranged from 0.81 to 0.97 and from 0.74 to 0.96 for TCA and TCW, respectively). This study provides new insights into the post-fire recovery dynamics at successional stages and driving factors. The proposed method could be an approach to model the recovery for the Mediterranean areas and help managers in determining which areas may not be able to recover naturally.

1. Introduction

Wildfires constitute one of the most widespread and important natural disturbances of forest ecosystems, playing a paramount role in the dynamics of the terrestrial system [1]. Forest fires impact at a wide range of scales causing ecological, economic and human health impacts [2,3]. Specifically in Europe, the Mediterranean region registers the highest number of fires and burned areas [4], with around 85% of the total burnt area [5].

Notwithstanding, Mediterranean ecosystems are adapted to fire recurrence as it constitutes the most important natural disturbance, driving vegetation dynamics [6]. Mediterranean species have developed post-fire ecological strategies including resprouting capacity, seed bank persistence and increased dispersal capacity [7,8]. Nevertheless, land use changes and the impacts of climate change may affect the dynamics of post-fire ecological succession in the immediate future [3,9]. Although large fire (i.e., ≥500 ha) occurrence for the European Mediterranean region does not show a strong increasing trend in the recent decades [5], climate change projections indicate an increase in the frequency and intensity of megafires, as a result of more extended and severe seasonal droughts [10], which will impact ecosystems’ species composition and functioning [3]. Forest ecosystems must adapt not only to changes in average climatic variables, but also to a wide variability with higher risk of extreme climatic events, such as prolonged droughts. Thus, forest management in European Mediterranean countries is challenging due to the vulnerability of natural regrowth capability of these ecosystems [11,12].

Time series of satellite data have long been used for retrospectively generating information on forest disturbance and recovery dynamics [13]. The opening of the Landsat archive in 2008, now available geometrically and radiometrically corrected, provided new opportunities for improved understanding of the mechanisms of forest changes [14,15]. Several studies have addressed the spatial and temporal analysis of post-fire vegetation dynamics through different forest ecosystems: Mediterranean [16,17], boreal [18,19], Siberian [20,21], temperate [22], tropical [23], savannah [24] or across different ecozones at the regional or national scale [25,26,27,28].

The use of the Landsat time series (LTS) for change detection has increased substantially in recent years as new methodological approaches have emerged [29]. Early approaches characterized post-fire recovery dynamics by applying linear regression functions to spectral trajectories obtained from the Landsat time series [30,31]. More recently, several change detection algorithms have been developed and widely used in analyzing forest changes, such as Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) [32] and Vegetation Change Tracker (VCT) [33], to provide change information on an annual time-scale [29]. Others include, Breaks For Additive Seasonal and Trend (BFAST) [34] and Continuous Monitoring of Forest Disturbance Algorithm (CMFDA) [35], which use a high frequency of time series. The trajectory-based segmentation algorithm LandTrendr enables the characterization of distinct subtrends within a simplified representation of the spectral trajectory, which provides the essential information needed to identify abrupt disturbances in forests (e.g., fire and harvest), as well as slowly evolving processes (e.g., regrowth and defoliation). The utility of LandTrendr has been demonstrated in different regions for assessing disturbance and recovery dynamics [25,28,36].

Several spectral measures can be derived from LTS and used as inputs for segmentation algorithms such as spectral indices or Tasseled Cap Transformations (TCT). Some spectral indices focused more on the red and near-infrared bands, making them sensitive to canopy greenness and photosynthetic activity, such as the Normalized Difference Vegetation Index (NDVI) employed for characterizing post-fire recovery [16,18,37], whereas other indices using the SWIR bands are more sensitive to vegetation moisture and forest structure [38], such as the Normalized Burn Ratio (NBR) [39], commonly used for recovery assessment [19,28,40]. TCT are created via linear transformations using defined coefficients [41] and have been widely used for studying forest changes [36,42,43,44]. TCT components correspond to the physical characteristics of vegetation: Brightness (TCB) is related to the pixel albedo of the land surface and values are typically high after a stand replacing disturbance; Greenness (TCG) is a contrast between the visible and near-infrared bands, being sensitive to green vegetation [19], and Wetness (TCW) is a contrast of the visible and near-infrared with the SWIR bands, making it sensitive to canopy moisture and structure [45]. Several metrics can be derived from TCT such as the TC Angle (TCA), which is related to the vegetation cover within the TCB-TCG spectral plane [46]. Considering that spectral indices are sensitive to different vegetation conditions, the use of TC components and derived metrics enables the characterization of different forest conditions [19,38,44,47].

Although the dynamic of post-fire vegetation recovery has been studied through different forest ecosystems, few studies have investigated the recovery driving factors in Mediterranean ecosystems [16,30,37,48] and fewer have focused on characterizing successional recovery stages [49,50], which are key to understanding forest changes for sustainable forest management. This study assesses the post-fire recovery dynamics at different stages in fire-prone Mediterranean pine forests. The specific objectives of this study were: (1) to identify the different post-fire recovery trajectories using temporal segmentation of LTS; (2) to analyze the recovery patterns for each trajectory group through stages and (3) to appraise the environmental and contextual drivers of the recovery process.

2. Materials and Methods

2.1. Study Area

This research was based on two large fires that occurred in the summer of 1994 (Figure 1): the Yeste Fire (7 August), which burned 11,685 ha of wooded area, and the Requena Fire (5 July), which burned 16,373 ha of wooded area. For this study, we selected sections that had neither burned in subsequent fires nor been reforested after the main disturbance of 1994 in order to ensure the analysis of natural recovery only.

Both study areas are located in the Southeast of the Iberian Peninsula, in the Mediterranean biogeographic region, which is characterized by mean annual precipitations of 600–700 mm with soil hydrological deficit in summer and mean annual temperatures around 15 °C. These areas were dominated by anthropogenic coniferous forests, mainly composed of species of the genus Pinus along with certain deciduous species of the genus Quercus, and, sclerophyll species such as Rosmarinus, Thymus or Juniperus species in the understory [51]. Due to the differences in post-fire ecological strategies to recover, we selected those patches dominated by Pinus halepensis and Pinus pinaster according to the Second National Forest Inventory of Spain (SNFI) [52]. Both species are obligate seeders since they have serotinous cones that enable the natural post-fire regeneration [7,53]. The post-fire recovery stages that can be identified in a Mediterranean pine forest 24 years after fire, range from the stand initiation (establishment phase including remnant pines, herbaceous and pine seeding processes) [8] to the stem exclusion (a young regrowth forest composed by shrubs and tree plantlets) [49,54,55]. Since competition between shrubs and trees starts immediately following fire [50,56], vegetation recovery in this study refers to both tree and shrub recovery.

2.2. Data

We downloaded the Landsat TM/ETM+/OLI images from the United States Geological Survey (USGS) Earth Explorer server [57] to build the time series covering the period 1990–2018, including 4 years pre-fire (Path/Row: 200/033, 199/032, 199/033). We selected images from Tier 1 Surface Reflectance products generated from the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) software [58] for TM and ETM + images and Landsat 8 Surface Reflectance Code (LaSRC) for the OLI dataset [59]. We prioritized scenes with less than 10% cloud-cover, within the summer period to minimize the effect of phenological changes. To delimit fire perimeters and to assess fire severity we used two Landsat 5 TM images acquired in 1994 for the Yeste Fire (22 July pre-fire, 23 August post-fire) and for the Requena Fire (29 June pre-fire, 16 August post-fire) (Figure 1). Fire perimeters were constructed by applying the USGS thresholds [39] to the differenced Normalized Burn Ratio (dNBR) [60].

Spatial reference to vegetation types was based upon the Forest Map of Spain, which was made between 1986 and 1997 through aerial photographs and field work [61]. Topographic variables were built from the LiDAR-based Digital Elevation Model (25-m spatial resolution) from the National Geographic Institute (IGN) of Spain [62]. As climatic information we used the Standardized Precipitation-Evapotranspiration Index (SPEI) [63], a multi-scalar drought index that calculates the effect of potential evapotranspiration (PET) on drought severity. Compared to other drought indices, the SPEI has the advantage of combining multi-scalar character with the capacity to include the effects of temperature variability [64]. Data are available for the entire time series (1990–2018) at 1-km resolution. We also downloaded the orthophotos at 0.5-m resolution from the Aerial Orthophotography National Plan, as reference data to assess recovery through time series (Years 2002, 2009, 2010, 2017 and 2018) [65].

2.3. Methods

A flowchart of the methodology is depicted in Figure 2. Firstly, annual composites were created for the available time series. Secondly, spectral indices and TCT was performed to delineate burned areas and derive trajectory metrics of forested pixels based on the LandTrendr segmentation algorithm [32]. Fitted trajectories were categorized according to the change magnitude and duration represented in the sequence of segments of each trajectory using a K-means algorithm. Finally, we analyzed the environmental and contextual drivers of the recovery process. The accuracy assessment of the recovery classes was carried out by visual assessment of vegetation cover in randomly distributed sample plots, in reference to high-resolution orthophotos.

2.3.1. Landsat Time Series

Summer time series was created for the 28-year time period using the closest cloud-free image to the mid of the summer season: Requena median Julian day 209 and Yeste median Julian day 218. The temporal window used to select the images spanned ± 38 days around the reference date (Figure 3) in order to ensure consistency through the time series.

We used the R package LandsatLinkr [66,67] to create annual cloud-free image stacks and calculate the TCT components using the coefficients defined for reflectance data [41]: Brightness (TCB), Greenness (TCG) and Wetness (TCW). Subsequently, the angular component of the TCT (TCA) was computed as follows [46]:

$$TCA={\tan }^{-1}\left(\frac{TCG}{TCB}\right)$$

(1)

2.3.2. Trajectory Segmentation and Clustering

To extract recovery trajectories, we applied the LandTrendr trajectory-based segmentation algorithm [32] to TCA and TCW time series. In this study we opted to use the TCA due to its relation with the percentage of vegetation cover in coniferous and mixed forests [46,47], and the TCW since it is sensitive to canopy moisture and structure [44,45]. LandTrendr goes through the time series and creates a fitted trajectory as a sequence of line segments for each pixel. Firstly, the vertices of each segment were established from an iterative regression process using Ordinary Least Squares (OLS) by estimating the years of change using the TCA and TCW time series. Then, the trajectories are iteratively simplified from a selection process using the angle criterion until a number of segments equal to or less than a user-defined threshold were obtained (segmentation process) [32]. We set this threshold to the maximum available (6) as we attempted to unravel multiple recovery trends. In a second step, the spectral values at candidate vertices were estimated (fitting process), generating a trajectory of connected segments for each pixel. The best model was chosen based on the p-value according to the F-statistics (p < 0.05). Further details on the segmentation process can be found in Kennedy et al. (2010) [32].

From the derived trajectories, we selected the change magnitude and duration of the segments, at the pixel level, as parameters for the classification. We calculated a recovery ratio (RR) for each segment (Equation (2)), which allows us to describe the recovery rates through time or alternatively for each successional stage defined by the number of segments (Figure 4).

$$\mathrm{RR}=\frac{{\mathrm{Magnitude}\mathrm{of}\mathrm{Change}}_x}{{\mathrm{Duration}\mathrm{of}\mathrm{Change}}_x}$$

(2)

We performed an unsupervised clustering method, since our first objective was to unravel the different vegetation recovery patterns. K-means clustering [68] is one of the most popular unsupervised machine learning algorithms and has been previously applied to summarize vegetation types and changes [69]. K-means allowed us to group the trajectories (defined by the RR of the segments) in k groups, minimizing the sum of the distances between each trajectory and the centroid of a given class. To define the optimal number of classes we employed the elbow method, which is based on the percentage of variance explained as a function of the number of clusters [70]. This process has been carried out with the scikit-learn library of Python [71]. The clusters obtained were overlapped to the LandTrendr outputs to characterize the categories.

2.3.3. Assessing Driving Factors of Vegetation Recovery

We aimed to explain the influence of fire severity, pre-fire conditions, topographic and climatic variables on recovery ratios through regression analysis. Six explanatory variables obtained from the Landsat imagery as well as auxiliary data sources (see Section 2.2) were derived to model post-fire recovery (

Table 1. List of variables used in the regression analysis.

Variable		Units	Description
Dependent	Recovery Ratiox (RR-TCAx RR-TCWx)	Z value	Represents the slope of the fitted trajectory at each segment
Explanatory
Pre-fire conditions	TCA90-93 or TCW90-93	Z value	TCA shows the percent vegetation cover and TCW the moisture and structure before the fire
Fire severity	dNBR	Values between −1 and 1	Represents the short-term post-fire effects on vegetation cover and structure. Severity thresholds proposed by the USGS [39]: Low: 0.1 ≤ dNBR < 0.27 Moderate-low: 0.27 ≤ dNBR < 0.44 Moderate-high: 0.44 ≤ dNBR < 0.66 High: ≥0.66
Topography	Elevation	Meters
	Slope	Percent
	Aspect	Values between 0 and 1	(TRASP) [72]. Values of 0 correspond to cooler, wetter north-northeastern aspects; values of 1 correspond to hotter, dryer south-southwestern aspects
Climatic Anomalies	Drought index	Z value	(SPEI) [63]. Positive values represent positive water balance and negative values indicate drought conditions

). For the variables that did not meet the assumptions of normality of the residuals and homogeneity of variance, we used log-transformed or rank-transformed data.

Previous studies have shown strong effects of fire severity on post-fire vegetation recovery [12,16,56]. Fire severity, defined as the degree of ecosystem change caused by a fire with respect to the pre-fire situation [73], was evaluated through the dNBR [60]:

$$dNBR=NB{R}_{pre-fire}-NB{R}_{post-fire}$$

(3)

$$NBR=\frac{\left(NIR-SWIR2\right)}{\left(SWIR2+NIR\right)}$$

(4)

Previous work has also addressed the important role of topography in explaining variations in forest establishment following fire [18,20,30]. Other studies have established the relevance of the post-fire climate [16,28,37] since it is related to water availability, and pre-fire vegetation conditions [20] due to its relationship with post-fire seed availability [12]. Regarding the drought index (SPEI), the 3 month-scale (cumulative from June, July and August) was selected since vegetation activity responds predominantly to short drought time-scales [74] and because maximum Pearson correlation coefficients between SPEI aggregated from summer season and TCA–TCW time series were recorded for both study areas.

Firstly, an exploratory regression analysis was carried out to diagnose the suitability of the selected variables. Due to the presence of spatial autocorrelation and heteroscedasticity in our data (the significance of Koenker statistics at the 95% confidence level), we executed a Geographically Weighted Regression (GWR) [75,76], a local regression model that considers spatial heterogeneity in data relationships. Model fitting was conducted using optimized adaptive bi-square kernel bandwidth (according to the corrected Akaike information criterion) [77].

2.3.4. Recovery Assessment

We evaluated the recovery through TCA and TCW trajectories in the short, mid and long-term as Key and Benson (2006) [39] proposed to assess burn severity: 2002 (8 years post-fire), 2009–2010 (15–16 years post-fire) and 2017–2018 (23–24 years post-fire), with 80% of the pre-fire value of TCA and TCW as the recovery threshold [19].

We carried out a stratified validation based on the recovery categories identified with a sample size of 500 plots randomly selected at each phase for each index (a total of 3000 reference plots). Vegetation recovery was evaluated by visual analysis of high resolution orthophotos. To facilitate the visual interpretation, we divided the 30-m pixel of Landsat images equally with a 6 m by 6 m grid (25 cells in one 30-m pixel) and overlaid the grids onto the orthophotos (Figure 5), similarly to Zhao et al. (2016) [78]. As an approximation to the pre-fire fractional cover due to the lack of pre-fire orthophotos, we established as reference the pre-fire fractional cover obtained from the SNFI (mean cover of 42.9% and 47.9% for Requena and Yeste, respectively). A disturbed pixel was considered to have recovered if the tree and shrub cover was at least 40% (i.e., 10/25; Figure 5a,b). Otherwise, that pixel was interpreted as having not recovered (Figure 5c).

Four accuracy metrics derived from confusion matrices were computed to validate our vegetation recovery classification (

Table 2. Confusion matrix example.

	Reference Data
Estimation	Recovered	Non-Recovered	Row Total
Recovered	P11	P12	P1+
Non-recovered	P21	P22	P2+
Col. total	P+1	P+2	N
OE	P21/P+1		(5)
CE	P12/P1+		(6)
OA	P11 + P22/N		(7)
DC	2P11/(P1+ + P+1)		(8)

), the omission error (OE; Equation (5)), the commission error (CE; Equation (6)), the overall accuracy (OA; Equation (7)) and the Dice coefficient (DC; Equation (8)) [79].

3. Results

3.1. Classification of Post-Fire Trajectories

Four different categories were identified for TCA (CR, CRSC, CRS and NCR) and five for TCW (CR, CR2, CRSC, CRSC2 and CRS). The main characteristics of the categories describing recovery are defined in

Table 3. Recovery categories definition.

Category	Acronym	Stages	Description
Continuous Recovery	CR CR2	1 2	Pixels show a continuous increase in the TCA and TCW values since the year of fire (CR) or since the first year post-fire (CR2).
Continuous Recovery with Slope Changes	CRSC CRSC2	4 3	Continuous recovery follows disturbance but slope changes occur through the time series. Changes occur at a different time for TCA and TCW (CRSC and CRSC2).
Continuous Recovery Stabilized	CRS	2	Continuous recovery, which slow down or stop 4–5 years after fire.
Non-continuous Recovery	NCR	3	Recovery process is interrupted in the mid-term followed by a second phase of continuous recovery (only found with TCA).

.

Recovery dynamics in TCA and TCW tended to be spatially clustered, indicating strong spatial effects (Figure 6). Moreover, there are differences in terms of the magnitude of change, as well as the year in which changes detected between TCA and TCW occurred, and both fires. For the Requena Fire we observed greater homogeneity, with CR as the main category, whereas for the Yeste Fire the CRSC category predominates.

The plots in Figure 7 show the mean fitted trajectories of recovery categories from the TCA and the TCW. In the pre-fire period, the TCA and TCW trajectories did not show significant changes, indicating relative stability in forest cover until the occurrence of fire. Nevertheless, the trend is negative in the case of TCA trajectories, which could be due to a loss of vegetation vigor before the fires occurred. The structure seems to have not changed in the pre-fire period according to stable values in TCW trajectories.

Although burned areas were generally characterized by an increase in spectral values after the fire events, the mean TCA and TCW trajectories showed differences in vegetation recovery between the two fires and across time series. TCA showed faster recovery, with high slopes in trajectories in the short-term, because it is associated with the percent of vegetation cover, whereas TCW showed slower recovery since it is related to the vegetation structure and moisture.

Regarding the TCA categories, the same trends were found in both fires and post-fire mean values tended to overtake pre-fire values in the mid-term. Even so, recovery rates were higher in Requena, as the post-fire mean values reached the pre-fire values before 2005, whereas in Yeste the pre-fire values were reached around 2010. CR pixels correspond to lower disturbance magnitude whereas in CRSC, CRS and NCR the magnitude of change is clearly higher. In the cases of CRSC and NCR the recovery rates were high even though the recovery dynamic was interrupted. CRS showed a particular trajectory as it stabilized in the mid-term and maintained the values across a subsequent phase of slight recovery.

According to the TCW categories, different trends were found between fires with the exception of CR and CRS pixels. In the case of the Yeste Fire the categories CR2 and CRSC continued to decrease until one year after the fire. In the case of Requena Fire, two different categories of CRSC were found since breakpoints in the recovery process occurred in different years. In spite of the trend toward pre-fire values none of the TCW trajectories reached pre-fire conditions after 24 years.

3.2. Assessing Drivers of Post-Fire Vegetation Recovery

We summarized the results from regression analysis for TCA and TCW in

Table 4. Modeling results for Tasseled Cap Angle (TCA) recovery categories at each stage.

Category	Variable	Stage 1		Stage 2		Stage 3		Stage 4
		Coefficient	Standard Error	Coefficient	Standard Error	Coefficient	Standard Error	Coefficient	Standard Error
CR	Intercept	0.662	0.145
	Pre-fire conditions	0.012	0.005
	Fire severity	0.061	0.005
	Elevation	−0.018	0.017
	Slope	0.003	0.005
	Aspect	−0,010	0.004
	Drought Index	0.279	0.492
		R2: 0.77; Adjusted R2: 0.76; AICc: 955.89
CRSC	Intercept	0.651	0.142	0.384	0.351	−0.847	1.421	−0.565	0.341
	Pre-fire conditions	0.075	0.043	0.055	0.068	0.058	0.277	0.081	0.066
	Fire severity	0.553	0.046	0.265	0.332	0.064	0.362	0.122	0.327
	Elevation	−0.261	0.069	−0.135	0.143	−0.162	0.586	−0.068	0.140
	Slope	−0.046	0.040	0.025	0.058	0.157	0.238	0.043	0.057
	Aspect	0.034	0.034	0.008	0.048	−0.008	0.197	−0.015	0.047
	Drought Index	0.758	0.154	1.143	0.058	2.525	0.130	0.788	0.047
		R2: 0.77; Adjusted R2: 0.76; AICc: 1656..55		R2: 0.81; Adjusted R2: 0.80; AICc: 4031.20		R2: 0.75; Adjusted R2: 0.74; AICc: 11055.78		R2: 0.73; Adjusted R2: 0.72; AICc: 4238.93
CRS	Intercept	−0.495	0.068	0.196	0.053
	Pre-fire conditions	0.093	0.025	0.024	0.009
	Fire severity	0.260	0.032	0.040	0.055
	Elevation	−0.132	0.064	−0.001	0.032
	Slope	−0.006	0.026	−0.004	0.012
	Aspect	−0.024	0.022	−0.002	0.008
	Drought Index	1.849	0.065	0.075	0.086
		R2: 0.88; Adjusted R2: 0.88; AICc: 1746.11		R2: 0.61; Adjusted R2: 0.58; AICc: 12831.31
NCR	Intercept	−0.052	0.101	0.344	0.349	−0.540	0.294
	Pre-fire conditions	0.079	0.032	0.033	0.064	0.048	0.054
	Fire severity	0.449	0.039	−0.826	0.353	0.582	0.300
	Elevation	−0.234	0.063	0.213	0.165	−0.201	0.141
	Slope	−0.018	0.036	0.051	0.070	−0.001	0.059
	Aspect	−0.041	0.029	0.021	0.051	−0.032	0.044
	Drought Index	1.613	0.110	1.065	0.060	0.845	0.231
		R2: 0.83; Adjusted R2: 0.82; AICc: 1527.51		R2: 0.78; Adjusted R2: 0.77; AICc: 7082.34		R2: 0.74; Adjusted R2: 0.72; AICc: 6085.56

and

Table 5. Modeling results for Tasseled Cap Wetness (TCW) recovery categories at each stage.

Category	Variable	Stage 1		Stage 2		Stage 3		Stage 4
		Coefficient	Standard Error	Coefficient	Standard Error	Coefficient	Standard Error	Coefficient	Standard Error
CR	Intercept	0.416	0.040
	Pre-fire conditions	−0.008	0.002
	Fire severity	0.023	0.002
	Elevation	−0.002	0.005
	Slope	−0.005	0.002
	Aspect	0.013	0.002
	Drought Index	0.121	0.446
		R2: 0.71; Adjusted R2: 0.69; AICc: 955.23
CR2	Intercept	0.025	0.101	−0.076	0.046
	Pre-fire conditions	−0.312	0.046	0.049	0.020
	Fire severity	0.106	0.046	0.006	0.019
	Elevation	−0.158	0.065	0.027	0.036
	Slope	0.072	0.036	0.004	0.015
	Aspect	−0.171	0.039	0.025	0.016
	Drought Index	0.734	0.156	−0.609	0.170
		R2: 0.80; Adjusted R2: 0.79; AICc: 7870.95		R2: 0.92; Adjusted R2: 0.92; AICc: 1227.11
CRSC	Intercept	−0.322	0.120	0.612	0.214	0.358	0.229	−0.119	0.167
	Pre-fire conditions	−0.116	0.051	−0.091	0.116	−0.008	0.126	0.024	0.092
	Fire severity	0.094	0.044	0.206	0.099	0.140	0.108	0.053	0.079
	Elevation	−0.040	0.072	−0.098	0.248	−0.005	0.270	−0.059	0.198
	Slope	0.041	0.038	0.004	0.091	0.021	0.099	0.014	0.073
	Aspect	−0.148	0.031	0.103	0.070	0.106	0.077	0.050	0.056
	Drought Index	0.122	0.126	0.437	0.088	0.737	0.058	0.640	0.052
		R2: 0.88; Adjusted R2: 0.87; AICc: 1959.93		R2: 0.67; Adjusted R2: 0.65; AICc: 15831.16		R2: 0.62; Adjusted R2: 0.60; AICc: 17032.23		R2: 0.56; Adjusted R2: 0.54; AICc: 13009.93
CRSC2	Intercept	0.765	0.065	0.324	0.198	0.099	0.047
	Pre-fire conditions	−0.061	0.041	0.085	0.113	−0.008	0.026
	Fire severity	0.045	0.039	0.139	0.106	0.018	0.025
	Elevation	−0.127	0.062	−0.166	0.247	−0.017	0.057
	Slope	0.014	0.030	0.030	0.087	−0.009	0.020
	Aspect	−0.046	0.027	0.195	0.070	0.025	0.016
	Drought Index	0.730	0.039	0.735	0.077	−0.150	0.075
		R2: 0.81; Adjusted R2: 0.80; AICc: 1064.20		R2: 0.69; Adjusted R2: 0.67; AICc: 11669.50		R2: 0.70; Adjusted R2: 0.69; AICc: 1451.47
CRS	Intercept	0.146	0.061	−0.030	0.074
	Pre-fire conditions	−0.047	0.028	−0.002	0.031
	Fire severity	0.014	0.028	0.045	0.031
	Elevation	−0.030	0.058	−0.005	0.093
	Slope	0.018	0.022	−0.007	0.026
	Aspect	−0.021	0.019	0.034	0.021
	Drought Index	0.244	0.095	0.012	0.121
		R2: 0.94; Adjusted R2: 0.93; AICc: 2685.86		R2: 0.78; Adjusted R2: 0.77; AICc: 1967.958

, respectively. All of the aforementioned variables for predicting the Recovery Ratio (RR-TCA_x and RR-TCW_x) at each stage were statistically significant at the 95% level. The successional stages correspond to the segments of the fitted trajectories from each category.

With regards to the relationship between the RR-TCA_x and the predictor variables, different responses were found between the categories but also common trends among recovery classes (Table 4). There was a positive influence of pre-fire conditions since a higher percent of vegetation cover prior to the fire will lead to a higher percent of vegetation cover also after the fire. In addition, climate was positively related to RR-TCA_x in all stages. Positive water balance resulted in a higher ratio of percent vegetation cover. The coefficients varied for the different stages since post-fire climate increases its explanatory power until stage 3 and drops in the long-term. In relation to severity, high recovery rates were related to high severities whereas low ratios were associated with low-burned pixels. Regarding topographic variables, coefficients varied according to the stage and category. Elevation showed a negative relation as a result of the limiting effects of temperature on vegetation growth, being more important in the short-term and decreasing through subsequent stages. The slope showed a weaker, negative influence as well as the aspect since cooler, wetter north-northeastern aspects (i.e., lower aspect values according to TRASP) were preferred.

The relationship between the RR-TCWx and the predictor variables showed different responses reflecting changes in forest structural complexity (Table 5). In this case, a negative relationship with pre-fire values was found because the greater complexity of the forest structure (i.e., higher TCW), the lower the recovery ratio, indicating slower recovery processes of remnant trees and seeding in contrast to the quicker recovery of shrubs. With regard to severity higher values also lead to higher recovery rates. Concerning topographic variables, the slightly positive relationship of the slope might be attributed to pine forest distribution preferably in the foothills. Identically to TCA, the elevation and aspect relationship tended to be negative although in the mid-term stages south-southwestern aspects lead to higher recovery rates. Post-fire climate conditions showed a positive or negative relationship depending on the stage. Enough available moisture post-fire is important for seed germination but abrupt changes in the climate conditions seem to reduce the recovery ratio in the long-term.

The relative importance of the explanatory variables for each trajectory category was assessed through t-statistics (Figure 8 and Figure 9). Post-fire climate in terms of drought had high predictive power in most of the stages and categories according to both TCA and TCW regression analysis. This power increased in the mid-term and long-term while all other variables decreased. Fire severity had the second largest power for explaining percent vegetation cover in the short-term, although its power diminished in the estimation of the recovery in terms of forest structure. Pre-fire conditions also had a higher importance in the short-term, since the seed bank and the seeding processes will depend on the pre-fire forest cover and structure. Topographic variables showed the lowest explanation power. Elevation had higher relative importance in the case of TCA compared to TCW, in which aspect was the most important topographic variable in all stages.

3.3. Recovery Estimation Assessment

According to the recovered and non-recovered definitions, all categories showed high accuracy, with OA values ranging from 0.7 to 0.94 (

Table 6. Errors metrics for recovery estimation according to TCA and TCW.

	Category	2002				2009–2010				2017–2018
		OA	DC	OE	CE	OA	DC	OE	CE	OA	DC	OE	CE
TCA	CR	0.86	0.91	0.07	0.11	0.90	0.94	0.04	0.07	0.95	0.97	0.00	0.05
	CRSC	0.70	0.84	0.14	0.19	0.82	0.90	0.09	0.11	0.87	0.93	0.07	0.07
	CRS	0.79	0.86	0.05	0.22	0.86	0.92	0.08	0.08	0.94	0.97	0.02	0.05
	NCR	0.74	0.81	0.16	0.20	0.88	0.93	0.10	0.04	0.92	0.96	0.03	0.05
TCW	CR	0.91	0.85	0.21	0.04	0.84	0.88	0.20	0.02	0.93	0.96	0.08	0.00
	CR2	0.82	0.74	0.36	0.13	0.86	0.85	0.21	0.08	0.89	0.93	0.13	0.00
	CRSC	0.82	0.78	0.33	0.06	0.87	0.89	0.15	0.01	0.82	0.91	0.16	0.01
	CRSC2	0.82	0.82	0.29	0.04	0.86	0.91	0.16	0.00	0.88	0.92	0.14	0.00
	CRS	0.86	0.87	0.17	0.10	0.77	0.93	0.12	0.01	0.89	0.93	0.10	0.04

OA—Overall Accuracy; DC—Dice coefficient; OE—omission error and CE—commission error.

). The OE and CE of the recovery categories varied, with the highest errors in the short-term (OE of 0.36 for TCW and CE of 0.22 for TCA). Post-fire recovery estimated by TCA shows more balanced errors, although with higher CE due to early successional recovery processes of herbs and shrubs. In contrast, recovery estimated by means of TCW shows higher OE since it is more sensitive to moisture and structure than early soil colonization of herbs. In both cases, more stable recovery classes (CR and CRS) showed higher accuracy than the more disrupted ones (CRSC). DC increased through the time series since forest cover and structure were more clearly defined.

4. Discussion

4.1. Post-Fire Recovery Trajectories from LTS

Characterizing post-fire recovery processes is challenging due to the variety of factors driving vegetation recovery, resulting in different recovery dynamics. Post-fire vegetation recovery estimated from spectral data is not a direct measure of actual forest regrowth. However, trends of forest recovery can be quantified by linking spectral change metrics with a reference dataset [50]. Here, we identified different recovery categories according to TCA and TCW recovery ratios at different stages through a 24-year-series in Mediterranean pine forests.

Time series analysis from Landsat data using the LandTrendr segmentation algorithm has been suitable to capture the different post-fire recovery trends. The trajectories extracted revealed continuous and non-continuous recovery processes, allowing us to identify slight changes in the slowly evolving recovery process. This was of great importance in determining slow but more stable recovery processes (CR and CR2) compared to other faster, but also interrupted, recovery processes (CRSC, CRSC2, CRS and NCR), indicating changes in greenery and forest structure throughout the recovering process. Other studies that also employed a Landsat trajectory-based approach with LandTrendr were able to identify different patterns of vegetative regrowth depending on the state, owner category and ecoregion in North America [25] or set recovery levels across sclerophyll forest in Australia [36]. Even though several studies have addressed the analysis of post-fire recovery trends, our findings highlight the importance of defining and grouping recovery patterns to facilitate the understanding of recovery processes.

TCA and TCW have been proved useful to characterize both vegetation cover (TCA) and forest structure (TCW). TCW trajectories were much more gradual, while TCA trajectories tended to saturate around 5 years post-disturbance. Frazier et al. (2015) [44] and Nguyen et al. (2018) [36] also concluded that the TCW contained more detail on the vegetation structure and regrowth in the regeneration processes since it is highly correlated to stand age and structural complexity in mature forest stands. The wetness values rise with an increasing amount of canopy [45], making it more suitable for analyzing mid and long-term recovery. Alternatively, initial increases in vegetation cover were well-characterized with shorter visible and near-infrared wavelengths as used in TCA [43,47,80]. In contrast to the canopy layer, recovery of the understory through both resprouting and seeding is quicker [8,54]. Thus, TCA was more suitable for tracking early stages, suggesting greater sensitivity toward shrub recovery and changes in vegetation condition rather than structure.

The Recovery Ratio (RR) varied across the stages and for each category (Figure 7). Generally, the recovery rates of TCA and TCW were greatest shortly after the fire (Stage 1 in CRSC, CRSC2, CRS and NCR), decreasing afterwards, due to the early post-fire colonization of annual herbs and shrubs [8,56]. In the following stages (2, 3 and 4), the RR according to both TCA and TCW was lower, associated with stem exclusion processes in a young regrowth forest characterized by intense competition among regenerated species [7,49,53]. Our results agree with other studies that also obtained higher recovery rates in the short-term according to spectral vegetation indices in Mediterranean forests [40] and NBR trajectories in pine, mixed conifer and conifer–oak forest [28].

Accordingly, the recovery time for the fitted mean trajectories also varied across categories. TCA categories tended to reach mean pre-fire values quicker, as it tends to saturate earlier (short-mid-term) due to the influence of herbaceous vegetation on TCA. Specifically, in the Requena Fire, the TCA trajectories tended to overtake pre-fire values in the long-term, possibly because these are fire-adapted forests in which fire creates favorable conditions for vegetation germination and regeneration [7]. Nevertheless, TCW trajectories did not reach the previous values 24 years after the fire, suggesting that burned areas did not recover the complexity of the pre-fire forest structure. This agrees with the interval of minimum 15 years to consider a Pinus halepensis forest recovered after a fire proposed by Eugenio et al. (2006) [54] since post-fire populations of seeder species does not overpass the reproductive juvenile phase up to 12–20 years after fire [81] and thus, the canopy seed bank is not completely fulfilled [12]. Some studies in Mediterranean ecosystems reported recovery times from remote sensing, which fit with our findings [82]. In this sense, Gouveia et al. (2010) [69] found recovery times of vegetation cover around 3–5 years according to NDVI (highly correlated with TCA), and Fernández-Manso et al. (2016) [17] estimated the time of vegetation cover with VRI between 7 and 20 years depending on fire severity level. However, the estimated recovery times with NBR, which is highly correlated with TCW, were generally longer compared to NDVI [40].

Some of the limitations for the modeling of post-fire vegetation recovery using optical data are related to saturation at high biomass levels. Previous studies found that the saturation of optical indices is reached after 20 years in Mediterranean environments [83]. Likewise, Schroeder et al. (2011) [80], Pickell et al. (2016) [19], Viana-Soto et al. (2017) [37] and Hislop et al. (2018) [38] found post-fire recovery of NDVI or TCA returning back to pre-fire levels rapidly (i.e., around 5–7 years post-fire). Structural information derived from airborne LiDAR data would enable a better characterization of the recovery trajectories and to improve our understanding post-fire vegetation recovery [84,85].

Accuracy Assessment of Post-Fire Recovery

Accuracy assessment of post-fire recovery is often avoided because historical reference datasets are scarce and field data is costly and time-consuming [14]. Here, we used a human interpretation approach to derive reference data of recovered and non-recovered areas, which has been widely employed to derive reference data for disturbance and recovery mapping [23,36,78,86].

All classes showed high accuracy with increasing OA and DC from the short-term to the long-term since forest cover and structure are more clearly defined 24 years post-fire. DC values were slightly lower in the long-term for TCW compared to TCA, which might indicate a recovery process that evolves to a secondary forest with higher shrub domain [12]. This could also be related to the uncertainty associated with signal sensitivity to changes in vegetation cover and biomass [87]. Moreover, the highest accuracy in the more stable categories (CR), in contrast to the more variable (CRSC), also highlights the difficulty in establishing the level of recovery for those more dynamic areas. Nguyen et al. (2018) [36] also pointed out the challenge of determining the post-disturbance recovery level, and identified recovery levels from NBR trajectories after fires according to whether pre-fire conditions were reached or not. Here, we did not distinguish among recovery levels as we were unable to accurately characterize the pre-fire forest structure due to the lack of high-resolution imagery, for which a more extensive reference dataset would also be needed.

The main source of error stemming from OE for both TCA and TCW in the recovered pixels. Our results showed that the TCW had higher accuracy for the non-recovered areas but omitted pixels that had already recovered. Zhao et al. (2016) [78] also found higher OE in the recovered class regarding post-fire and post-harvest forest recovery. Similarly, DeVries et al. (2015) [23] assessed post-harvest regrowth in tropical forests obtaining lower CE in the regrowth class but the highest OE. On the other hand, TCA showed higher CE, as it is more sensitive to fast detection of early recovering processes but also tends to saturate earlier as Schroeder et al. (2011) [80] reported.

4.2. Assessment of Post-Fire Recovery Drivers

Regression modeling of TCA and TCW recovery ratios for the recovery categories identified showed a varied influence of environmental factors, fire severity and pre-fire conditions. The results indicated post-fire climate as one of the most important factors for vegetation recovery in Mediterranean pine forests in Spain. Likewise, Meng et al. (2015) [16], Liu (2016) [88] and Viana-Soto et al. (2017) [37] found that climate conditions in the first post-fire seasons were critical for predicting short-term recovery. Tree regeneration after disturbances in Mediterranean ecosystems could be limited under post-fire drought events since droughts constrain seedling establishment and growth [7,53]. Bright et al. (2019) [28] also reported that post-fire climate explained substantial variation of the degree to which vegetation greenness recovered after a fire. Further analysis between trajectories and post-fire climate revealed that stages of recovery slowdown and even breakpoints coincided with negative SPEI values (i.e., dry or very dry periods). The year of the fires was followed by a slightly humid period, which supported the recovery. However, 5 years after the fires, a new drought event interrupted the recovery, as can be observed very clearly in the stabilization of the recovery in the CRS category from 1999–2000, not only in the TCA trajectories but also in TCW trajectories. The effect of this drought event was also noticeable in the categories of CRSC and CRSC2. Furthermore, the impact of post-fire climate on the recovery process was also shown in the mid-term and long-term. In the case of the Yeste Fire, an extreme drought event in 2005 coincided with the breakpoints in NCR and CRSC, whereas in the Requena fire this event was not as severe as it was in 2012, coinciding with the breakpoints in NCR and CRSC2.

Fire severity was also a key factor in the short-term recovery of pine forests. Some studies also found that fire severity was decisive in recovery dynamics both in mixed-conifer forests [16] and pine forests [18] due to its relation to pine seedling densities after fire, depending on pre-fire vegetation composition, seedling mortality and reestablishment processes [8,53,56]. In this study, the areas were burned at high severity (dNBR > 0.66) as forested areas tended to experience a higher severity compared to herbaceous and shrublands. Moderate-high severity was only found in those pixels of CR, which showed a slower but stable recovery trend. In agreement with the results reported by Shvetsov et al. (2019) [21], a positive relationship between recovery rate and fire severity was found, since recovery rates were higher for the higher severity areas than for high-moderate severity (corresponding to successful recovery). Moderate fire severity sites might result in higher soil organic matter mineralization, and thus in higher post-fire soil fertility that produces faster growth in pine seedlings [89]. Bright et al. (2019) [28] also found that areas burned at higher severities recovered at faster rates, possibly because they are fire-adapted forest in which fire creates favorable conditions for vegetation germination and regeneration. Nevertheless, some studies in Mediterranean pine forest reported that conversion from forest to shrubland occurred in the most xeric sites (south-facing areas) [48] or in those areas with a high severity [12]. In this sense, Baudena et al. (2019) [90] predicted that future potential increases in aridity may drive these fire-prone ecosystems past a tipping point, after which closed forest structure would be replaced by open shrublands.

Topographic variables can also influence post-fire vegetation recovery through its effects on local microclimate, soil and hydrological processes [18,20]. Wittenberg et al. (2007) [91] and Ireland and Petropulous (2015) [18] found that north facing aspects exhibit higher rates of vegetation recovery compared to south facing aspects as we found in the short-term recovery in CRS, NCR and CR categories according to TCA, but also in CRSC according to TCW. The negative influence of elevation was also detected in Mediterranean pine forest [37] and red fir forests [16] that might be attributed to the decreased temperature with elevation. In contrast, Chu et al. (2017) [20] and Shvetsov et al. (2019) [21] reported that topographic variables were the least important factors in explaining the regeneration rate in Siberian forests. In our study, recovery in relation to topographic position did not show any clear pattern. This could be due to the fact that most of the pixels were located at either upper or mid-slope in the foothills and very few were bottom slopes.

Although we found that post-fire climate was the most important variable in explaining post-fire recovery in the mid and long-term, other variables not included could be influencing the recovery process. Further analysis would consider the influence of historical management legacies as well as the distance to seed banks in the post-fire recovery patterns.

5. Conclusions

Time series analysis from a temporal segmentation approach allowed us to unravel and characterize different post-fire recovery trajectories. Although several studies have addressed the estimation of post-fire recovery rates, fewer have been done in defining and characterizing the differences among recovery trends in Mediterranean pine forests. Here, we identified different recovery categories according to TCA trajectories and TCW trajectories, which enabled us to define slow but more stable recovery processes (CR and CR2) compared to other faster but also interrupted recovery processes (CRSC, CRSC2, CRS and NCR). The appraisal of the environmental and contextual drivers of the recovery process showed that fire severity is important to predict the RR in the short-term but post-fire climate in terms of drought better explained the RR in the mid and long-term.

The thermophilous pine forests are the most affected by wildfires in Europe. Increased wildfire activity is expected to continue under warmer and drier conditions, making post-fire vegetation recovery of concern to researchers and forest managers. Since these forests may not be allowed the time to develop into a mature forest that would be able to recover rapidly, the resilience of these ecosystems will therefore be significantly reduced. Hence, a better understanding of fire regimes and forest recovery patterns in different environmental and climatic conditions is needed for developing forest management strategies that enhance forest resilience.