Neural-Network Time-Series Analysis of MODIS EVI for Post-Fire Vegetation Regrowth

The time-series analysis of multi-temporal satellite data is widely used for vegetation regrowth after a wildfire event. Comparisons between preand post-fire conditions are the main method used to monitor ecosystem recovery. In the present study, we estimated wildfire disturbance by comparing actual post-fire time series of Moderate Resolution Imaging Spectroradiometer (MODIS) enhanced vegetation index (EVI) and simulated MODIS EVI based on an artificial neural network assuming no wildfire occurrence. Then, we calculated the similarity of these responses for all sampling sites by applying a dynamic time warping technique. Finally, we applied multidimensional scaling to the warping distances and an optimal fuzzy clustering to identify unique patterns in vegetation recovery. According to the results, artificial neural networks performed adequately, while dynamic time warping and the proposed multidimensional scaling along with the optimal fuzzy clustering provided consistent results regarding vegetation response. For the first two years after the wildfire, medium-highto high-severity burnt sites were dominated by oaks at elevations greater than 200 m, and presented a clustered (predominant) response of revegetation compared to other sites.


Introduction
Monitoring vegetation regrowth following a wildfire is essential for forest managers to assess the regeneration capacity of burnt lands [1].Forest ecologists are particularly interested in the effects of fire on seedling emergence, survival, and growth of the main species involved in plant regeneration and the changes in the relative abundance of species [2].Understanding regeneration dynamics are crucial for Mediterranean ecosystems due to repeated wildfires, followed by human activities on the affected areas that limit the ecosystem's capacity to naturally regenerate.Therefore, research efforts focus on understanding the link between burn severity and other first-order fire effects on vegetation recovery.The ability to accurately estimate vegetation recovery potential can aid in fire rehabilitation efforts to focus efforts on burnt areas where human intervention is required to prevent undesired vegetation succession.
A significant number of studies explored the science of remote sensing to extract data that can be used for modeling vegetation regrowth.Several studies focus on species-specific vegetation regrowth in different types of ecosystems such as conifers [3,4], mixed conifer/oaks [5], chaparral/maquis forests [6][7][8], and subalpine forests [9].Most published research uses satellite images from the Landsat series [4,6,7] and Moderate Resolution Imaging Spectroradiometer (MODIS) [8, 10,11] sensors due to their high temporal resolution and public access to archived data; however, radar and multisensor approaches were also applied successfully [5,12,13].
Different approaches are also used for time-series analysis of remote-sensing data.Some methods use temporal trajectory diagrams and compare ground-truth sampling data or estimates of burn severity with unburnt control sites [16].These comparisons are either visual or statistical between burnt and unburnt sites using the average differences [17].However, the same EVI does not imply the presence of the same vegetation type, especially in Mediterranean mixed conditions and vice versa, i.e., the same vegetation type does not mean equal EVI.Furthermore, EVI also depends on vegetation status, which means variations in solar radiation and weather conditions (temperature and precipitation) influence the evapotranspiration [18].Due to the above uncertainties, similarity of the burnt and an unburnt control site cannot be assured.Other studies estimate the temporal trend of post-fire vegetation regrowth by fitting indices using logarithmic and polynomial models and extracting the recovery rate by interpreting the slope coefficient [19] or the Theil-Sen slope estimator [20].Pickell et al. [21] and White et al. [22] developed an index named "years to recovery" (Y2R) to examine when a pixel will attain 80% of the mean spectral value it had two years prior to wildfire.Moreover, parametric (multiple linear regression (MLR)) and non-parametric tests (classification and regression trees (CARTs)) were also applied to analyze the effects of fire severity, terrain restoration, and pre-fire fuel treatments on post-fire vegetation response.The temporal analysis of each study varied from a few months [23][24][25] to several years [21,26,27].
This study examines whether burn severity, vegetation type, and elevation influence post-fire vegetation regrowth by combining a time series of MODIS-based enhanced vegetation index with post-fire ground-truth estimation for a burnt area in northern Greece.The EVI time series covers a post-wildfire period (2011-2017) and was combined with pre-fire data to examine "what if" scenarios of the potential vegetation condition in the absence of wildfire.The difference between the real and predicted post-fire time series was calculated, and we clustered the sampling points based on these differences to distinguish the effects of (i) burn severity, (ii) vegetation type, and (iii) elevation on vegetation recovery for a short-(two years) and a long-term (seven years) period.This study incorporates an artificial neural network (ANN) as a nonlinear autoregressive (NAR) model approximator to analyze the time-series prediction of the post-fire EVI ("what if" scenario).We applied the dynamic time warping (DTW) algorithm to estimate the similarity between the real and the predicted post-fire time series.We used multidimensional scaling to transform the DTW distances into points in a low-dimensional space and used optimal fuzzy clustering to partition the above points into meaningful groups.

Study Area and Field Data Collection
Our study area is in the northeast (NE) region of Greece (Evros), mostly covered by dense forests and agricultural lands.The Evros wildfire started on 24 August 2011 and, under strong northeastern winds, rapidly expanded southwest (Figure 1) [28].The fire was contained after three days, resulting in over 5900 ha of burnt conifer forests, oak stands, shrublands, grasslands, and agricultural areas in low-elevation areas with gentle slopes (between 120 and 385 m from 5 to 25 • , respectively).The vegetation of the study area mainly consists of conifer stands (Pinus brutia) with understory vegetation composed of shrubs and broadleaf trees (Phillyrea media, Quercus coccifera, Laurus nobilis, Acer spp., Juniperus oxycedrus, Erica malipuliflora, and Arbutus andrachne).A secondary vegetation type is oak-dominated forest (Quercus pubescens and Q. conferta) or mixed with conifers.During the last 30 years, there was a net increase of 33% in forest coverage from forest clearings and agricultural lands [29].A large part of the historical burnt area is covered by conifer reforestations established 50 years ago by the Forest Service.The area was dominated by shrubs mixed with oak which were removed by mechanical means to be replaced by conifer species planted in human-made terraces.Fuel breaks, agricultural lands, and vegetation-free zones helped fire containment efforts and the preservation of "green" forested enclaves, creating a mosaic of burnt and unburnt patches.The Evros fire burned into a peripheral part of an important protected area for bird species (i.e., Dadia National Park) where 36 out of 38 European raptor species reside, some of them very rare (Aquila heliacal and Aquila pomarina), along with 219 other bird species, 40 reptile species, and 48 mammal species [30][31][32][33][34].
ISPRS Int.J. Geo-Inf.2018, 7, 420 3 of 30 vegetation of the study area mainly consists of conifer stands (Pinus brutia) with understory vegetation composed of shrubs and broadleaf trees (Phillyrea media, Quercus coccifera, Laurus nobilis, Acer spp., Juniperus oxycedrus, Erica malipuliflora, and Arbutus andrachne).A secondary vegetation type is oak-dominated forest (Quercus pubescens and Q. conferta) or mixed with conifers.During the last 30 years, there was a net increase of 33% in forest coverage from forest clearings and agricultural lands [29].A large part of the historical burnt area is covered by conifer reforestations established 50 years ago by the Forest Service.The area was dominated by shrubs mixed with oak which were removed by mechanical means to be replaced by conifer species planted in human-made terraces.Fuel breaks, agricultural lands, and vegetation-free zones helped fire containment efforts and the preservation of "green" forested enclaves, creating a mosaic of burnt and unburnt patches.The Evros fire burned into a peripheral part of an important protected area for bird species (i.e., Dadia National Park) where 36 out of 38 European raptor species reside, some of them very rare (Aquila heliacal and Aquila pomarina), along with 219 other bird species, 40 reptile species, and 48 mammal species [30][31][32][33][34].Following the wildfire in 2011, 27 sampling field plots were established to estimate burn severity and measure key post-fire stand variables (Figure 1), including tree species composition, and an initial estimation of tree survival rate, diameter at breast height, tree height, and crown base height, as well as scorch height, and canopy cover.Within the selection of sampling field points, we assured that all the possible combinations of burn severity, vegetation cover, and elevation would be covered while all the sites were accessible by the road network.That was assured by conducting stratified random sampling based on an initial mapping of fire severity through an early assessment of difference NBR (dNBR).It should be noticed that the first sampling was conducted three weeks after the fire; hence, regeneration or immediate change in growing stages was not observed.To measure fire severity, we applied the composite burn index (CBI) [14] and the Ryan and Noste [35] methods (Figure 2a).CBI entails a relatively large plot, independent severity ratings for individual strata, and a synoptic rating for the whole plot area [14].The Ryan and Noste [35] method estimates flame length Following the wildfire in 2011, 27 sampling field plots were established to estimate burn severity and measure key post-fire stand variables (Figure 1), including tree species composition, and an initial estimation of tree survival rate, diameter at breast height, tree height, and crown base height, as well as scorch height, and canopy cover.Within the selection of sampling field points, we assured that all the possible combinations of burn severity, vegetation cover, and elevation would be covered while all the sites were accessible by the road network.That was assured by conducting stratified random sampling based on an initial mapping of fire severity through an early assessment of difference NBR (dNBR).It should be noticed that the first sampling was conducted three weeks after the fire; hence, regeneration or immediate change in growing stages was not observed.To measure fire severity, we applied the composite burn index (CBI) [14] and the Ryan and Noste [35] methods (Figure 2a).CBI entails a relatively large plot, independent severity ratings for individual strata, and a synoptic rating for the whole plot area [14].The Ryan and Noste [35] method estimates flame length classes, derived from direct observations or inferred from post-burn observations and reconstruction of the fire environment, in addition to the depth of char classification, derived from post-burn observations of the extent to which fuels were burnt, particularly on the soil surface.One year after the fire (September 2012), new fieldwork was conducted on the existed field plot locations, where we measured regeneration (number of seedlings and re-sprouted shrubs), late survival rate, and soil condition (Figure 2b).classes, derived from direct observations or inferred from post-burn observations and reconstruction of the fire environment, in addition to the depth of char classification, derived from post-burn observations of the extent to which fuels were burnt, particularly on the soil surface.One year after the fire (September 2012), new fieldwork was conducted on the existed field plot locations, where we measured regeneration (number of seedlings and re-sprouted shrubs), late survival rate, and soil condition (Figure 2b).

Satellite Data Retrieval and Preprocessing
EVI is a modified NDVI that includes adjustments for canopy background and aerosol effects, and it is widely used in phenological studies [36][37][38][39].For time-series modeling, it is required to have high-temporal-resolution periodic data (i.e., one or two images per month), accessing a significant amount of historic data.The most widely used sensors that fulfill these prerequisites are the MODIS, Advanced Very High Resolution Radiometer (AVHRR), and Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER).From these three sensors, only MODIS provides a blue channel [40], and the calculation of the EVI can be applied using Equation (1) [41].
2.5 6 7.5 1 where  is the reflectance for each corresponding channel.We retrieved EVI (MYD13Q1 product- In Figure 3, we show the EVI for the month of September one year before the wildfire and for the following six years.

Satellite Data Retrieval and Preprocessing
EVI is a modified NDVI that includes adjustments for canopy background and aerosol effects, and it is widely used in phenological studies [36][37][38][39].For time-series modeling, it is required to have high-temporal-resolution periodic data (i.e., one or two images per month), accessing a significant amount of historic data.The most widely used sensors that fulfill these prerequisites are the MODIS, Advanced Very High Resolution Radiometer (AVHRR), and Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER).From these three sensors, only MODIS provides a blue channel [40], and the calculation of the EVI can be applied using Equation (1) [41].
where ρ is the reflectance for each corresponding channel.We retrieved EVI (MYD13Q1 product- In Figure 3, we show the EVI for the month of September one year before the wildfire and for the following six years.

Time-Series Prediction of Post-Fire EVI
To evaluate the influence of wildfire on post-fire vegetation recovery, we estimated an EVI time series for a scenario where the 2011 wildfire did not occur.Various time-series models are applied in remote-sensing vegetation indices such as the linear model autoregressive integrated moving average (ARIMA), or variations developed for dealing with non-stationary data, such as the seasonal ARIMA (SARIMA) [43].EVI is a nonlinear, nonstationary, and seasonal time series, and the ANN is preferable to other nonlinear time series such as bilinear models, threshold autoregressive models, and regression trees [44].In our time-series prediction, a nonlinear autoregressive model (NAR) was applied.In this case, the data have only one series and the future value y(t) is predicted using d past values according to Equation (2).
where f is a nonlinear known function that was estimated with a multilayer perceptron (MLP).For the training of the ANN, the Levenberg Marquardt (LM) algorithm was applied, while the performance was evaluated with the mean square error (MSE).The hyperbolic tangent sigmoid function transferred the output signal from the input to the hidden layer, while a linear function was used for the activation of the output layer.LM is one of the most efficient algorithms for small and medium-sized networks with increased convergence speed [45].Before the training, the input data were smoothed by applying the Savitzky-Golay filter [46,47], where a polynomial was fitted through a fixed number of points.Due to the fact that pixels of maximum value composite vegetation indices from MODIS may contain clouds and other disturbances, using these noisy data will affect the analysis.To cope with this problem, previous efforts proposed the Savitzky-Golay filter to smooth the noise from NDVI time series caused primarily by cloud contamination and atmospheric variability [48].After various tries, a third-degree polynomial was used to compute the smoothed value based on the 10 adjacent data points.Then, the 171 pre-fire input data were divided into two continuous subsets (training and validation), including 80% and 20% of the original dataset, respectively.For each of the 27 sampling points, an ANN was prepared for the prediction of the EVI future values.Training of the ANNs was based on the assumption that the future value of EVI is predicted using the values of the previous 12 months (one value per month).Consequently, the feedback delays were set to 23 with a step of two.An open-loop network was trained after 10 random initializations for one to four hidden nodes, and the network with the best validation performance was saved.Then, a feedback connection was included in the best network for predicting the 151 post-fire EVI values in the absence of wildfire.

Post-Fire Response Comparison
The influence of the wildfire was calculated by subtracting the actual time series of the EVI from the simulated EVI that was estimated with the ANN.The resulting errors (i.e., differences) for the sampling points generated a new time series.Our aim was to compare these errors across all points to evaluate the response and to examine whether any clusters were present.The comparisons were performed by applying the dynamic time warping (DTW) similarity algorithm [49].With this technique, the best alignment between two time series was determined by handling time deformations and speeds in time-dependent data (Figure 4).In comparison to Euclidean distance, DTW takes into consideration the misalignment of peaks and dips between two sequences; therefore, the distance is computed based on the actual corresponding points using Equation (3).
where S = {s 1 , s 2 , ...s i , ..., s n }, T = {t 1 , t 2 , ...t i , ..., t m }, ∆(w k ) is the sample error, δ ST is the DTW distance between S and T time series, and W = w 1 , w 2 , ..., w k is the warping path that aligns the elements of S and T by minimizing the distance between them.DTW is rarely used for remote-sensing time-series applications [50][51][52].In the present study, DTW was calculated for two time periods-two years and seven years post-fire-in order to study the short-and mid-term response of the vegetation.That is, this analysis refers to a two-and a seven-year period, and not to vegetation status after two or seven years.
ISPRS Int.J. Geo-Inf.2018, 7, 420 7 of 30 the elements of S and T by minimizing the distance between them.DTW is rarely used for remotesensing time-series applications [50][51][52].In the present study, DTW was calculated for two time periods-two years and seven years post-fire-in order to study the short-and mid-term response of the vegetation.That is, this analysis refers to a two-and a seven-year period, and not to vegetation status after two or seven years.

Clustering of DTW Distances
Multidimensional scaling was used to transform the set of DTW distances into a set of points in a low-dimensional space, which were then clustered using optimal fuzzy clustering.Assume that we are given a set of M pairwise distances l or entities.The multidimensional scaling attempts to transform the above distances into points in a low-dimensional space [53].This mapping is established in a certain way, according to which the abovementioned distances are translated into distances between the points lying in the lowdimensional space [54].This procedure enables the use of standard techniques, such as cluster analysis, to discover the underlying distribution of the original data.The dimension of the lowdimensional space is denoted as p .A typical value for this parameter, also used here, is In this paper, the aforementioned transformation was carried out using a multidimensional scaling approach called metric least-squares scaling [56].The points that lie in the p -dimensional transformed space are denoted as the objective function of the metric least squares scaling can be calculated using Equation ( 5) [56]. where is the matrix of the dissimilarity degrees, and the matrix of the points' coordinates in the transformed Euclidean space.The above function can be minimized with respect to the elements of X using Equation ( 6) [56].

Clustering of DTW Distances
Multidimensional scaling was used to transform the set of DTW distances into a set of points in a low-dimensional space, which were then clustered using optimal fuzzy clustering.Assume that we are given a set of M pairwise distances δ l ( = 1, 2, ..., M; l = 1, 2, ..., M; l = ) between objects or entities.The multidimensional scaling attempts to transform the above distances into points in a low-dimensional space [53].This mapping is established in a certain way, according to which the abovementioned distances are translated into distances between the points lying in the low-dimensional space [54].This procedure enables the use of standard techniques, such as cluster analysis, to discover the underlying distribution of the original data.The dimension of the low-dimensional space is denoted as p.A typical value for this parameter, also used here, is p = 2 [53 -55].
In this paper, the aforementioned transformation was carried out using a multidimensional scaling approach called metric least-squares scaling [56].The points that lie in the p-dimensional transformed space are denoted as χ = χ 1 χ 2 ... χ p with 1 ≤ ≤ M. By defining the Euclidean distance, the objective function of the metric least squares scaling can be calculated using Equation ( 5) [56].
where ∆ = [δ l ] is the matrix of the dissimilarity degrees, and X = χ j the matrix of the points' coordinates in the transformed Euclidean space.The above function can be minimized with respect to the elements of X using Equation ( 6) [56].
where t is the iteration number, η is a factor that takes small positive values and, To summarize, the above process created a set of points {χ 1 , χ 2 , ..., χ M } in a p-dimensional Euclidean space, each of which corresponded to a certain classifier.
In the next step, optimal fuzzy clustering was applied to group the points into a number of clusters c that were well separated and compact.Given the dataset, X = {χ 1 , χ 2 , ..., χ M }, the objective was to minimize the objective function [57] under the constraint where {v 1 , v 2 , ..., v c } is the set of the cluster centers, u i (χ ) is the membership degree of the vector χ to the i-th fuzzy cluster, and m ∈ (1, ∞) the fuzziness parameter.The cluster centers and membership degrees that solve the abovementioned optimization problems are respectively given as and Equations ( 11) and ( 12) were iteratively applied until convergence.To perform optimal fuzzy clustering and to generate well separated and compact clusters, the Xie-Beni index was applied [58]: The above index was applied with respect to different numbers of clusters varying between 2 ≤ c ≤ c max , where c max << M. The best result corresponded to the number of clusters that yielded the minimum value of S XB .
The workflow of the present study is outlined in Figure 5.
The above index was applied with respect to different numbers of clusters varying between The workflow of the present study is outlined in Figure 5.

Neural Network Training and Validation
Table 1 presents the results of the CBI estimation and the EVI time series training for the 27 field sampling sites.The MSE for the validation datasets was quite low for all points, which meant that the ANN successfully predicted the EVI time sequence for all sites with MSEs ranging between 0.00025 and 0.00170.Linear correlation between MSE and the explanatory variables (CBI, elevation, and landcover type) revealed that MSE was not correlated with any of these variables, meaning that the success of ANN training did not depend on burn severity, elevation, or vegetation.
Furthermore, the constructed networks required one to four hidden nodes for successful training with good generalization and without any over-fitting (Table 1).The R 2 values for the linear regression between the actual and the predicted outputs were also quite high (>0.75)for all sites, except for one (identifier (ID) 11), meaning that the predicted and actual outputs were correlated.
Training success was also confirmed through visual inspection of the output training plots (Figure A1). Figure 6 presents the output training plots for four representative points of the four different CBI classes.The training and the validation outputs (blue and green crosses, respectively) are laid on the actual time series line (blue line), while the predicted time series post-fire (red line) shows a similar pattern to that pre-fire.By subtracting the actual time series after the fire (purple line), we computed the "error", i.e., the influence of the fire, represented by the green line.The low CBI point (ID 16) presents a clearer agreement between actual and predicted series; thus, the error line fluctuated around zero.On the contrary, the high CBI point (ID 8) shows a high difference

Neural Network Training and Validation
Table 1 presents the results of the CBI estimation and the EVI time series training for the 27 field sampling sites.The MSE for the validation datasets was quite low for all points, which meant that the ANN successfully predicted the EVI time sequence for all sites with MSEs ranging between 0.00025 and 0.00170.Linear correlation between MSE and the explanatory variables (CBI, elevation, and land-cover type) revealed that MSE was not correlated with any of these variables, meaning that the success of ANN training did not depend on burn severity, elevation, or vegetation.
Furthermore, the constructed networks required one to four hidden nodes for successful training with good generalization and without any over-fitting (Table 1).The R 2 values for the linear regression between the actual and the predicted outputs were also quite high (>0.75)for all sites, except for one (identifier (ID) 11), meaning that the predicted and actual outputs were correlated.
Training success was also confirmed through visual inspection of the output training plots (Figure A1). Figure 6 presents the output training plots for four representative points of the four different CBI classes.The training and the validation outputs (blue and green crosses, respectively) are laid on the actual time series line (blue line), while the predicted time series post-fire (red line) shows a similar pattern to that pre-fire.By subtracting the actual time series after the fire (purple line), we computed the "error", i.e., the influence of the fire, represented by the green line.The low CBI point (ID 16) presents a clearer agreement between actual and predicted series; thus, the error line fluctuated around zero.On the contrary, the high CBI point (ID 8) shows a high difference between actual and predicted outputs for the first 50 time steps, i.e., approximately two years.We also note that the difference between simulated and actual time series had some peaks and dips due to the time lag between the two series.Figure 7 shows the actual versus predicted values for the same four representative points.The comparison showed that the observed versus predicted output values of EVI approximated a 1:1 relationship, an indicator that the ANN models had good predictive capacity for almost all 27 sites (Figure A2).
Next, we compared the influence of the wildfire across sampling sites.This was accomplished by calculating the similarity measure DTW for each pair of the 27 sites.The results are shown in Table A1 for the first two years after the wildfire, and in Table A2 for the seven years following the wildfire.Low DTW values indicated that these sites had similar responses.For example, the lowest DTW value was observed between ID 0 and ID 18 (0.471) two years after the wildfire, meaning that the difference between the post-fire simulated and actual EVI for site ID 0 was similar to the difference between the post-fire simulated and actual EVI for site ID 18.It should be noted that the minimum and maximum values were not presented at the same pairs of points across two and seven years, i.e., for sampling sites that might respond similarly for the first two years, their responses differentiate within a seven-year period.On the other hand, sites with different burn severity presented different responses for the first two years (high DTW); however, within a seven-year evaluation, their responses became more similar (low DTW).This was because vegetation could react differently in the first two years, while different management practices, weather conditions, and other environmental factors differentiate the long-and short-term responses.

Multidimensional Scaling Clustering of DTW Distances
The site with ID 11 was an outlier according to the primary tries of the clustering, which was also verified by the low R 2 of the ANN training; therefore, it was eliminated from the clustering process.The fuzziness parameter was set to 1.5 m  . The estimation of the Xie-Beni index revealed that the optimal number of clusters was three for both time periods (Table 2); i.e., despite the

Multidimensional Scaling Clustering of DTW Distances
The site with ID 11 was an outlier according to the primary tries of the clustering, which was also verified by the low R 2 of the ANN training; therefore, it was eliminated from the clustering process.The fuzziness parameter was set to m = 1.5.The estimation of the Xie-Beni index revealed that the optimal number of clusters was three for both time periods (Table 2); i.e., despite the increased diversity of the properties (burn severity, vegetation, and elevation) of the sampling sites, their response to wildfire did not show the same variety (as also interpreted from Figures 8 and 9 that present the three clusters for the two study periods).More specifically, for the first two years after the wildfire, cluster C1 (red) was well separated from the other two.Moderate-high and high CBI sites covered by oak or other mixed types with oak located at elevations >200 m dominated this cluster, meaning that they represented a response that differentiated them from the other sites.On the other hand, the response within seven years seemed to be more similar for all points, whereby only two points belonged to a third cluster (Figure 9).However, cluster C1 (red) and cluster C2 (green) included sites that were separated.Specifically, cluster C1 was dominated by sites with low and moderate-low CBI located above 200 m.On the other hand, cluster C2 included sites mostly with moderate-high and high CBI located at elevations of 100 m to 300 m.Vegetation seemed not to influence the above separation.In Figure 10, we can see the number of sampling sites for each cluster for the three properties we examined (CBI, land cover, and elevation).It was expected that unclustered sites would follow the same distribution in each histogram based on the proportion of each category to the total number of points.On the other hand, the absence of a category from a cluster or a distribution of a category that does not follow the distribution of the rest of the categories may signify a possible clustered response to wildfire.For example, in Figure 10a, we noticed that cluster C1 did not include any sites with low CBI, while Figure 10c shows that oak sites had a different distribution than conifer and conifer mixed with oak.More specifically, for the first two years after the wildfire, cluster C1 (red) was well separated from the other two.Moderate-high and high CBI sites covered by oak or other mixed types with oak located at elevations >200 m dominated this cluster, meaning that they represented a response that differentiated them from the other sites.On the other hand, the response within seven years seemed to be more similar for all points, whereby only two points belonged to a third cluster (Figure 9).However, cluster C1 (red) and cluster C2 (green) included sites that were separated.Specifically, cluster C1 was dominated by sites with low and moderate-low CBI located above 200 m.On the other hand, cluster C2 included sites mostly with moderate-high and high CBI located at elevations of 100 m to 300 m.Vegetation seemed not to influence the above separation.In Figure 10, we can see the number of sampling sites for each cluster for the three properties we examined (CBI, land cover, and elevation).It was expected that unclustered sites would follow the same distribution in each histogram based on the proportion of each category to the total number of points.On the other hand, the absence of a category from a cluster or a distribution of a category that does not follow the distribution of the rest of the categories may signify a possible clustered response to wildfire.For example, in Figure 10a, we noticed that cluster C1 did not include any sites with low CBI, while Figure 10c shows that oak sites had a different distribution than conifer and conifer mixed with oak.

Discussion and Conclusions
Wildfires have severe a posteriori impacts on tree cover, especially in the forested areas of Mediterranean ecosystems.The recovery of burnt areas facilitated by post-fire forest management practices should be guided by burn severity estimations, pre-existing vegetation types, and other environmental factors to ecologically restore the landscape.Former studies showed that the use of remote-sensing indices provides a cost-effective way of monitoring vegetation regrowth after a wildfire event [59].In this research, we studied whether the various combinations of burn severity, vegetation, and elevation present a clustered response after a wildfire.We relied on EVI retrieved by MODIS, and recovery was analyzed for two time periods: two years and seven years after the event.
Our results demonstrated that, during the first two years, there were dissimilarities in the response based on CBI, land cover, and elevation.Oaks burnt with moderate to high severity at higher elevations responded differently during the first two years compared to other combinations of burn severity, elevation, and land cover.This is consistent with what was found in previous studies, where environmental factors related to elevation, such as precipitation, seemed to control vegetation recovery, especially in the short term.High precipitation increases soil moisture that helps the regrowth of vegetation.Furthermore, oak forests seemed to respond differently compared to other sites and this was also confirmed by recent studies where in a pine-oak forest, oaks regenerated

Discussion and Conclusions
Wildfires have severe a posteriori impacts on tree cover, especially in the forested areas of Mediterranean ecosystems.The recovery of burnt areas facilitated by post-fire forest management practices should be guided by burn severity estimations, pre-existing vegetation types, and other environmental factors to ecologically restore the landscape.Former studies showed that the use of remote-sensing indices provides a cost-effective way of monitoring vegetation regrowth after a wildfire event [59].In this research, we studied whether the various combinations of burn severity, vegetation, and elevation present a clustered response after a wildfire.We relied on EVI retrieved by MODIS, and recovery was analyzed for two time periods: two years and seven years after the event.
Our results demonstrated that, during the first two years, there were dissimilarities in the response based on CBI, land cover, and elevation.Oaks burnt with moderate to high severity at higher elevations responded differently during the first two years compared to other combinations of burn severity, elevation, and land cover.This is consistent with what was found in previous studies, where environmental factors related to elevation, such as precipitation, seemed to control vegetation recovery, especially in the short term.High precipitation increases soil moisture that helps the regrowth of vegetation.Furthermore, oak forests seemed to respond differently compared to other sites and this was also confirmed by recent studies where in a pine-oak forest, oaks regenerated vigorously [60].
In our study area, field re-sampling one year after the fire confirmed our current findings, whereby oaks and shrubs started regenerating vigorously with sprouting (Figure 11).and the fewer short Mediterranean shrubs were present on the understory, the lower the fire intensity was.
With respect to the EVI simulation and cluster analysis, the applied ANN approach used to forecast EVI performed quite well.Overall, the methods utilized in this study can enhance vegetation recovery estimation and assist forest managers in decision-making based on vegetation type, elevation, and burn severity.Prioritizing management options with a long-term wildfire risk mitigation perspective can reduce future catastrophic events.In the case of the Evros wildfire, human efforts to change the predominant vegetation types (oak and shrub) to a productive coniferdominated forest failed, since fire rapidly spread through the revegetated terraces and burned there with higher severity.On the other hand, within a seven-year period after the wildfire, the recovery was not clustered to the same extent.The analysis revealed that vegetation type which existed before the fire event did not influence the amount of the vegetation after the wildfire, as depicted by the EVI time series.However, low-burn-severity sites located at higher elevations responded differently compared with high-burn-severity sites at low elevations.Once again, environmental factors such as rainfall are probably the main cause for these findings.From Figures 6 and A1, it is evident that the EVI values (purple line) after three years tend to approach the maximum values of pre-fire EVI (blue line), regardless of the severity.These results align well with previous studies where MODIS EVI had a major reduction mainly in the most severe burns, while, after 5-8 years, the summer EVI for all severity classes recovered to within 90-108% of pre-fire levels [61].
From Figures 8 and 9, it is evident that the pixels vegetated with conifer mixed with oak did not present any clustering in the recovery, unlike unmixed conifer.Moreover, the mixture conditions did not provide any difficulties in training the ANN (see Table 1).However, Table 1 shows that the unmixed pixels burned with moderate to high severity.Only one sampling point of the unmixed conifer (site with ID 13) burned with low severity.On the other hand, the mixed pixels presented an equal distribution regarding the burn severity.A similar conclusion was reached by Thompson et al. [62] who found that areas that were salvage-logged and planted (unmixed pixels) burned more severely than similar unmanaged areas (mixed pixels).In addition, mixed Mediterranean shrubs with conifers burn at a moderate speed, usually with high intensity.Half of the mixed conifer/oak sampling points burned with low intensity, while the rest burned with high intensity.This was regulated by the number of shrubs on the understory, the age of the stand, the presence of ladder fuel, and the proportions of oak trees in the stands.We noticed that stands with more dominant oak trees managed to survive, since fire intensity and spread usually slowed down when it burned in broadleaf oak litter.The same rule applied to mixed broadleaf with oak stands; the more oak trees and the fewer short Mediterranean shrubs were present on the understory, the lower the fire intensity was.
With respect to the EVI simulation and cluster analysis, the applied ANN approach used to forecast EVI performed quite well.Overall, the methods utilized in this study can enhance vegetation recovery estimation and assist forest managers in decision-making based on vegetation type, elevation, and burn severity.Prioritizing management options with a long-term wildfire risk mitigation perspective can reduce future catastrophic events.In the case of the Evros wildfire, human efforts to change the predominant vegetation types (oak and shrub) to a productive conifer-dominated forest failed, since fire rapidly spread through the revegetated terraces and burned there with higher severity.Future work will include new field measurements to explore how vegetation management and climatic conditions influence actual restoration outcomes and dominant species several years after the wildfire event.Other sensors with high spatial and temporal resolution (e.g., SENTINEL 2A and 2B) will be evaluated for vegetation response following a wildfire, given the extended data archive that is underway after their recent launch.Finally, new time-series models will be explored, such as the long short-term memory (LSTM) networks.Recent studies show that LSTM networks present a robust performance in time series due to the fact that they are capable of learning long-term dependencies [63].

Figure 1 .
Figure 1.Study area and field plots.

Figure 1 .
Figure 1.Study area and field plots.

Figure 2 .
Figure 2. Fieldwork at the study area with (a) high-severity burnt coniferous forest identified after the wildfire, and (b) clear-cut salvage of the dead trees and regrowth views during a first-year revisit.

Figure 2 .
Figure 2. Fieldwork at the study area with (a) high-severity burnt coniferous forest identified after the wildfire, and (b) clear-cut salvage of the dead trees and regrowth views during a first-year revisit.

Figure 3 .
Figure 3. Enhanced vegetation index (EVI) of the study area from the Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua sensor (MYD13Q1) for the month of September, for the years 2010-2017.The Evros wildfire started at the end of August 2011.

Figure 3 .
Figure 3. Enhanced vegetation index (EVI) of the study area from the Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua sensor (MYD13Q1) for the month of September, for the years 2010-2017.The Evros wildfire started at the end of August 2011.

Figure 4 .
Figure 4. Dynamic time warping (DTW) versus Euclidean distance for two arbitrary time series.

Figure 4 .
Figure 4. Dynamic time warping (DTW) versus Euclidean distance for two arbitrary time series.
best result corresponded to the number of clusters that yielded the minimum value of XB S .

Figure 5 .
Figure 5. Flowchart of the analysis conducted in the present study.

Figure 5 .
Figure 5. Flowchart of the analysis conducted in the present study.

Figure 6 .
Figure 6.Output of the artificial neural network (ANN) training plots for sampling sites identifier (ID) 16 (low CBI), ID 3 (moderate-low CBI), ID 0 (moderate-high CBI), and ID 8 (high CBI).The green line reflects the difference between the actual and simulated post-fire enhanced vegetation index (EVI) time series (TS).

Figure 8 .
Figure 8. Clustering of the dynamic time warping (DTW) distances across the 26 sample points for the first two years after the wildfire labeled by (a) CBI, (b) land-cover type, and (c) elevation.(χ1 and χ2 are the dimensions of the two-dimensional transformed space of the DTW distances obtained from the metric least-squares scaling).Axes are dimensionless.

Figure 8 .
Figure 8. Clustering of the dynamic time warping (DTW) distances across the 26 sample points for the first two years after the wildfire labeled by (a) CBI, (b) land-cover type, and (c) elevation.(χ1 and χ2 are the dimensions of the two-dimensional transformed space of the DTW distances obtained from the metric least-squares scaling).Axes are dimensionless.

Figure 9 .
Figure 9. Clustering of the DTW distances across the 26 sample points for the seven-year period after the wildfire labeled by (a) CBI, (b) land-cover type, and (c) elevation.(χ1 and χ2 are the dimensions of the two-dimensional transformed space of the DTW distances obtained from the metric leastsquares scaling).Axes are dimensionless.

Figure 9 .
Figure 9. Clustering of the DTW distances across the 26 sample points for the seven-year period after the wildfire labeled by (a) CBI, (b) land-cover type, and (c) elevation.(χ1 and χ2 are the dimensions of the two-dimensional transformed space of the DTW distances obtained from the metric least-squares scaling).Axes are dimensionless.

Figure 10 .
Figure 10.Distribution of the number of sampling sites for each cluster according to (a) CBI for the first two years, (b) CBI for the first seven years, (c) land cover for the first two years, (d) land cover for the first seven years, (e) elevation for the first two years, and (f) elevation for the first seven years.

Figure 10 .
Figure 10.Distribution of the number of sampling sites for each cluster according to (a) CBI for the first two years, (b) CBI for the first seven years, (c) land cover for the first two years, (d) land cover for the first seven years, (e) elevation for the first two years, and (f) elevation for the first seven years.

Figure 11 .
Figure 11.A first-year revisit of the study area identifying regeneration outcomes at (a) plot ID 15 and (b) plot ID 23.

Figure 11 .
Figure 11.A first-year revisit of the study area identifying regeneration outcomes at (a) plot ID 15 and (b) plot ID 23.

Figure A1 .
Figure A1.Output of the ANN training plots for the 27 sampling points.
Version 6 of MODIS Aqua sensor) by submitting the 27 field sampling points through the Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) platform [42].This product is a 16-day maximum value composite (MVC) EVI, with a spatial resolution of 250 m, for the period of 13 March 2004 to 26 February 2018 (i.e., approximately seven years before and after the wildfire).During quality control, bad quality values according to the supplied metadata were replaced with interpolated values.Therefore, for each sampling point, a time series of 322 EVI values was available.For the pre-fire period (13 March 2004 to 24 August 2011) and the post-fire period (24 August 2011 to 26 February 2018), 171 and 151 EVI values were retrieved, respectively, i.e., one EVI value per 16 days.

Table 1 .
Composite burn index (CBI) estimation and training results for the 27 field sampling sites.ID-identifier; MSE-mean square error.

Table 2 .
Number of optimal clusters of the sampling sites for the post-fire study periods according to the Xie-Beni index.

Table A1 .
Similarity measure (dynamic time warping (DTW)) of the response for each pair of the 27 points for the first two years.

Table A1 .
Similarity measure (dynamic time warping (DTW)) of the response for each pair of the 27 points for the first two years.