Short-term temporal and spatial analysis for post-fire vegetation regrowth characterization and mapping in a Mediterranean ecosystem using optical and SAR image time-series

Abstract In the present study, the temporal and spatial dynamics of the post-fire recovery of different Mediterranean vegetation types during the three years after the fire event were analyzed, according to different fire severity categories, integrating the use of Synthetic Aperture Satellite Radar (SAR) (Sentinel-1) and optical (Sentinel-2) image time series. The results showed that Mediterranean forest species and shrub/herbaceous species are adapted to fire, with high efficiency in restoring the vegetation cover. Differently, the ecological vulnerability of non-native eucalyptus plantations was found in a lower recovery trend during the observation period. The use of optical short-wave infrared (SWIR) and SAR C-band-based data revealed that some ecological characteristics, such as the woody biomass and structure, recovered at slower rates, comparing to those suggested by using near-infrared (NIR) and red-edge data. An optimized burn recovery ratio (BRR) was proposed to estimate and map the spatial distribution of the degree of vegetation recovery.


Introduction
Mediterranean ecosystems are critical at the European level due to their high biodiversity and intense interaction with human activities (Aragones et al. 2019;Moreira et al. 2020). The typical climatic regime of the Mediterranean region, characterized by long periods of summer drought, entails an increasing wildfire risk in terms of frequency, extension, and severity (Chuvieco 2009;Moreira et al. 2020). This latter factor, defined as the degree of environmental chemical-physical alterations, decomposition, and loss of above/belowground organic matter caused by fire (Key and Benson 2006;Keeley 2009), plays a critical role in the sustainability of Mediterranean habitats, influencing the competitive interactions between species and their post-fire recovery capability Frazier et al. 2015;Fern andez-Garc ıa et al. 2018;Morresi et al. 2019;Lacouture et al. 2020). Moreover, the widespread accumulation of woody fuel, consequent to the abandonment of the semi-natural and agricultural areas that have been occurring in Mediterranean territories, causes suitable conditions for the ignition and the progress of such high intensity/severity fires (Xofis et al. 2022). After a forest fire, the spatial and temporal analysis of both the fire severity and the vegetation response and recovery is a crucial step for assessing the fire impact on ecosystems and their capacity to repristinate the ecosystem services partially lost due to the fire (Ryu et al. 2018;Semeraro et al. 2019;Huang et al. 2020). The post-fire recovery of vegetation cover structure promotes the restoration of numerous ecosystem services, such as the carbon sequestration induced by the regenerating process of forest vegetation, which mitigates the carbon emissions to the atmosphere due to fire (Frolking et al. 2009;Morresi et al. 2019;Semeraro et al. 2019;Huang et al. 2020).
The efficiency of remotely sensed data applied to wildfire assessment has increased significantly in the last decade, thanks to the availability of satellite imagery at increasing spatial, spectral, and temporal resolutions (Gitas et al. 2012;Chu and Guo 2013;Mallinis et al. 2018;Chuvieco et al. 2019;Garc ıa-Llamas et al. 2019). Remote sensing techniques based on time-series of optical vegetation indices (VIs) have been widely used for postfire analysis and monitoring (Gouveia et al. 2010;Gitas et al. 2012;Chu and Guo 2013;Frazier et al. 2015;Chompuchan and Lin 2017;Meng et al. 2018;Ryu et al. 2018;Morresi et al. 2019;Semeraro et al. 2019;Huang et al. 2020). VIs maximize the sensitivity to plant biophysical characteristics and reduce the effects of atmosphere, surface topography and soil variability (Morresi et al. 2019;Xofis et al. 2022). The normalized difference vegetation index (NDVI) is widely employed to detect and examine the post-fire vegetation recovery using a multi-temporal approach (Goetz et al. 2005;Mitri and Gitas 2013;;Polychronaki et al. 2014;Fern andez-Garc ıa et al. 2018;Ryu et al. 2018;Morresi et al. 2019;Semeraro et al. 2019;Lacouture et al. 2020). Other VIs based on shortwave infrared (SWIR) bands are also applied for long-term post-fire vegetation monitoring (Epting et al. 2005;Chen et al. 2011;Kennedy et al. 2012;Chompuchan and Lin 2017;Ryu et al. 2018;Morresi et al. 2019;Semeraro et al. 2019;Grabska et al. 2020). Semeraro et al. (2019) showed that integrating NIR and SWIR data improved vegetation water content retrieval. Morresi et al. (2019) also used the NDVI and SWIR-based indices (normalized burn ratio, NBR and normalized difference moisture index, NDMI) to analyze the post-fire recovery dynamics in Mediterranean landscapes; the latter indices are more sensitive to those purposes.
Despite their proven efficiency in fire effects analysis, optical data have some limitations, namely the presence of clouds (Minchella et al. 2009;Morresi et al. 2019;Huang et al. 2020;Lacouture et al. 2020). Also, the spectral confusion between cloud shadows and burned areas creates considerable difficulties in discriminating one from the other (Pereira et al. 1999;Chuvieco et al. 2005). Furthermore, NIR-based VIs can be affected by earlier saturation already at relative low conditions of vigorous photosynthetic activity and growth of the leaf structure due to their high sensitivity to the chlorophyll content and their high correlation with the leaf area index (LAI) (Wang et al. 2005;Minchella et al. 2009;Chompuchan and Lin 2017;Huang et al. 2020). Moreover, the optical waves do not penetrate the canopy, only providing information on the regeneration of the most superficial vegetation layers (Chompuchan and Lin 2017). In the context of post-fire monitoring, this translates into a fast recovery of the VIs values, close to those of pre-fire conditions, corresponding to an overestimation of full recovery of the ecosystem and a not entirely realistic relationship between fire severity and post-fire regrowth dynamics (Wang et al. 2005;Tanase et al. 2011;Meng et al. 2015Meng et al. , 2018Zhou et al. 2019). Some studies (e.g., Ryu et al. 2018;Morresi et al. 2019), analysing long-term post-fire spectral dynamics, observed that SWIR-based VIs express a more gradual temporal recovery rate than NDVI. The latter was 1-5 years earlier due to its greater sensitivity to photosynthetically active vegetation, combining the red and NIR bands. However, other scholars stated that spectral optical signals remain coherent only when restricted to the first decade after the disturbance Kennedy et al. 2012;Frazier et al. 2015).
In this context, active synthetic aperture radar (SAR) can integrate optical information in vegetation recovery analysis and mapping (Minchella et al. 2009;Tanase et al. 2015;Martins et al. 2016;Chen et al. 2018). Its high sensitivity to the structural properties of the vegetation, with a generally linear correlation between backscatter and vegetation biomass (Quegan et al. 2000;Saatchi et al. 2012;Martins et al. 2016;Yu and Saatchi 2016;Chen et al. 2019;Saatchi 2019), and its capabilities for all-weather and solar radiation independency, make the SAR backscatter complementary information with optical data (Minchella et al. 2009;Tanase et al. 2011;Polychronaki et al. 2014;Zhou et al. 2019). SAR data, however, has its challenges. The complex interactions between the backscatter and scattering components of the soil and vegetation affected by the fire are influenced by intrinsic SAR sensor/signal parameters (e.g., wavelength, polarization, incidence angles, look direction), structural/geometrical (objects structure, amount of scattering elements, surface roughness, geometry and topography of study area) and dielectrics characteristics of the affected surfaces, as well as environmental conditions (e.g., soil moisture, rain, wind) Tanase, Santoro, Wegm€ uller, et al. 2010;Tanase et al. 2014;Chen et al. 2018;Ban et al. 2020;). Among these factors, the wavelength is the one that most influences the ability of waves to penetrate vegetation cover, and thus the type and amount of information that can be derived about the impact of disturbance on it ). SAR shorter wavelengths (e.g., C-band, 3.8-7.5 cm) mainly interact with vegetation scattering elements such as leaves/ needles, twigs and small branches, or herbaceous vegetation; this reduces the capacity to penetrate the regrowing dense forest canopy, meanwhile becoming less sensitive to structural modifications of vegetation strata (Paloscia et al. 1999;Minchella et al. 2009;Tanase et al. 2011;Santi et al. 2017;Chen et al. 2018;Saatchi 2019). Moreover, the different polarizations of the SAR signal affect the interaction with the forest surface.
Immediately after a disturbing event and for the first year after, the scattering effect of the damaged vegetation structure is lacking/decreasing. At the same time, the contribution to the back diffusion by the humidity and the roughness of the exposed soil is higher. In Mediterranean ecosystems, this generally results in a lowering of the cross-polarized signal, interacting with multiple scattering within the forest canopies (volume scattering), and an increase in the co-polarized signal, interacting with small branches, stems and, principally, the ground surface (direct and specular backscatter) (Richards 2009b;Imperatore et al. 2017;Saatchi 2019;. The backscatter, indeed, typically increasing with forest biomass, has been found more directly correlated to above-ground biomass at cross-polarization than co-polarizations (Saatchi et al. 2012;Yu and Saatchi 2016;Saatchi 2019). Both polarizations can be decisive in detecting better the different effects of fire on vegetation (Tanase et al. 2014;Chen et al. 2018). SAR polarimetric indices were generally employed for environment monitoring (Kim et al. 2014;Periasamy 2018;Nasirzadehdizaji et al. 2019;Pipia et al. 2019;Saatchi 2019;Mandal et al. 2020;De Luca, Silva, Oom et al. 2021;. This work builds on previous research based on Sentinel-1 (S1) adapted dual-polarimetric SAR indices (dual-polarization SAR vegetation index, DPSVI; the radar vegetation index, RVI) applied to burned areas detection , with the addition of the radar forest degradation index (RFDI), for burn severity estimation and mapping (De Luca, Silva, Oom, et al. 2021).
Among the several space agencies providing satellites platforms operating both optical and SAR sensors at different spatial, temporal and spectral resolutions , the Copernicus missions by European Space Agency (ESA) include in its fleet both C-band band (centre wavelength of 5.6 cm) SAR (S1), with both cross-(VH) and co-(VV) polarizations, and multispectral (Sentinel-2, S2) sensors, each of which consists of two polar-orbiting satellites (S1A/B and S2A/B, respectively) (ESA Sentinel Homepage 2022). The high spatial and temporal resolution makes the Copernicus constellation particularly suitable for mapping, quantitative-qualitative characterization, and temporal monitoring of the effects of fire on ecosystems; its free availability, moreover, is a fundamental attribute in risk management and monitoring framework Gitas et al. 2012;Verhegghen et al. 2016;Martinis et al. 2017;Tanase et al. 2020;De Luca, Silva, Oom, et al. 2021;.
Although numerous studies are using SAR C-band S1 or/and optical S2 imagery for burned area detection (Donezar et al. 2019;Roteta et al. 2019;Ban et al. 2020;Carreiras et al. 2020;Pulvirenti et al. 2020;Tanase et al. 2020; or burn severity estimation (Fern andez-Manso et al. 2016;Mallinis et al. 2018;Amos et al. 2019;Quintano et al. 2019;Morresi et al. 2022), the contributions concerning their combined use in assessing and monitoring the temporal response of ecosystems to fire effects and the subsequent recovery patterns are very scarce (Evangelides and Nobajas 2020;Han et al. 2021;Zhang et al. 2021); even less in Mediterranean ecosystems.
Understanding the spectral interactions between these patterns and the main factors that influence the damage and recovery processes is also important. Depending on fire severity, type of vegetation and climate conditions, recovery processes can be very heterogeneous, with large changes in forest structure and species composition (Morresi et al. 2019;Lacouture et al. 2020). Considering the post-fire recovery as a homogeneous and predictable process is almost unrealistic, as stated by Morresi et al. (2019). Understanding the quantitative relationship between post-fire vegetation recovery and fire severity allows the assessment of the temporal effects of fire on ecosystem characteristics such as biodiversity, evapotranspiration, carbon cycling, soil chemical and physical properties Chuvieco 2009;Kasischke et al. 2011;Chuvieco et al. 2014;H€ ausler et al. 2018;Meng et al. 2018). Some authors evaluated post-fire recovery processes in relation to the fire severity category (Fern andez-Manso et al. 2016;Martins et al. 2016;Viana-Soto et al. 2017;Meng et al. 2018;Ryu et al. 2018). However, few studies have analysed both S1 and S2 data on different fire severity categories and different types of forest cover in the Mediterranean region.
The spatial assessment of the post-fire recovery temporal rates was also explored by several authors proposing and applying recovery indices to both optical (Lin et al. 2004;Chou et al. 2009;Chompuchan and Lin 2017;Meng et al. 2018;Ryu et al. 2018;Morresi et al. 2019) and SAR (Minchella et al. 2009). Some of these indices considered the effects of delayed vegetation mortality (Chompuchan and Lin 2017), of annually varying meteorological effects (Ryu et al. 2018), or inter-annual modifications of natural phenological cycles (Morresi et al. 2019) by adding approximate parameters based on the generalization of the surrounding unburned vegetation. However, further investigations are needed given the increased availability of better resolution satellite images, such as S1 and S2, and performant open-source prediction algorithms based on machine learning.
The present study aims to: Assess post-fire forest recovery dynamics and their spatial patterns in a Mediterranean ecosystem, using S1 and S2 spectral vegetation indices; Assess the relationship between the recovery rate, the type of vegetation, and the fire severity level, by also taking into consideration the climatic conditions; Estimate and map the spatial distribution of recovery degree through the burn recovery ratio (BBR), based on pre-and post-fire conditions and optimized through machine learning regressors.

Study area
The study area is the Serra de Monchique mountain chain, located in the Algarve region, south of Portugal (37 18'N; 08 30'W), part of which is included in the European Natura 2000 network as a Special Area of Conservation (SAC) (Natura 2000 Site Code: PTCON0037). The climate is typically Mediterranean, with hot, dry summers and mild humid winters, with oceanic influences given by the relative proximity of the Atlantic Ocean and the heights of the mountain range that intercept the humidity (Martins et al. 2015). Considering the surface occupied, the forest cover is mainly composed of Eucalyptus plantations (Eucalyptus globulus, Labill. 1800) and cork oak forests (Quercus suber L.), part of which consists of the typical semi-natural agro-forestry system (montado in Portuguese), and part consisting of more close forest stands. Typical meso-Mediterranean forest ecosystem with other Mediterranean oaks species (e.g., Quercus ilex L.) and other secondary Mediterranean autochthonous broad-leaves trees are also present. A small part of the forest cover is also composed of isolated areas of Mediterranean conifers (Pinus pinea L., Pinus pinaster Aiton.) (Alves et al. 2007;Catry et al. 2015;San-Miguel-Ayanz et al. 2016;H€ ausler et al. 2018;De Luca et al. 2022). However, the largest part of the study area is covered by heathlands and pastures, composed of typical herbaceous and sclerophyllous shrubby flora alternating with agricultural fields surrounding small urban areas (Mitchell et al. 2009;San-Miguel-Ayanz et al. 2016;De Luca et al. 2022).

Description of the fire event
The fire event (Figure 1) occurred in August 2018 from the 3 rd to the 10 th , affecting 268.9 km 2 , almost entirely represented by herbaceous, shrub and forest areas . For most of the burnt surface (> 50%), the event was of high severity where, regardless of the type of vegetation, only the residues of the burned biomass (ash and coal) remained on the ground, and the bare ground was left exposed. Where the fire had moderate-high severity, the vegetation structure was, however, affected at various levels, with the consumption of the lower layers (grass and shrubs) and a predominant crown fire occurrence, which did not destroy all the canopy structure. In moderate and low severity categories, representing less than 15% of the burned surface, the fire affected the vegetation partially, with part of the canopy killed by heat proximity and the other amount remaining alive (De Luca, Silva, Oom, et al. 2021).

Dataset and pre-processing
Sentinel-1 (S1) dataset and pre-processing The SAR dataset was composed of S1-A/B Level-1 high-resolution ground range detected (GRDH) time-series, acquired in interferometric wide (IW) mode, dual-polarized available: co-polarized VV and cross-polarized VH. The time series covered a period of four years (April 2017-Jun 2021) and, considering that the fire event time occurred in the first days of August 2018, the timeframe was split into a pre-fire (from April 2017 to July 2018) and post-fire (from the second half of August 2018 to June 2021) periods. Moreover, the Long-Term Access policy of the ESA (Copernicus Long Term Archive Access 2022) makes very time-expensive the massive download of images from the official Copernicus Open Access Hub platform. Therefore, the S1 images were downloaded using the Alaska Satellite Facility (ASF) interface (ASF 2022), which also provides the python code for bulk-downloading. The S1 time series comprised 273 images, including ascending (51 pre-fire and 88 post-fire) and descending (49 pre-fire and 85 post-fire) flight paths. The pre-processing of the S1 dataset was carried out using the Sentinel-1 Toolboxes, implemented in the SNAP v.8.0.3 open-source software (ESA SNAP Homepage 2022), and performed via the SNAP-Python interface (Snappy), the access provider to SNAP Java API (ESA SNAP Cookbook 2022). Starting by applying the auto-downloaded orbit information file and the removal of the thermal noise, the SAR data pre-processing involved the radiometric calibration to beta (b0) noughts backscatter standard conventions (Small 2011) and the radiometric terrain correction (RTC) process. RTC consists of the radiometric terrain flattening and the geometric terrain correction of the images using a digital elevation model (DEM) to reduce the geometric and radiometric distortions due to the rough surface topography. The shuttle radar topography mission (SRTM) DEM at 1 arc-second spatial resolution (Farr et al. 2007) was resampled using the bilinear interpolation method (Mendes et al. 2019;; no pixel resampling was instead necessary to the GRDH image products, already available in 10 m x 10 m resampled pixel spacing (ground range x azimuth) (ESA Sentinel-1 User Guide 2016). Subsequently, the stack of all the time series was carried out separately for each of the two flight paths. The geolocation of a master image (automatically chosen by the model among the time series) was adopted in this phase. An image Refined-Lee speckle filter (Lee and Pottier 2009), with a 7 Â 7 pixel window size, was applied to reduce the first speckle-noise effects. Finally, a backscatter monthly time average was computed for each month. The backscatter time average, besides reducing the massive amount of the images compositing the time series, further minimize the adverse effects of speckle noise and environmental variables affecting the SAR signal (Tanase et al. 2015;Zhang et al. 2019;Lapini et al. 2020;).
Sentinel-2 (S2) dataset and pre-processing The optical time series was composed of 253 (60 pre-fire and 193 post-fire) S2-A/B Level-2A (Bottom-Of-Atmosphere, BOA) multispectral images. These images cover the same period of the S1 dataset (April 2017-Jun 2021), excluding the first half of August 2018 (when the fire occurred). Due to the same problem concerning the oldest-acquisition image accessibility from the official Copernicus database (section 2.3.1), the S2 dataset was downloaded from the Google Earth Engine (GEE) collections database and executed through the GEE Python API (Google Earth Engine Guides 2022). The GEE Python API was also employed for S2 image pre-processing before the download, including resampling all S2 bands to 10 m of GSD using the nearest neighbour resampling algorithm. Additionally, the S2-Cloud Probability was used to mask each pixel of the time-series images by cloudiness probability (scaled from 0 to 100). The S2-Cloud Probability mask is available in the GEE Data catalog (GEE Data Catalog: Copernicus S2-Cloud Probability 2022), sampled to 10 m of GSD, generated using the automatic pixel-based sentinel2cloud-detector package (s2cloudless) (Sentinel Hub's cloud detector repository 2022), developed by Sentinel Hub's research team (Sentinel Hub Homepage 2022) and based on the LightGBM machine learning library (LightGBM documentation 2022). Higher values of cloud probability have a higher ability to detect dense clouds or highly reflective surfaces, but the omission of less dense clouds could be equally high. Lower values allow the possibility of detecting all the clouds on the scene; however, they increase the risk of more frequent commission errors due to the confusion between clouds and medium-high reflectance surfaces, masking them from the resulting image. To this end, a pixel probability threshold of 10 was set based on previous experiences in the same study area (De Luca et al. 2022). The masked pixels are replaced by applying a temporal linear interpolation involving all the time-series images, differentiating pre-fire and post-fire sub-sets. After the cloud masking, the images were averaged each month using the criterion adopted for the S1 time series.

Land use land cover (LULC) and fire severity mapping
The temporal analysis of post-fire vegetation recovery was applied to three main LULC classes representing the vegetation of the study area: the eucalyptus (Euc) plantations; the autochthonous forest (AuFor), both natural and the semi-natural, consisting of two dominant species, Quercus suber and Q. ilex, and other secondary broadleaves; heathlands, shrublands or pastures vegetation (Pas/Shr). The pine class was excluded in this study due to its limited representativeness in the study area. The definition of the primary reference LULC classes was carried out using a classified LULC map retrieved from a supervised machine learning-based classification processing developed by De Luca et al. (2022) and based on the combined use of both SAR S1 and optical S2 data. The LULC map has a spatial resolution of 10 m x 10 m, and its overall accuracy is higher than 90%.
The burned vegetation was monitored based on the fire severity. The spatial distribution of fire severity was retrieved from De Luca et al. (2021), obtained applying the random forest (RF) machine learning model on a dataset constituted by an S1 þ S2 dataset and derived indices, and trained using a set of field measurements of the composite burn index (CBI) (Key and Benson 2006). Based on the CBI sampling protocol, six conventional fire severity categories were detected (Key and Benson 2006), and five of them were considered in the further analyses: a) unburned soil/rock (not taken into consideration); b) unburned vegetation; c) low severity: low impact of the fire, which was mainly kept at ground level with low levels of alteration of the shrub and/or tree cover; d) moderate severity: level of alteration of the vegetation higher than in low severity category, with the fire reaching the lower layers of the forest canopy, resulting in a mixture of scorch and green vegetation; e) moderate-high severity: predominance of burnt vegetation with a high percentage of tree foliage affected by scorch, and a part of the woody components of the canopy partially or totally consumed by fire; f) high severity: the short vegetation is consumed, as well as most of the tree's foliage, the surface is mainly covered with ash and charcoal.

Climate variables
Two datasets with monthly climate variables series were retrieved, from April 2017 to January 2021, to analyse the relationship with the temporal profiles of vegetation affected by fire: monthly rainfall (mm) and average monthly temperature ( C). After analysing the meteorological stations' spatial distribution, the data were retrieved from 16 stations distributed over the study area (Sistema Nacional de Informações de Recursos H ıdricos, SNIRH 2022).
Where t represents the corresponding image date (months) constituting the time series (April 2017 -June 2021), and VV and VH represent the respective single-polarized backscatters. Besides the SAR VIs, computed separately for orbit path (RVI_As, RVI_Ds, RFDI_As, RFDI_Ds, DPSVI_As, DPSVI_Ds), the backscatter time-series for the single coand cross-polarization (VH_As, VH_Ds, VV_As, VV_Ds) were also involved to compone the final SAR dataset.
A similar procedure was carried out for each month image of the S2 time-series, for which five spectral indices were calculated: the NDVI (Eq. 4), the green normalized vegetation index (GNDVI) (Eq. 5), the normalized red-edge vegetation index (NDRE) (Eq. 6), the normal burn index (NBR) (Eq. 7) and the normalized difference water index (NDWI) (Eq. 8).
B3, B4, B5, B6, B8, B11 and B12 represent the S2 bands conventionally named (ESA Sentinel Homepage 2022). The final optical dataset was composed of the monthly time series of the VIs above.

ROIs collection and temporal profiles extraction
The temporal analysis of vegetation dynamics affected by the fire was divided according to the LULC class and the fire severity category. A series of square 3 Â 3 pixel regions of interest (ROIs) were collected to fulfil these criteria. Each ROI was homogeneous regarding LULC class and fire severity; the whole dataset covered the spatial distribution of the land cover and fire severity. The ROIs were retrieved by visual assessment of LULC and fire severity maps and the Esri ArcGIS World Imagery high-resolution satellite map (Esri ArcGIS World Imagery 2022). A total of 700 ROIs were collected and distributed: 50 ROIs for each fire severity category for both forest LULC classes and the high severity category of Past/Shr class; 25 ROIs for each low, moderate and moderate-high fire severity category for Past/Shr LULC. The choice of 25 ROI was due to the difficulty of finding homogeneous areas of Past/Shr vegetation belonging to these three severity categories. Using the mean value of ROIs allows to examine the overall change occurring inside these sampling units through the entire observation period and compare the different profiles, avoiding the influence of the values of single pixels (Lacouture et al. 2020).

Correlation between temporal profiles and climate variables
The monthly temporal profiles of rainfall and temperature were constructed and compared to the temporal profiles of each S1 and S2 spectral index to investigate the influence that the climate variables had on the dynamics of post-fire vegetation recovery. To define the relationship with the vegetation dynamics of post-fire recovery, we implemented a Pearson's correlation analysis (r) between the temporal profiles of each S1 and S2 index and the climate variables.

The burn recovery ratio (BRR)
The spatial distribution of the degree of vegetation recovery was estimated and mapped using the burn recovery ratio (BRR) (Eq. 9), proposed by Chompuchan and Lin (2017), in turn inspired by the recovery index of Chou et al. (2009) and Lin et al. (2004), which represents the ratio between the value of a given VI at the time of assessment (VI post-fire ) and its value at the reference pre-fire time (VI R_pre-fire ).
As further explained below, VI post-fire and VI R_pre-fire values falling inside the burned perimeter were predicted by the mean of a regression algorithm to minimize the interannual biases.
Due to the natural inter-annual variations that may have occurred in unburned vegetation, if an image taken in the pre-fire period was used for BRR calculation, it would lead to biases for each year of its estimation. Morresi et al. (2019), for example, exploiting the concept already used in Miller and Thode (2007) for the computation of the relative NBR, integrated their recovery index with a coefficient calculated by averaging the difference between pre-fire median and annual post-fire VI values of unburned vegetation to account the inter-annual phenological changes. In this study, to cope with this aspect, the unburned reference values for each index and each year of BRR assessment were estimated using the RF machine learning regression model (Breiman 2001). To this end, 250 additional sample ROIs (ROIs BRR ) (3 Â 3 pixels) were selected for each of the three LULC classes in areas not affected by the fire (outside the fire perimeter). Their pixel values were used as reference unburned training points for the RF model. For each year of calculation of the BRR, the respective ROIs BRR pixels of the same year were used as trainers; meanwhile, all the pixels falling within the burned area were reconstructed by interpolating them with the individual pre-fire image values (June 2018). This month was chosen since it was the closest to the fire event to also be present for all the observed following years (2019, 2020 and 2021), implicating that the conditions of the vegetation were more similar to those at the time of the fire event, as already considered by De Luca et al. (2022). Furthermore, the cloud-free seasonal period was optimal to avoid using those pixels reconstituted from cloud gaps (Section 2.3.2) and the rain noise on the SAR signal. Thus, the new unburned pixels for the BRR calculation have been predicted. The RF regressor parameters values were set using those already calculated in De Luca et al. (2021,2022) and optimized for the study area through an exhaustive grid search approach.

Vis temporal profiles
The temporal profiles of S1 and S2 spectral indices (Figures 2-4) show the mean reflectance value of the pixels retrieved from ROIs, separated by LULC class and fire severity category. The time series, divided between the pre-fire period (April 2017 -Jul 2018) and the post-fire period (second half of August 2018 -Jun 2021), is composed of images representing the monthly averages of the observed timeframe. Different patterns are distinguishable in each temporal profile, with a common breakpoint denoting the fire event in August 2018. An evident and expected difference between the two types of sensors is also found. Other noticeable differences in the LULC class, fire severity category, seasonal variability, and long-term trend are also observed. Large variability in the short-term post-fire behaviour (within the first year after the fire) is evident in both S1 and S2 indices but is significantly attenuated during the second and the third post-fire year in the S2 time series.

Optical Sentinel-2 profiles
As expected, for each LULC class, similar behaviour is observed between the different temporal profiles divided by fire severity category in the time frame preceding the fire event. Such profiles thus outline LULC-specific spectral signatures, an aspect much more evident in optical-based profiles than in SAR-based profiles (Figures 2-4).
After the fire event, different recovery patterns are recognizable according to the severity category in all the optical indices. Two common patterns are observed among all the profiles. A pronounced drop of the Vis values, with a magnitude related to the fire severity, in the period immediately following the fire event and maintained for about a year; a recovery of the pre-fire levels of spectral response (positive increasing trend) after the first year.
Observing the unburned category, the three vegetation classes had specific seasonal patterns during the observed period, with slight variations in the temporal dynamics between the different spectral indices in specific years. The highest index values occur in the same period for all the LULC classes: the autumn-winter period with peaks observed in December in almost all periods of all spectral indices, with some variation found in the two SWIR-based indices. The AuFor class showed less pronounced seasonal fluctuations and lower amplitude of the curve, opposite to the Euc profiles characterized by higher amplitude. Keeping the minimum values higher than the other classes, the AuFor secondary growth phase takes place in the earlier summer period between June and July, attributable to the growth of new leaves, with the only exception for the year 2019, in which there is no summer increase in the values of all the analysed spectral indices. This phase is anticipated by a minimum value occurring during the second half of May (NDVI, GNDVI, NDRE) or April (NBR, NDWI). The Pas/Shr class, generally characterised by lower VIs values compared to the forest classes, presents its single and regular period of annual growth that goes from September-October (minimum value), when the autumnal rains break the summer drought, through the whole winter, reaching the maximum peaks in December-January, and until they start to dry again at the end of spring.
The AuFor and Past/Shr post-fire profiles relative to the low severity category profiles present values very close to the reference unburned vegetation profiles in all the spectral indices. The intra-annual variability of phenology is clear in the time series even after the fire occurrence. The periods of growth and decrease are comparable between the two severity categories along almost the entire post-fire period. The impact of the low severity fire on the surface occupied by eucalyptus trees led to a lowering of the values of all the calculated optical indices. However, the magnitude of the recovery is the same as the corresponding unburned profile. Noticeable is the higher amplitude of the Euc low-severity profile compared to those falling into the other severity categories of the same LULC class. After a detachment the first year after the fire event, the moderate-high and high fire severity profiles showed similar values and followed the cyclical trend already observed for the unburnt category. However, the spectral response of the high fire severity category showed a steeper decrease immediately after the fire, more pronounced in SWIR-related VIs. Concerning the recovery degree observed in the temporal profiles, a complete recovery seems to be recorded for AuFor and Pas/Shr classes in all the VIs, except for the SWIR-based indices, highlighting a slower recovery of AuFor vegetation. Unlike what was observed for the other two classes, Euc profiles failed to restore the respective pre-fire Vis values within the analysed post-fire period.

SAR Sentinel-1 profiles
The comparison of SAR-based profiles denotes significant differences in the interactions between SAR indices and each LULC class, including those related to the ascending or descending orbit paths (Figure 3). The characterization of specific pre-fire seasonality and post-fire recovery patterns is complex in most profiles, mainly in the case of the singlepolarized VV backscatter signal (Figure 4). Indeed, the fire event is visually perceptible, observing the profiles. A distinction between pre-and post-fire patterns is observed for all the LULC and fire severity combinations except the VV-related profiles. Furthermore, copolarized profiles tend to slightly increase values in the post-fire period in all severity categories, but with an amplitude of the curve fluctuations that increase with fire severity. The profiles of dual-polarized SAR indices (Figure 3) are similar to those of the optical indices, exhibiting a better relationship between the fire severity category and the SAR signal ( Figure 2). The RFDI represents an exception, which shows an expected increase after the fire. The post-fire response of the SAR indices is more evident for high and moderate-high severity categories than for moderate or low severity and stronger for Euc and Pas/Shr. Indeed, the SAR response of vegetation affected by low severity fire does not present distinguishable traits from the unburned vegetation profile in most cases. Conversely, Euc low severity profiles differ from unburned and high severity profiles during the first post-fire year. The common characteristic of all SAR profiles for all LULC classes and fire severity categories is a slower recovery trend than optical profiles. In the case of Euc, even after three years, the indices' values do not reach those of pre-fire situations ( Figure  3). In the case of pre-fire or unburned vegetation, a general slight decrease is visible over the whole period when DPSVI and RVI are used (resulting in an increasing trend for RFDI), as observed in optical indices profiles. For the Euc class, the decreasing trend can be observable until November 2019 -February 2020, after which a slight increase is detectable. The same behaviour appears in ascending and descending VH backscatter profiles for all LULC classes and ascending VV backscatter for AuFor.

Climate variables and relationship with satellite temporal profiles
The Bagnouls-Gaussen Thermo-pluviometric diagram in Figure 5 shows the temporal profiles (April 2017 -January 2021) of the monthly rainfall and monthly mean temperatures, highlighting the Mediterranean climate of the study area (mild and rainy autumn, winter and early spring seasons, dry summers). Notable decreases in rainfall are observed in December 2018 (12.6 mm), March 2019 (14.8 mm) and February 2020 (3.96). Another exceptionality is given by the maximum peak recorded in March 2018 (266.7 mm); other maximum peaks were observed in October 2018 (122.5 mm), December 2019 (129.84 mm) and November 2020 (112.2 mm). The seasonal thermal variation is not high, and this is due to the oceanic influence that mitigates the average temperature during the year.
The Pearson's correlation coefficient between climate variables and temporal trends of both SAR and optical data was computed to investigate their linear relationship, considering the fire severity category, the LULC class and period considered (pre-or post-fire). The resulting Pearson's coefficient values are reported in the heatmaps of Figures 6-8. For optical datasets (Figure 6), a general positive correlation between rainfall and spectral indices can be observed for most LULC classes and fire severity categories, higher for Pas/Shr and the post-fire period. The NDRE index of non-forest vegetation was most Figure 6. Heatmaps reporting the Pearson's correlation coefficient performed between optical vegetation indices (NDVI, GNDVI, NDRE, NBR and NDWI) and climate variables monthly rainfall (Rain) and monthly mean temperature (Temperature). The correlation was carried out separately between the pre-and post-fire period and fire severity categories. The heat map colour palette gradient describes the Pearson's values range from lower (dar blue) to higher (yellow-light green).

Figure 7.
Heatmaps reporting the Pearson's correlation coefficient performed between SAR single polarizations (RVI_As, RVI_Ds, RFDI_As, RFDI_Ds, DPSVI_As, DPSVI_Ds) and climate variables monthly rainfall (Rain) and monthly mean temperature (Temperature). The correlation was carried out separately between the pre-and post-fire period and fire severity categories. The heat map colour palette gradient describes the Pearson's values range from lower (dar blue) to higher (yellow-light green). correlated with the rainfall in the pre-fire period, with values not lower than þ0.48. The degree of correlation for the post-fire period decreases with increasing fire severity, with the minimum value recorded for the NBR profile used AuFor class affected by high severity (þ0.087). The AuFor class mainly showed a lower correlation with rainfall than the other two LULC classes. Oppositely, most of the correlations with temperature were negative, with values that reached À0.89 (Pas/Shr-GNDVI-moderate severity-pre-fire). The AuFor class had some exceptions, presenting positive correlations for the SWIR-based indices (NBR and NDWI), with a higher positive correlation (þ0.38 to þ0.46) for low fire severity. Figure 8. Heatmaps reporting the Pearson's correlation coefficient performed between SAR dual-polarized vegetation indices (VH_As, VH_Ds, VV_As, VV_Ds) and climate variables monthly rainfall (Rain) and monthly mean temperature (Temperature). The correlation was carried out separately between the pre-and post-fire period and fire severity categories. The heat map colour palette gradient describes the Pearson's values range from lower (dar blue) to higher (yellow-light green).
Concerning SAR datasets (Figure 7), the RVI and DPSVI indices reached a similar correlation in their respective orbit path. For AuFor class, these indices had a higher correlation with the rainfall variable in the pre-fire condition, showing values between þ0.14 (RVI ascending) and þ0.67 (DPSVI ascending). This class reached a higher correlation in the two indices in post-fire conditions. For Ecu and Pas/Shr the highest positive and/or negative correlation was observed in pre-fire conditions, ranging from À0.47 (RVI descending) to þ0.5 (DPSVI ascending). During the post-fire period, both Ecu and Pas/ Shr achieved values not exceeding þ0.2 in any case. The RFDI showed values of opposite correlations to RVI for both orbit paths in Euc and Pas/Shr classes. Considering the AuFor class, the RFDI reached a similar and opposite to the DPSVI. A higher correlation was found between single-polarization profiles and the rain variable. Noticeable is the high level of r reached by VV_Ds in all three LULCs during the post-fire period (0.24 < r < 0.75). Except for some cases (e.g., VV_Asc for Pas/Shr; VH_Ds for Euc), the other single-polarized combinations followed this trend in both pre-fire and postfire periods.
A divergence was observed between dual-polarized indices and single-polarized backscatter and between ascending and descending orbit paths when correlated with the monthly mean temperature. After the fire event, the unburned and low severity categories correlated to the pre-fire period using dual-polarized VIs in descending mode. With ascending VIs, a slight divergence was observable for the low severity category. The other fire severity categories registered similar behaviour, except AuFor, which changed signs from negative (pre-fire) to positive (post-fire) or vice versa (RVI_Ds). Generally, pre-fire condition guaranteed highest correlation coefficient (e.g., r ¼ 0.81, RVI_As, Euc; 0.77, DPSVI_As, Euc; 0.61, DPSVI_Ds, Pas/Shr). Dissimilar behaviours were found using single-polarized backscatter profiles. The unburned vegetation did not maintain its correlation level with the climatic variable in all the cases. The VH_As was an exception, although this data presented the highest continuity between pre-fire and post-fire values in almost all the categories. The VV_Ds showed the highest negative correlation with temperature in the post-fire period, totally contrasting with the values expressed before the fire event.
Spatial distribution of recovery rate: the BRR Figure 9 shows the spatial distribution of the general recovery rate of vegetation for the three VIs NDVI, NBR and RVI_As. Areas in dark red (BRR < 0) denote slow/absent recovery, with values of the VIs lower than in the pre-fire conditions; conversely, where the BRR was more than one (light green areas), the recovery condition was higher than the pre-fire condition. Considering the different VIs, the three maps showed noticeable differences, with BRR NBR reporting recovery levels much lower than BRR NDVI . This disparity is kept during subsequent years, although the spatial pattern of re-greenness is quite similar, with a recovery distribution starting from the area around the unburned island (Monchique town) located on the west side of the study area toward the easter and north part/side of the study area. Other high-recovery areas are patchy, and spread throughout the 2021 map. Concerning SAR BRR, the recovery level has higher absolute values starting from 2019, and its spatial distribution is more homogeneous across the area. However, comparing the two highest recovery categories during the temporal progression from 2019 to 2021, the SAR shows a lower and more gradual regrowth rate than BRR NDVI . Some dark red areas are viewable in the year after the fire, recognizable by a sudden decrease in BRR values compared to the previous year, resulting from other fires' recurrence after August 2018. Although present on the map (Figure 9), these fires have been detected in advance and excluded from any calculation (ROIs BRR and pixel distribution). Figure 10 reports the proportion of the area occupied by each BRR category for each LULC class and year, assessed for the NBR, NDVI and RVI_As indices (BRR NDRE and BRR RVI_Ds are shown in Figure 1S of Supplementary Material since they displayed similar values of BRR NDVI and BRR RVI_As , respectively). Figures 9 and 10 show that the BRR NBR points out lower degrees of recovery for all the LULC classes. The substantial difference is denoted in the percentage values of pixel Figure 9. Burn recovery ratio (BRR) was separately calculated for NDVI (BRR NDVI ; first row), NBR (BRR NBR ; second row) and RVI_As (BRR RVI_As ; third row). The columns separate the three post-fire dates of observation: June 2019 (first column), 2020 (second column), and 2021 (third column). The colour palette represents a recovery gradient from very low (dark red) to very high (light green). distribution in the two lowest recovery classes (very low and low), for Euc (40.88% and 27.36%) and Pas/Shr (55.50% and 16.44%) classes on the first date after the fire (Jun 2019). For AuFor class, the pixels falling in the Very low category were only 4.15%, but Figure 10. Distribution of the number of pixels (%) falling in every burn recovery ratio (BRR) category, separately calculated for: vegetation index, NDVI (BRR NDVI ; first row), NBR (BRR NBR ; second row) and RVI_As (BRR RVI_As ; third row); date of observation, June 2019 (internal ring), June 2020 (medium ring), 2021 (external ring); LULC class, AuFor (first column), Euc (second column) and Pas/Shr (third column). The colour palette represents a recovery gradient from very low (dark red) to very high (light green). the 23.25% resulted in low recovery on the same date. The other optical and SAR VIs did not classify pixels into the very low BRR category, while the number of pixels labelled as low recovery was less than 1%. However, as expected, in all the cases, a positive recovery trend was detected; indeed, a shift in the distribution of pixels from lower to higher recovery categories when increasing years after the fire is observable. Concerning SARbased BRR RVI_As , although it appears that most of the pixels fell within higher recovery levels already starting from the first year after the fire, the relative increase in the number of pixels within the same BRR category for each observation year is significantly lower than that observed in the optical data. Looking at the categories of high and very high recovery, respectively, the AuFor went from 10.99% (2019) to 59.26% (2021) and from 2% (2019)

Optical time-series
This study reveals a strong relationship between vegetation recovery rate and fire severity category (Figure 2), as expected based on existing literature (Meng et al. 2018;Ryu et al. 2018;Morresi et al. 2019). In higher severity categories the maintenance of lower values of the profiles during most of the observation period is due to the double effect of the loss of photosynthetic vegetation and the presence of charcoal and ash (Pereira et al. 1999). Each LULC class showed a specific temporal pattern. Typically, tree species' predisposition to fire resistance and/or adaptation affects the recovery rate and patterns ( Catry et al. 2012;Filipe X Catry, Pausas, et al. 2013;Chompuchan and Lin 2017;Gouveia et al. 2010 ), especially for low and moderate fire severity. The seasonal stability of AuFor vegetation, characterizing the unburned vegetation, and the highest annual variability of Pas/ Shr, with the relative differences detected by NIR-or SWIR-based VIs, comply with what was observed by Soares et al. (2022) concerning the phenology adaptation strategies of Portuguese cork oak ecosystems.
Generally, a gradual spectral post-fire recovery is displayed in almost all cases at the end of the observing period, mainly when the species are characterized by a high regeneration capacity and vegetative growth (Semeraro et al. 2019). In the final monitoring year, the optical profiles no longer present dynamic patterns and differences between the various severity categories. Similar behaviour of recovery patterns of the photosynthetic activity and the ecosystem's physiological cycles after a few years from the fire event was reported in other studies (Chompuchan and Lin 2017) using optical profiles based on NIR and SWIR bands. As expected, the lower and moderate severity categories faster achieved advanced recovery stages in AuFor and Pas/Shr, near or equal to the initial unburned VIs values. This is supposed to be due to vegetation generally remaining totally or partially unburned, contributing positively to the optical reflectance. Meng et al. (2018) observed that a low level of fire-inducted damage does not induce such a significant reconstruction action by tree species to be detected by the VIs. This statement could justify the Euc low severity profile behaviour, characterized by a parallel pattern to the unburned profile; however, further and longer-time investigations are needed to explain these aspects.
Additionally,  demonstrated in their study that although Eucalyptus spp. has a general fire-resistant behaviour, it presents higher fire-susceptivity when placed in artificial plantations. It should be considered the result reported by Hausler et al. (2018), where it was observed that four years after the fire event, many eucalyptus surfaces have not yet fully recovered the levels of evapotranspiration, thus suggesting that it may take longer for a full recovery. These results highlight the higher efficiency of native species in restoring the ecological equilibrium. These observations are corroborated by the spatial distribution of recovery rate, represented by the BRR, as better illustrated in Section 4.5.
Regarding the interpretation of forest LULC classes profiles, a specific influence of the underlying vegetation should be considered (Meng et al. 2018). Where plant structural layers are vertically superimposed, gaps between the dominant foliage that leave the underlying vegetation uncovered make it very difficult to distinguish the unique spectral signatures of each vegetation layer by optical sensors. (Lacouture et al. 2020). This is commonly observable in moderately burnt areas, where the gaps in the foliage caused by the fire expose the underlying vegetation that has regrowth in the meantime (Polychronaki et al. 2014;Meng et al. 2018). This behaviour can be ascribed to the ability of optical sensors to detect only the external reflective surface of objects and to the saturation of the optical VIs at relatively low levels of LAI, with an increased rate that could be induced by the rapid colonization of herbaceous and shrubby vegetation. (Frazier et al. 2015;Zhao et al. 2016).
The gradual declension of the VIs profiles perceptible in the unburned forest classes (AuFor and Euc), opposite to the increasing trend associated with a frequent and growing fire disturbance regime, was already observed and described in other studies (Goetz et al. 2005;Ryu et al. 2018;Semeraro et al. 2019). The stasis in the physiological trend may be due to a combination of different factors. The reaching of the point of maximum vegetation increase in which the levels of competition between plants are lowered, and the rainfall decrease, lengthening the drought period ( Figure 5).
Although annual optical fluctuations seem more influenced by season-specific phenological dynamics than by rainfall's direct effect, it should be considered that plant phenology is closely related to the climatic variables, and this is reflected in all-optical VIs (Song and Woodcock 2003;Fern andez-Garc ıa et al. 2018;Poon and Kinoshita 2018;Ryu et al. 2018;Morresi et al. 2019). The correlation indices displayed in the heatmaps ( Figure  6) show that precipitation positively impacts physiological activity, especially for Pas/Shr and Euc vegetation. Significant drops in optical VIs occurred in these two LULCs classes between the end of June and August, corresponding to a sharp decrease in total rainfall during these months, highlighting that some vegetation types may be more sensitive to a meteorological parameter than another. The lower correlation with rainfall regimes of AuFor describes the low susceptibility to water stress that characterized Mediterranean corks, corroborating what was found by the recent study of Soares et al. (2022) and what was observed by Vidal-Macua (2017). The former study observed a high sensitivity of herbaceous species and shrubs belongings to the cork oak ecosystems to high temperatures (negatively) and precipitations (positively), supporting our findings (Figure 6). Concerning the temperatures, the greater susceptibility to high summer temperatures is confirmed, concomitant with a phase of vegetative stasis.
Focusing on the observation, the NIR-and RedEdge-based VIs showed faster recovery trends than SWIR-based VIs. This confirms what other scholars (Wang et al. 2005;Tanase et al. 2011;Frazier et al. 2015;Zhao et al. 2016;Ryu et al. 2018;Morresi et al. 2019) stated concerning, on the one hand, the greater sensitivity of the combination of the Red and NIR bands to photosynthetically active vegetation and the LAI, induced by the rapid regrowth of vegetation canopy; on the other hand, the complementary SWIR-based VIs sensitivity the immediate damage caused by the fire on the forest cover and to the gradual development of the canopy cover, which might take even a few years longer (1-4) than the NIR-based ones, as observed in other studies.

SAR time-series
The comparison of the different profiles shows how the interactions of the SAR signal, both individual polarized backscatter (Figure 4) and dual-polarimetric VIs (Figure 3), are not so easy to interpret as in the case of optical reflectance. Generally, our results show that the cross-polarized backscatter and dual-polarized index (RVI and DPSVI) values of the burned areas initiate an incremental trend phase, representing the beginning of the vegetation recovery phase, starting from the first year after the fire event in a way that is directly proportional to the fire severity category and as a function of the type of vegetation (LULC class). This was after an initial decrease immediately after the fire event, which was also proportional to the fire severity. However, the recovery dynamics are different between the type of VI used, and the characterization of specific phenological patterns is complex. This was expected given the greater sensitivity of the SAR signal to local environmental parameters. In particular, surface and soil moisture can heavily influence SAR recovery backscatter Tanase et al. 2011;Hachani et al. 2019;Zhou et al. 2019), while less influential seems to be the biochemical composition of the exposed soil (Minchella et al. 2009). Nonetheless, geomorphology and topography can directly or indirectly affect hydrological processes and, consequentially, vegetation recovery dynamics (Vidal-Macua et al. 2017). The drainage or retention capacity of the rock basement or the slope's orientation towards the sun's direct radiation determines local ground water balance and evapotranspiration Viana-Soto et al. 2017;Christopoulou et al. 2019).
The indirect effects that fires of this size and severity can cause surrounding unburned vegetation should also be considered. The micro-climatic conditions and the hydrological balance are altered by the destruction of the proximal coverings, whose mitigating and balancing properties are well known (shading, wind repair, humidity balance, rainfall interception, etc.). The artificial ecosystems (such as eucalyptus plantations) used for biomass production are more susceptible to these alterations, both for their ecological vulnerability and higher nutritional needs (Catry et al. 2015;H€ ausler et al. 2018). This could partly explain the decline in the unburned SAR curve of the Euc. Differentially, native species in a natural or semi-natural environment are predisposed to adaptive strategies to environmental adversities (Catry, Pausas, et al. 2013;Soares et al. 2022).
The results of the interactions between SAR polarization and surface dynamics can have different meanings depending on the observation biome Chen et al. 2018). For this reason, coupling climate variables profiles ( Figure 5) to spectral SAR response is helpful. The interception of rain by vegetation, especially by the forest cover, plays a fundamental role in the water balance of the ecosystem, resulting in a difference between the total rainfall recorded and that reached the ground, which would greatly influence the SAR signal (Frison et al. 2018;Ban et al. 2020). Although few intra-seasonal dynamic fluctuations were found in the SAR profiles, many of them (e.g., all the dualpolarization combinations: RVI, DPSVI, and RFDI ascending for AuFor and Pas/Shr; RVI, DPSVI and RFDI descending for AuFor class) showed evident inter-seasonal patterns driven by humidity conditions, in particular by rain ( Figure 5) when compared to the relative plot ( Figure 4). Noticeable are the effects of the exceptional amount of rain that fell in March 2018 and of that rainfall that fell in November-December 2017 on the response of the SAR signal, observable in all combinations. The peak of the SAR curve relative to June 2017, noticeable in most of the profiles, is difficult to interpret; this behaviour, among other things, is not found in the patterns following the fire event both for the unburned vegetation and that affected by low/moderate severities. A notable detail is a minimum value that is observed in many SAR profiles (VH, DPSVI, RFDI, RVI_As/Ds for AuFor and Euc; VV_Ds for all the LULCs; VV_As for AuFor and Euc) and almost all related severity categories, and all the relative categories of severity in the months of January-February 2020. Observing the plots of the climatic variables ( Figure 5), this can be correlated to an exceptional decrease in rainfall (3.11 mm in February 2020) that occurred in conjunction with exceptionally mild temperature levels (14.35 C in February 2020: þ1.75 C compared to February 2019þ 7.85 C compared to February 2018). However, the co-polarised signal seems to have been perceptibly driven by the rainfall in all the possible path-severity-LULC combinations, and more than the other profiles, with parallelism between curve fluctuations and a general increasing trend after the fire event. Ban et al. (2020) explained that the rain events after the fire and intercept bare soil, encountering little or no interception by the destroyed vegetation, influence the co-polarized signal inducing its increase.
Furthermore, the same authors observed that the VV backscatter values were generally higher than VH, an aspect also found in the profiles of the current study. Minchella et al. (2009) also observed increasing VV backscatter using C-band in Mediterranean forest ecosystems affected by the fire. The VV dynamic patterns of burned vegetation were driven by soil moisture fluctuations during the monitored five years after the fire event, especially at high fire severity. Periods at a lower soil dielectric constant, such as summer, resulted in lower VV backscatter values. Indeed, rapid re-growth of forest cover increases the similarity to herbaceous vegetation patterns (Minchella et al. 2009). These observations could explain the presence of high amplitudes (high fluctuations) during the different seasons of VV profiles, in some cases proportional to the fire severity category. However, given the relatively short observation period (three years), the attenuation effects given by the increasing vegetation recovery were not observed.

SAR polarization dependency
The abovementioned aspects concern the more general topic of differences between crossand co-polarization in interacting with the surfaces covered by forest. Immediately after a severe fire event and for the first post-fire year, oppositely to the perceptible increase already described for VV backscatter, a noticeable decrease in the VH backscatter is observable (Figure 4). This opposite behaviour of the two polarizations accords with other studies (Imperatore et al. 2017;Mari et al. 2017;. The burn of stems and large branches reduces the volumetric backscattering contribution of these structural vegetation components (i.e., scatterers), to which the cross-polarized signal is sensitive. Consequently, as microwave penetration through the damaged canopy increases, the effect of surface and double bounce dispersion on total backscatter is more prominent, induced by the greater proportion of exposed soil moisture components and the underlying soil (to which the co-polarized signal is sensitive) Tanase et al. 2011;Martins et al. 2016;Chen et al. 2018). In unburned forested areas, the high reflection pattern of VV co-polarized signal is instead associated with backscattering returned from vertical stand largest trees and trunks (Martins et al. 2016).
As expected, except for the VV_Ds profile for Euc class, no ability to distinguish between different severity categories was observed in all co-polarized combinations. When the severity of the fire that affected the forest vegetation was moderate, the decrease in VH is more gradual as the trees most affected are those of the intermediate vertical layer, which contribute less to the total AGB than the specimens of the dominant layer. Regarding the VH profile for the Pas/Shr class, it is more difficult to interpret its patterns, especially for the higher severity categories.
Due to their biases, a clear conclusion cannot be obtained when only using singlepolarization temporal backscattering, despite the higher overcome of cross-polarized already observed in other studies Martins et al. 2016). Agree with Chen et al. (2018). Instead, dual-polarized VIs improved the delineation of the proportionality between fire severity and time profiles, demonstrating their attitude to decrease the cross-and co-polarized biases.

SAR orbit path dependency
The separate analysis of the orbit paths made it possible to delineate the divergences arising from the sensor geometry on the temporal profiles. This was due to the look direction of the SAR, which causes the signal to interact geometrically differently with the local structural orientation of the regrowing forest (M. Tanase et al. 2011). A fundamental role is played by the effects of the terrain, which determines the local angle of incidence of the microwave beam (Gimeno and San-Miguel-Ayanz 2004;Kurum 2015). The presence of reliefs can cause radar shadow effects, for which some areas facing the side opposite the radar beam are not detected by one of the two paths (Richards 2009a), as well as directly influencing the local soil water balance Viana-Soto et al. 2017;Christopoulou et al. 2019). This confirms the need to use both orbit paths to construct complete information about what is observed. Some authors have highlighted the role of both orbits in increasing the detection capabilities of burned areas (Donezar et al. 2019; or the LULC classification (Sayedain et al. 2020). A further observation is that the differences between orbit paths are noticeable when the single polarizations are compared. Integrating these two in the dual-polarized indexes seems to attenuate the look direction effects.

Spatial distribution of vegetation recovery rate, the BRR
The spatial distribution of the vegetation recovery categories, represented by the optical BRR (Figures 9 and 10), complies with what emerges from the analysis of the temporal profiles. The BRR results ( Figure 10) demonstrated that the surface occupied by native forest reaching the highest recovery level was more than Euc vegetation at the end of the observed period (2021). Even observing the map of Figure 9, it is clear that the west-east development track, spatially and temporally traced by the vegetation regrowth, follows the disposition of the areas occupied by AuFor vegetation when compared to the LULC map presented in De Luca et al. (2022). When optical VIs based on NIR (BRR NDVI ) and RedEdge (BRR NDRE, in supplementary material) are used, most pixels fell into the higher recovery categories since the first year after the fire event than when using the index NBR. Furthermore, the SWIR-based index includes most of the pixels within the lowest recovery classes (Very low and Low) in the first two post-fire years. The other optical indices did not detect pixels in the very low category, while less than 1% of the pixels fell in the low category. The season of occurrence of the fire also influences the characteristics with which it affects the vegetable fuel. A higher humidity of the fuel could increase the variability of the fire and, equivalently, of the VIs responses (Lacouture et al. 2020), especially of SWIR-based ones. Although, in the present study, having the fire occurred during the dry season (August 2018), it was characterized by homogeneous severity and intensity over large and contiguous areas (De Luca, Silva, Oom, et al. 2021).
The relative increase of BRR RVI recovery in the subsequent year was noticeable compared to the respective optical indices. The BRR relative to the RVI index showed that most of the pixels already fell into higher recovery categories in the first year after the fire. However, the spatial distributional track of the recovery adheres to that of optical BRRs. Anyway, if the relative change over time is analyzed and compared to the BRR NDVI , a more gradual recovery is detectable for the BRR RVI_As , as already found in the temporal profile. However, the more complex interpretation of the SAR makes it more challenging to characterise specific spatial patterns. Factors such as the different influence exerted by the remaining charcoal on the fire-affected surfaces, the high sensitivity to humidity, and the wavelength used suggest that a further examination is needed concerning the efficiency of using the SAR-based BRR.

Final considerations
Considering the findings of the present study, as well as the consulted literature, we found that the optical S2 and SAR S1 VIs appeared to have a higher aptitude for monitoring vegetation regrowth, which can be used effectively as complementary information to assess and monitor the short-term response of ecosystems to fire. Some limitations persist, such as the predisposition of optical data to reach saturation, partially compensated by the weak SAR C-band penetration capacity through the forest canopy. Further literature examination and considerations concerning the SAR wavelength dependency are reported in the supplementary material. Integrating multi-frequency SAR (L-band and/or P-band), VHR data, or LiDAR information can optimize the capture of the spectral and structural properties of the regrowth vegetation and quantify the understory seedling recovery rate more precisely (Meng et al. 2018).
Besides those already argued, many factors influence the SAR results or help their interpretation. The behaviour of the SAR depends very much on the environmental conditions, and further investigations will have to be carried out. In this study, two climatic variables are taken into consideration. However, we would like to suggest other parameters that should be investigated and compared, including pedological and lithological characteristics, soil moisture, topography, and territory geomorphology. Besides the heterogeneity of LULCs analyzed (especially AuFor and Pas/Shr), the topographical variability of the Mediterranean ecosystem strongly influences the microwave signal Tanase, Santoro, Wegm€ uller, et al. 2010;Tanase et al. 2011), as highlighted by the differences resulted from the dependent analysis of the orbital path. The use of averaged ROIs value for both the optical and SAR dataset, added to the speckle noise filter applied to SAR, attenuated the variability of each curve (Minchella et al. 2009;Tanase et al. 2011;Frazier et al. 2015;Martins et al. 2016). Several studies investigated how topographical and geomorphological variables positively or negatively influence local vegetation recovery dynamics and population post-fire biodiversity and how they are affected in turn Diakakis et al. 2017;Vidal-Macua et al. 2017;Viana-Soto et al. 2017;Christopoulou et al. 2019). Slope and aspect, for example, determine solar wave incidence radiation, which largely influences both water/moisture balance, evapotranspiration, and photosynthetic activity (R€ oder et al. 2008). On the other hand, slope and lithology affect erodibility, water drainage/accumulation, as well as root penetration capacity (Christopoulou et al. 2019). More considerations based on current literature are reported in the supplementary material.

Conclusions
The results reported in this study show that Mediterranean ecosystems respond rapidly to disturbances, initiating effective restoration processes. However, earlier recovery of unburned values is attributable to a premature saturation affecting both NIR-based indices and SAR C-band wavelength. Regrowth trends are observable from the first months after the fire event reaching an apparent almost complete recovery occurring by the three years of analysis in optical and SAR VIs. The vegetation type influences the time and the magnitude of recovery temporal activity in terms of spectral response, with a noticeable difference between native and non-native forest vegetation. However, the proportionality between recovery pattern and fire severity categories was kept. The higher degree of recovery occurred for the autochthonous forest class (AuFor), followed by Pas/Shr, reflected the higher adaptability of these ecosystems to stress regimes and efficiency in restoring the ecological equilibrium, directly proportional to their more differentiated structure and biodiversity compared to the area occupied by eucalyptus plantations. These observations were corroborated by the BRR, representing the spatial distribution of recovery rate, whose outcomes demonstrated that the surface occupied by highly and very highly recovered AuFor vegetation was more than that occupied by Euc vegetation during the three years of observation; even more evident when optical indices were used (BRR NBR , BRR NDVI and BRR NDRE ). In order to account for the natural phenological temporal effects, the BRR calculation was optimized using the RF machine learning regressor with which the hypothetical unburned conditions of the fire-affected area for each time of estimation of the index were predicted and reconstructed. However, although optical BRR has been demonstrated to be effective, the SAR-based needs further studies for its interpretation.
A high proportionality was found between the fire severity and the recovery profiles during the first two post-fire years, following an equally proportionated sharp decrease in optical and SAR values. Some exceptions persisted, such as for the SAR co-polarized (VV) profiles, due to their high moisture and exposed soil dependency. In fact, from the comparison and correlation of climate variables and temporal profiles, it has been clear that the precipitation events directly affected the C-band SAR profiles, especially in the categories of higher severity, where the interception of the vegetation cover guarantees the precipitation attenuation had been lacking. This denotes the menace of fire events toward the hydrological balance of the soil.
Considering these aspects, the management of SAR images required greater attention due to the higher presence of outliers and speckle noise caused by the intrinsic nature of the data. Using long time-series, the biases are more evident; however, several pre-processing processes were used to optimize the workflow and mitigate these issues, such as the speckle filter and the time average. Optical-related biases were also finely addressed, such as eliminating the clouds through the new product provided by Sentinel Team (s2cloudless) and subsequent linear interpolation to reconstruct the fill gaps.
Generally, this study demonstrated that the combined use of different sensors is essential to correctly delineate the dynamics that occur in Mediterranean fire-frequented habitats, compensating for the limitations of the single sensor, especially when a small temporal scale is needed. Moreover, the integration of free and open-source analysis software with equally free-available high temporal and spatial resolution data enables the accessibility by a wider audience involved in the forest risks monitoring framework.
However, the dynamics of post-fire Mediterranean vegetation must be further examined in long-term monitoring protocols in this and other study areas to assess the complete response even to delayed effects. Moreover, additional indicators and sensors may be necessary to determine which combination of temporal patterns best reflects the real post-fire dynamics in the Mediterranean ecosystems and their chemical, physiological and structural features. Focusing on SAR data, medium-long term monitoring may require the integration of multifrequency techniques with longer wavelengths (L-, P-band), able to penetrate further into the regenerated canopy, thus enabling in better understanding of forest recovery processes. The imminent availability of new data types (e.g., ESA BIOMAASS mission, Le Toan et al. 2011) and the development of cloud platforms (e.g., multi-mission algorithm and analysis platform, MAAP 2022) will optimize the performance and the already high inter-compatibility of these resources.
Further investigation could also involve the machine learning regression models to predict temporal and spatial recovery patterns, basing the regression on the values of the recovery metrics calculated. In this regard, we believe our results will efficiently provide helpful information.

Disclosure statement
No potential conflict of interest was reported by the authors.

Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request. These data were derived from the following resources available in the public domain, ESA Sentinel Homepage, https://sentinel.esa.int/web/sentinel/home.