Use of Sentinel-2 Satellite Data for Windthrows Monitoring and Delimiting: The Case of “Vaia” Storm in Friuli Venezia Giulia Region (North-Eastern Italy)

: On the 29th of October 2018, a storm named “Vaia” hit North-Eastern Italy, causing the loss of 8 million m 3 of standing trees and creating serious damage to the forested areas, with many economic and ecological implications. This event brought up the necessity of a standard procedure for windthrow detection and monitoring based on satellite data as an alternative to foresters’ ﬁeldwork. The proposed methodology was applied in Carnic Alps (Friuli Venezia Giulia, NE Italy) in natural stands dominated by Picea abies and Abies alba . We used images from the Sentinel-2 mission: 1) to test vegetation indices performance in monitoring the vegetation dynamics in the short period after the storm, and 2) to create a windthrow map for the whole Friuli Venezia Giulia region. Results showed that windthrows in forests have a signiﬁcant inﬂuence on visible and short-wave infrared (SWIR) spectral bands of Sentinel-2, both in the short and the long-term timeframes. NDWI8A and NDWI were the best indices for windthrow detection (R 2 = 0.80 and 0.77, respectively) and NDVI, PSRI, SAVI and GNDVI had an overall good performance in spotting wind-damaged areas (R 2 = 0.60–0.76). Moreover, these indices allowed to monitor post-Vaia forest die-off and showed a dynamic recovery process in cleaned sites. The NDWI8A index, employed in the vegetation index differencing (VID) change detection technique, delimited damaged areas comparable to the estimations provided by Regional Forest System (2545 ha and 3183 ha, respectively). Damaged forests detected by NDWI8A VID ranged from 500 m to 1500 m a.s.l., mainly covering steep slopes in the south and east aspects (42% and 25%, respectively). Our results suggested that the NDWI8A VID method may be a cost-effective and accurate way to produce windthrow maps, which could limit the risks associated with ﬁeldwork and may provide a valuable tool to plan tree removal interventions in a more efﬁcient way.


Introduction
The recent growth of satellite technology in the field of Earth observation allows the measurement of electromagnetic radiation with a higher spatial, temporal and radiometric resolution, being potentially useful to detect environmental disturbances caused by natural hazards [1,2]. This feature is particularly relevant considering future scenarios where extreme natural events and their impacts are expected to become more frequent and intense [3]. As a term of reference, more than half of the total recorded damages to European forested areas has been caused by extraordinary events of extremely intense winds as a consequence of the ongoing climate change, whose associated damage levels are steeply increasing [4].
In the future, there are several factors that could increase forest vulnerability to storms: among these, increased intensity of heavy storms and winds, changes to storm tracks, higher air temperature and precipitations during winter and, finally, an increased pathogen activity triggered by higher temperatures [4]. During the last centuries, moderate biotic and abiotic disturbances played a key role in European forest ecosystem dynamics, stimulating regeneration and ecological succession processes and having positive effects on biological diversity [5]. However, high-intensity disturbances can be catastrophic and have complex repercussions on the biosphere, such as diminishing ecosystem functions, reducing biodiversity and altering land-atmosphere feedbacks through changes in forest structure and its physiologic processes [6]. The susceptibility of individual trees and tree stands to wind damage, and its mechanism, is determined by the properties of the wind climate (wind speed, gust speed, duration), forest structure (e.g., fragmentation and age), tree/stand characteristics (e.g., species, height, diameter at breast height, crown and rooting characteristics and stand density) and site conditions (e.g., soil type and topography) [7][8][9]. Furthermore, several plant species prone to wind damages (i.e., endowed of shallow root systems, low height-to-diameter ratio and dense heavy crowns) were historically widely planted in Europe, thus increasing forest storm susceptibility, especially where heavy rain saturates soils reducing anchorage [4]. Forest resilience and resistance to windstorms have been further reduced by silvicultural practises which, by maximising production and minimize management, created dense monospecific and mono-stratified populations, even in age and composition [10].
From Saturday the 27th to Tuesday the 30th of October 2018, North Italy was subjected to the effects of a deep meteorological depression called "Vaia". Strong linear geostrophic winds were triggered by about 20 hPa pressure difference between the Western and Eastern Alps. Especially on the 29th, Vaia caused strong scirocco (on Eastern Alps) and libeccio (on Tirrenic coasts) winds that, combined with local winds and turbulences on the mountain valleys, generated gusts up to 200 km/h. The morphology and orography of the Alps amplified the phenomenon, channelling the wind in the valleys, thus increasing its speed [11], and compressing downslope currents. Eight million cubic meters of standing trees, mainly represented by Picea abies (L.) H. Karst., and Abies alba Mill. species were windthrown, with relevant implications for the environment, the economy, and the society [12]. In many Italian northern Regions, whole landscapes were catastrophically damaged together with the underpinned biodiversity and important forest-related ecosystem services were suddenly lost, including protection against landslides, avalanches and floods, and cultural and esthetical values. In most cases, the perimeters of impacted areas were manually registered on the field or derived from helicopter inspections [12]. Friuli Venezia Giulia, the most northeastern Region of Italy, manually digitalized single perimeters of all the areas subjected to the storm damage. At the national level, the most affected areas cover a total of 42525 ha, accounting for 8.5 million m 3 of timber. The damages in the Friuli Venezia Giulia Region were estimated to be around 3340 ha, accounting for 723 thousand m 3 of timber.
Nowadays, there is a strong need to better understand and identify the mechanisms of extreme events like Vaia, to analyse their impacts, and plan future natural risk management [2]. In this context, the emphasis is placed on satellite-derived sensing because these data provide a low-cost synoptic view, which is fundamental for risk monitoring and early hazard detection. Many projects and organizations rely on satellite data such as the European Copernicus EMS project (Emergency Management Service), which is devoted to early warning improvement, early impact assessment, information sharing and early emergency coordination [13]. In particular, the Sentinel-2 satellite mission provides free access to multispectral images. Sentinel-2 products are affected by vegetation components' capacity to absorb or reflect sunlight, thus providing important information on forest status after storm events, assessing the changes in many biophysical parameters. In particular, Vegetative Indices (VIs) enable the acquisition of ecological information from satellite data, through the algebraic combination of imagery bands (radiance values at different wavelengths intervals), enhancing the spectral information thus increasing the separability of the classes of objects of interest [14,15].
Regarding the damage estimation caused by the Vaia event, the Copernicus Emergency Mapping system was activated along with the national estimation operations, but it had been able to map just 4000 ha, less than 10% of the total damaged areas [16]. According to Chirici et al. [12], Sentinel-2 images have suffered shadows and snow presence, as well as cloud cover that occurred in the weeks after the event. More than 50% of the concerned surface still did not have a cloud-free image at the end of January 2019, 3 months after the Vaia storm.
Despite not providing immediate data for early damage detection, Sentinel-2 could provide precious information for post-Vaia ecological dynamics comprehension and monitoring, and can increase the accuracy of available maps in an efficient tree removal planning perspective. In this study, by using different image processing and analysis software, we aimed to: (1) determine the best band transformations (vegetation indices) to remotely separate damaged and un-damaged forest stands after a storm and to monitor the vegetation dynamics in the short period after the storm. (2) employing multispectral data in an image differencing change detection technique for windthrow detection in the whole Friuli Venezia Giulia Rregion.

Study Sites
This study was performed in Carnic Alps, Friuli Venezia Giulia (Italy). The study area is located in Sutrio municipality, where large areal windthrow was recorded. The study site covers the south-eastern side of Zoncolan Mount (46 • 30 08"N, 12 • 55 34"E), and it is located in an area ranging in altitude from 750 m a.s.l. to 1250 m a.s.l. An additional site with the same S-E aspect has been chosen in Paularo municipality at 1250 m a.s.l. (46 • 31 51"N, 13 • 07 07"E). The climate is cold-temperate with mild summers and cold, but not particularly harsh, winters. Mean annual rainfall exceeds 2000 mm, mostly occurring in spring and autumn. In the study areas, winds are influenced by land morphology: valleys orientations, depth, width, aspect and solar radiation [17]. In this context, each mountain site differs in wind circulation characteristics even from geographically close areas. The synoptic winds triggered by storm Vaia reached 139 km/h gusts on Zoncolan summit, with mainly south-east and south-west mean direction [18]. The vegetation here was composed of coniferous forests dominated by Picea abies and Abies alba. Fagus sylvatica (L.) was also present but mixed with the previously mentioned main forest types.
The Vaia storm caused a heterogeneous damage pattern in Italian Eastern Alps temperate forests, resulting in scattered damaged areas with single and small groups of fallen trees or entire large areas with a uniform complete tree loss. Four sites, characterized by a pre-storm homogeneous vegetation cover of Picea abies and Abies alba have been selected ( Figure 1): two were rather undamaged after Vaia, with an observed tree loss < 5%, and were referred to as control sites (C1 in Paularo, C2 in Sutrio), whereas L1, L2, both located in Sutrio municipality, had a uniform complete loss and were referred as windthrown sites. Each site covers a surface of approximately 1 ha, with a slope higher than 40% and south or east aspect. Slope, elevation and aspect maps were extracted from an 90 m resolution digital elevation model (DEM) processed using the "r.slope.aspect" function in GRASS 7.2 [19]. In "r.slope.aspect", the aspect output raster map indicates the direction that slopes are facing counterclockwise from East: 90 • is North, 180 • is West, 270 • is South, 360 • is East. We then reclassified slope and aspect maps in 4 percentage classes as follows: 0-5%, 6-20%, 21-40% and > 40% for slope, N (45 • -135 • ), E (0 • -45 • and ≥ 315 • ), S (225 • -315 • ) and W (135 • -225 • ) for aspect. The FVG 2017 Habitat Map [20] was used to select forest areas dominated by Picea abies and Abies alba. To avoid areas including deciduous species (i.e., Fagus sylvatica), pre-storm Sentinel-2 images acquired during the autumn season were used to visually check for the eventual presence of spots with colours typical of senescent leaves of deciduous species.  Table 1.

Remote Sensing Data
Sentinel-2 Level 2A (bottom of atmosphere) multispectral images with tile 33TUM were downloaded from Copernicus Open Access Hub [21] for the pre and post-Vaia dates showed in Table 1. Sentinel-2 images availability during late autumn and winter seasons was strongly reduced because of cloud cover, shadows, and snow presence. Accordingly, for the months immediately after the storm (from November to March), it was possible to retrieve only one image in November (15 November 2018) and one in March (5 March 2019).
In this study, 10 Sentinel-2 bands extend from the visible and near-infrared (VNIR) to the short-wave infrared (SWIR) were analysed in GRASS 7.2 [19], featuring: • four bands at 10 m: B02, blue (490 nm); B03, green (560 nm); B04, red (665 nm) and B08, near-infrared (NIR) (842 nm). , PSRI, SLAVI) were calculated using "r.mapcalc" function. Geographic boundaries for the sampling sites were set with the "g.region" and "r.mask" functions. The spatial resolution of the products is 10 m for GNDVI and NDVI. For the remaining indices, spatial resolution was set at 10 m by resampling lower resolution band pixels with the Nearest Neighbour algorithm with "r.resample" function in order to lose as little resolution as possible, considered the reduced extension of sampled sites. These indices take advantage of SWIR, NIR and RED, which have been proven highly effective in the detection of storm damages [22]. The indices used in this study are reported in Table 2. Table 2. Vegetation indices formula and description.

Reference Description
Green Normalized Difference Vegetation Index GNDVI = B08−B03 B08+B03 [23] GNDVI was suggested to avoid the NDVI saturation effect in dense canopy structures in order to correctly infer chlorophyll canopy content [24]. It uses the green band (B3, 0.60 µm) as a replacement for the red one. The values range from −1 to 1, with growing values as green biomass increases.
IRECI is one of the best indices for the estimation of leaf area index (LAI), namely the one-sided green leaf area per unit ground surface area. It has been also used to estimate canopy and leaves chlorophyll content. IRECI makes use of all Red-Edge (RE) Sentinel bands (B5, B6 and B7, respectively~0.7 µm, 0.74 µm and 0.78 µm, the latter corresponding to maximum reflectance) and the red band (B4,~0.65 µm, minimum reflectance) to characterize the RE slope. By the use of the B5/B6 ratio, IRECI does not put heavy emphasis on the red, which helps to avoid saturation at high chlorophyll concentrations, while utilizing the strong contrast of the B7 and B4 difference which is sensitive to LAI.
NDREDI is a very recent index evolved in a study for windthrow detection in Germany from PSRI. Both NDREDI and PSRI, by applying the RE and blue spectral bands, proved to be suitable measures for the detection of windthrows.
NDVI ranges between −1 and 1. Naked soils have values around 0, negative values are retrieved from water surfaces while positive values are found in vegetated areas, the higher values the greater is the biomass. It has been used for estimation of Leaf Area Index, biomass, chlorophyll concentration in leaves and plant productivity [28]. Table 2. Cont.

Index Reference Description
Normalized Difference Water Index 11 NDWI = B08−B11 B08+B11 [29] NDWI combines NIR (~0.85 µm) and SWIR (~1.6 µm) bands for water stress monitoring in plants. SWIR reflectance reflects changes in both the vegetation water content and the spongy mesophyll structure in vegetation canopies. NIR reflectance is affected by leaf internal structure and leaf dry matter content but not by water content. The combination of NIR and SWIR removes variations induced by leaf internal structure and leaf dry matter content, improving the accuracy in retrieving the vegetation water content. Due to the high sensitivity of NDWI and NDWI8A to changes in water content of leaves, this index is efficient in detecting damage, monitoring drought periods and estimating yield depletion [30]. Index values vary from −1 to 1, they are negative when the vegetation is dry or on bare soils, while they are positive when the vegetation is green and rich in water [31,32].
PSRI was introduced to estimate the onset, the stage, relative rates and kinetics of the senescence process. An increase of reflectance between 0.55 and 0.74 µm (sensed by B4 and B6) accompanied senescence-induced degradation of chlorophyll, whereas in the range 0.4-0.5 µm (B2) it remains low due to retention of carotenoids. The index is sensitive to the carotenoids and chlorophyll ratio and it has been used as a quantitative measure of leaf senescence.
Specific Leaf Area Vegetation Index SLAVI = B08 B04+B12 [34] SLAVI was developed to estimate the specific leaf area (SLA), namely the one side area of the leaf divided by its dry weight. It is based on red and near-infrared reflections, which have a high correlation with net photosynthesis, net primary production and leaf area index [35], and on SWIR reflection, which guarantees leaves water content estimation. SLAVI is negatively correlated with cellulose and lignin content. It is also sensitive to the volume and height of vegetation [36].

Statistical Analyses
Pixel values for each multispectral index were plotted for each study site (L1, L2, C1 and C2) and each date through a boxplot using "ggplot2" R package [37]. In addition, each Sentinel-2 band (B02, B03, B04, B05, B06, B07, B08, B8A, B11 and B12) mean value and standard deviation for wind-damaged and control sites were plotted for each date to create spectral signatures. Generalize Least Square (GLS) models were used to test whether remote sensing indices were different between windthrown and control sites per each date, using gls function in "nlme" package [38]. Specifically, a GLS model was calculated for each remote sensing index (response variable) as a function of the study site type (fixed effect, 2 levels: control vs. windthrown). Temporal dependence among observations was incorporated through an autocorrelation structure of order 1. Model performance was evaluated computing the variance explained by the fixed effects using the function r2beta [39], with Nakagawa and Schielzeth method, available in R package "r2glmm" [40]. Differences between areas for the same date were assessed performing a GLS model for each index (response variable) as a function of area, date and their first order interaction using gls function, as described above. For significant models, posthoc pairwise comparisons were evaluated by "emmeans" R package [41]. Differences in band values for wind-damaged and control sites were assessed performing a GLS model for each band (response variable) as a function of site type, date and their first order interaction with gls function as described above. For significant models, posthoc pairwise comparisons of fixed effects interaction were evaluated as detailed above. Statistical analyses were performed using R 3.6.1 [42].

Change Detection
One of the most widely used mathematical operators to highlight image differences and detect changes is the image differencing method [43]. The method simply subtracts preimage from post-image pixel values or vice versa. Despite the simplicity, image differencing methods have been widely used in previous studies devoted to wind-damage detection with considerable success [44][45][46]. In this study, the best VI identified by GLS models analysis (NDWI8A, see Results section) was employed in a (VID) technique in order to detect windthrown areas in whole Friuli Venezia Giulia Region Alpine forests. A new image called NDW I8A di f f (20 m spatial resolution) was simply derived from the difference of NDWI8A raster values before Vaia and after Vaia: Summer and early autumn are recommended as the most suitable times for forest evaluation [22]. The 21st of April 2018 and the 21st of April 2019 were proven to be the best images for NDWI8A difference calculation as they were not affected by seasonal variations. Images acquired in different periods strongly differ by shadows presence and illumination on the ground. Moreover, the comparable period acquisition guarantees low phenological differences. The choice of spring season avoided the inconveniences related to the yearly variable water availability for forests at high altitudes. It has been observed that, at altitudes where snow is accumulated during winter (above approximatively 1500 m a.s.l.), the autumn vegetation water content, detected through the NDWI index, strongly depends on previous winter snow accumulation. However, albeit in a limited manner, spring images were not exempt from issues: most areas above 1500 m a.s.l. were still covered with snow in April, thus not included in the study. NDWI8A diff values were masked using the Scene Classification Layers (SCL), retrieved in Sentinel-2 products, and the 2017 Regional Habitat Map [20]. In particular, NDWI8Apre was filtered for the SCL vegetation category while NDWI8Apost for both vegetation and not vegetated categories in order to include forest areas that turned not-vegetated because of the storm. NDWI8A diff was then calculated from filtered images and further masked using the areas indicated as wooded in the 2017 Regional Habitat Map. Finally, NDWI8Adiff was filtered for values > 0.15, identified as the best threshold for windthrows detection, and converted from a raster product into a vectorial polygon. Finally, the obtained windthrown vectorial products were filtered for surfaces bigger than 800 m 2 (2 pixels) in order to reduce the "noise" associated with scattered and spatially limited outputs.

Descriptive Statistics on Regional Windthrows
Windthrows detected in Friuli Venezia Giulia Alpine forests with NDWI8A vegetation index differencing method were compared to reference Regional Forest Service (RFS) data. Descriptive statistics on the coverage of windthrown areas identified by the VID were extracted according to the elevation (4 classes: < 500 m, 500-1000 m, 1000-1500 m, >1500 m a.s.l.), aspect (4 classes: North, East, South, West), slope (4 classes: 0-5%, 6-20%, 21-40%, >40%) and extension (4 classes: 0.8-2 ha, 2-10 ha, 10-25 ha, > 25 ha) and compared to the ones extracted from the RFS windthrows perimeters. Moreover, GIS-provided geospatial data statistics were also used to assess the degree of similarity between the reference RFS data and the novel windthrows map obtained with the NDWI8A VID method in RFS perimeters. Statistics on damaged forests by municipality for both NDWI8A VID and RFS were extracted and the % error of NDWI8A VID method with respect to RFS estimations was calculated as:

Map Validation
A reference map of Vaia windstorm impacts was generated based on field survey data provided by the Regional Forest System (RFS). The fieldwork took place in all Friuli Venezia Giulia Region after Vaia storm. Although the accuracy of this data is limited by surveyors' possibility to access and spot damaged areas, the RFS reference map provides a comprehensive collection of available data of Vaia-caused disturbance in forests across a large portion of the Regional territory, thus providing a useful tool for comparison with remote sensing data. For both RFS and NDWI8A VID maps, validation was obtained by field assessment along the valley of Tagliamento river, where, according to the reference material, most damages occurred. Field survey observations were substantiated by georeferenced pictures, which offered a descriptive insight into Vaia disturbances impacts, also providing a qualitative check of windthrow maps.

Temporal Profiles of Vegetation Indices
Boxplot of temporal profiles of the NDWI8A and PSRI indices both for windthrown and control study sites showed that indices values after Vaia in L1 and L2 windthrown sites are significantly lower for NDWI8A and higher for PSRI than in C1 and C2 control sites, differently, being almost comparable before Vaia (Figures 2 and 3). The remaining indices temporal profiles are shown in Supplementary Materials (Figures S1-S6). GLS model outputs for each index as a function of area, date and their interaction are shown in Supplementary Materials (Table S1). For each of the eight indices considered in this study (except IRECI and NDREDI), the values of the indices remained quite constant in control sites during all the period of investigation, taking into account seasonal variations. In NDVI, NDWI, NDWI8A, SLAVI, GNDVI, IRECI indices, L1 and L2 sites pixel values were comparable to control sites before the Vaia event, but they showed significantly lower values after it. Notably, the PSRI index showed an inverse pattern compared to the other indices. The index NDREDI showed differences between controls and windthrown sites on the first date after Vaia (2018.11), having comparable values in windthrown and control sites in the remaining dates ( Figure S6). An increase in control sites' variability was detected after Vaia for all indices, less evidently in NDWI8A compared to other indices.  Table 1. Letters denote statistical significance between control and wind-damaged sites index values (p-value < 0.05). The dash line separates pre and post-Vaia dates. Table 3 showed GLS model outputs; NDWI8A was the best index for windthrows detection accounting for 80% of variance explained (see R 2 in Table 3). NDWI, NDVI, PSRI, SAVI and GNDVI had an overall good performance in spotting wind damages (spanning from 0.60 to 0.77 of R 2 ). For IRECI the effect of the treatment explained 32% of data variability, while for NDREDI just 3%. In contrast, PSRI showed a significant decrease after the Vaia storm compared to the previous period.  Table 1. Letters denote statistical significance between control and wind-damaged sites index values (p-value < 0.05). The dash line separates pre and post-Vaia dates. Table 3. Generalized Least Square (GLS) models coefficients (intercept and treatment), standard error (SE), p value and variance explained (R 2 ) for NDWI8A, NDWI, NDVI, PSRI, SLAVI, GNDVI, IRECI and NDREDI indices (response variables) in L1, L2, C1, C2, as a function of site type (fixed effect; control vs. windthrown), * for p ≤ 0.05, ** for p ≤ 0.01 and *** for p ≤ 0.001. Indices are ordered according to their windthrows detecting performance (variance explained).

Spectral Signatures
Controls and Wind-damaged sites spectral signatures for each date after the Vaia event are shown in Figure 4. GLS models' outputs for each band are shown in Supplementary Materials (Table S2)   In November 2018, control and wind-damaged signatures followed a very similar pattern. Indeed, both showed two absorption peaks on blue and red regions (B02 and B04, respectively), an increasing reflectance from red to near-infrared (NIR) (from B04 to B8A) and lower reflectance values in the short-wave infrared (SWIR) region (B11 and B12). Control and wind-damaged signatures overlapped from red-edge to NIR (from B06 to B8A), but wind-damaged sites had significantly higher reflectance in the remaining bands. Control sites mean reflectance in the visible (form B02 to B04) spanned from 0.5% to 0.8%, in contrasts, damaged sites reflectance ranged between 1.3% to 3.5%. The reflectance of the first band of red-edge (B05) was higher in wind-damaged sites (5.8%), rather than in controls (3.8%). Differently from wind-damaged sites, controls had an evident absorption in the SWIR region. Accordingly, SWIR1 and SWIR2 reflectance was 11.2% and 6.8% in wind-damaged sites and 5.4% and 2.8% in controls.
In March and April windthrown sites spectral signatures, we observed increasing reflectance values in the B04 red band and in B11 and B12 SWIR ones. In March B04, B11 and B12 reflectance in windthrown sites (4.6%, 12.4% and 7.5%) was almost 3 times higher than controls reflectance (1.7%, 4.6%, 2.3%). In April, B04 and B12 reflectance in wind-damaged sites was 4 times that of the control. In April 2019, the highest reflectance for wind-damaged sites was registered in the SWIR region (B11, 23.5%) rather than in the NIR (B8A, 23,1%). Moreover, in wind-damaged site signatures, the absorption peak in the red was no more visible as reflectance grew linearly from B02 to B04. In the same period a slight increase in the signatures distance between the red-edge and the NIR overlapping region was observed; in April 2019, B05, B06, B07, B08, B8A mean reflectance in windthrown sites (12.1%, 17.1%, 19.6%, 21.9%, 23.1%) was higher than controls one (5.4%, 15.2%, 18.5%, 20%, 20.4%).
From July 2019, a significantly higher reflectance between the red-edge and NIR region (B05-B8A) was also registered in wind-damaged sites rather than in controls. In September 2019, the gap between the two signatures in the NIR had a further growth; B8 and B8A reflectance was 25.3% and 16% in wind-damaged sites and 6.4% and 2.9% in controls. In July and September 2019, a slight absorption in the red region reappeared, as an inflection point is visible in B03.
For both methods, wind-damaged areas were arranged in comparable classes of elevation. According to NDWI8A VID and RFS methods, the highest damage was detected from 1000 to 1500 m a.s.l. (49% and 66%, respectively), followed by 500-1000 m a.s.l. (33% for NDWI8A VID and 43% for RFS). Interestingly, at higher elevations (>1500 m a.s.l), NDWI8A VID detected higher damage compared to RFS (6% vs. 0.57%, respectively). Both NDWI8A VID and RFS methods had similar classes of slope aspects. However, the NDWI8A VID method had higher damage in north exposed slopes detrimental to south ones, when compared to RFS. Most damages were registered in the south (42% for NDWI8A VID, 52% for RFS) and east slopes (25% and 23%, respectively). For slope, again, damaged areas identified by the two methods are arranged in comparable percentual classes. More than 90% of the damage was detected in areas whose inclination is bigger than 40%. Finally, the methods showed discrepancies in classes of damages extension. According to the NDWI8A VID method half of the damaged surface is made up of a multitude of small areas (0.8-2 ha). On the contrary, RFS detected half of the damage in large areas (bigger than 25 ha).
Of the total 2545 ha delimited areas of NDWI8A VID, 661 ha were included in RFS estimations (26%), while the remaining 1884 ha (74%) fell outside from RFS perimeters. NDWI8A VID % error with respect to RFS estimations was 65.60%. Table S3 reports NDWI8A VID and RFS estimations in RFS delimited areas and NDWI8A VID % error as compared to RFS estimations for each municipality.

Map Validation
The Friuli Venezia Giulia windthrown areas, obtained through the NDWI8A VID technique, were visually compared to the ones provided by the Regional Forest Service (RFS) and overall accuracy was assessed by field surveys. The field survey was conducted on the 1st of February 2020 in Tagliamento Valley, in the municipalities of Villa Santina, Ampezzo and Forni di Sotto, where most damages were registered according to RFS. As shown by field surveys pictures available in Supplementary Materials (Figure S7), the VID method successfully delimited areal windthrows, being in most cases more precise than foresters' estimation and able to identify areas that were not included in RFS perimeters.

Vegetation Indices Temporal Profiles
In our study, NDWI and NDWI8A were identified as the best indices to separate winddamaged from undamaged sites, suggesting windthrow-induced water stress (R 2 = 0.80 and 0.77, respectively). Windthrow inevitably disrupted the soil-plant-atmosphere water continuum, decoupling the root system from soil or interrupting water movement through the xylem by damages at stem level [47]. The catastrophic failure of the water transport system led to water deficit in damaged plants immediately after the Vaia storm, and it was successfully detected using NDWI and NDWI8A indices right after the storm, through the abrupt rise of SWIR reflectance (Figure 2, Supplementary Table S1). Moreover, for most plants, canopy water content could not have recovered to control levels, as SWIR reflectance did not change in the following 11 months (Figure 4). This suggests that Vaia-event could have triggered extensive branch desiccation or, in the worst cases, tree mortality (forest die-off), as windthrows likely negatively affected the canopy water content both in the immediate and in the long-term period. This hypothesis is furtherly confirmed by the lower values of NDVI, which is highly correlated to plant productivity and vitality [44], calculated in damaged sites after the storm, Plant water status can be estimated by several physiological parameters, such as: (1) leaf water potential (i.e., water-free energy in plant organs), (2) leaf water content (usually expressed as the volume of water in leaves relative to the fully hydrated status) and (3) equivalent water thickness (EWT, i.e., the total water mass in the foliage per ground surface area) [48,49]. According to Hunt Jr et al. [50], EWT can be estimated by spectral reflectance at shortwave-infrared (SWIR) wavelengths, which are sensitive to water content. Moreover, indices calculated from a normalized difference of NIR and SWIR bands, such as NDWI, were found to be strongly related to canopy EWT using various sensors such as: Landsat TM, ASTER, AWiFS and SPOT [48,51]. Recently, Marusig et al. [52] found strong correlations between both midday leaf water potential and leaf water content and Sentinel-2 derived NDWI. Finally, special attention should be paid to NDWI8A which, despite the relatively low spatial resolution of 20 m, was the best index for this study purpose. This result confirms the usefulness of the narrow B8A at 865 nm in the NIR, designed to avoid contamination from water vapour through a high spectral resolution and yet successfully being able to represent the NIR plateau for vegetation.
Despite computing shorter wavelength bands (B8 and B4), the NDVI index successfully detected windthrows events (Table 3: R 2 = 0.76), as previously shown in other studies [26,45,53]. NDVI is the most consolidated index used to study the spectral behaviour of vegetation and it has been widely used as a proxy for green biomass and plant productivity [28,44]. The increase in reflectance registered in the red and blue regions (B02 and B04) in wind-damaged sites (Figure 4) suggests that windthrow inevitably affected also plants photosynthetic efficiency. Usually, under water-limited conditions, plants tend to prevent water loss by decreasing stomatal aperture (or by complete stomatal closure), which in turn also causes a decrease in the photosynthetic rates and allows plants to avoid excessive photosynthetic pigments excitation [54]. In our case, windthrow-induced water transport interruption triggered plant mortality thus affecting the photosynthetic activity from the outset. GNDVI was less efficient than NDVI in separating wind-damaged from undamaged classes (Table 3: R 2 = 0.60). The GNDVI was suggested to avoid the NDVI saturation effect in dense canopy structures and it uses the green band (B3) as a replacement for the red one (B04) [24]. However, the use of B03 rather than B04 has lowered GNDVI sensitivity to Vaia-induced changes in forests' photosynthetic efficiency.
In this sense, an interesting output is shown by PSRI performance in separating winddamaged and control sites ( Table 3: R 2 = 0.69). This index is sensitive to photosynthetic pigments relative content (both chlorophyll and carotenoids) and it has been widely used as a quantitative measure of leaf senescence [55,56]. When a plant starts dying, canopy discoloration occurs and immediate reduction in chlorophyll rather than in carotenoid concentrations decreases the spectral absorbance in the red region [33]. Moreover, reflectance in the red region increases significantly due to the production of anthocyanins, triggered by water stress [57]. The increase of reflectance in the red region rather than in the blue one proved that PSRI could be successfully used to monitor tree vitality and to determine the onset of wind-induced mortality. This also explains the index inverse pattern with respect to all the other indices, in PSRI low values were found in control sites while high values were associated with wind-damaged ones (Figure 3).
Despite it computes B12 band, which should correlate to changes in canopy water content, the SLAVI index was not as effective as NDWI and NDWI8A in detecting windthrows (Table 3: R 2 = 0.65). It is probable that the sum of red and SWIR bands, a proxy for photosynthetic efficiency and water content respectively, has increased the sensibility of the index in healthy forests. High seasonal fluctuations in SLAVI temporal profiles were registered in control values rather than wind-damaged sites, which maintain very low values ( Figure S3). NIR reflectance provides information about the cell structure of the plant and is widely used in vegetation monitoring [28]. However, in contrast with outputs in Einzmann [26], NIR reflectance did not change immediately after Vaia, indicating that the leaf structure of the damaged trees was not affected by the storm (Figure 4). This might be due to the fact that conifers usually occupy the conservative region of the Leaf Economic Spectrum (LES) [58], having higher LMA (leaf mass per area), greater lamina thickness, content of mechanical tissues, and thickness and sclerification of mesophyll cell walls [59], which guarantee longer leaf lifespan, with respect to broadleaves. Thus, the immediate post-storm period is likely excessively early to detect any leaves' structural deterioration. Besides the NIR, B06 and B07 red-edge bands of post-storm images were also not important for windthrow characterization (Figure 4). An exception is applied for the B05 band, whose values were significantly higher immediately after Vaia, thanks to its close relationship to the red region rather than NIR. In our study, IRECI and NDREDI, both using red-edge bands, were the least effective indices for wind-damage detection ( Table 3). The model including IRECI explained the 32% of data variance while regarding NDREDI, the model explained a minimal part of data variance (3%), proving to be incapable of discriminating damaged from undamaged stands. These results collide with those in Einzmann [26], where the red-edge was indicated as the most important channel for windthrow identification. These differences could be due to the fact that damaged sites in Einzamann [26] might have a higher proportion of deciduous species (e.g., Fagus sylvatica), or to the different season when the storm occurred (March for Storm Niklas while October for Vaia), and not least the different pixel resolution (5 and 10 m in Einzmann vs. our study, respectively). In this study, PSRI was the only index that includes a RE band to be suitable for the detection of stressed or damaged vegetation. However, as already mentioned, its suitability for windthrows detection was mainly attributed to changes in blue and red regions reflectance rather than in RE.
Together, spectral signatures ( Figure 4) and LMMs outputs (Table 3) showed that visible (B02, B03, B04), the first RE (B05) and short-wave infrared (B11, B12) were the most important bands to detect windthrows, as their values differed between damaged and un-damaged sites, even immediately after the storm. On the other hand, remaining RE (B6, B7) and NIR (B8 and B8A) band values did not differ between damaged and un-damaged sites. Indices which include one or more of these variable bands (B02, B03, B04, B05, B11, B12), together with one band which do not change from damage to un-damaged classes (B6, B7, B8, B8A), so-called "invariant", have been proved to be successful in identifying damaged sites. Similar results were recently obtained by Stych et al. [22], as they concluded that indices that take advantage of red, NIR and SWIR bands are the most sensitive and effective for the detection of wind disturbances in forests.

Spectral Reflectance of Vegetation after Windthrow
NIR reflectance did not vary between windthrown and control sites immediately after Vaia hit. A significant increase of B8 and B8A values in wind-damaged sites spectral signatures was registered 6 months after the storm, which corresponded to the resumption of the forest vegetative cycle (Figure 4). In the field of remote sensing of vegetation, it is common knowledge that NIR reflectance is generally lower for conifer crowns than broadleaf crowns due to different structural features which affect light scattering [60][61][62][63]. In view of the increase of reflectance in the NIR registered in spring (Figure 4), we suggest that windthrown sites may have undergone a secondary succession process. Consequently, this process might have caused a rapid recolonization of grasses and other herbaceous plants, several shrubs and many rapid-growing pioneer species. This recolonization has been facilitated by the tree removal which occurred in L1 (first two weeks of April 2019) and L2 (August 2019) batches. This was particularly evident in site L1, which was the first area subjected to batch cleaning. On 24 September 2019, the inspection in site L1 revealed that many typical light-requiring species were present in post-disturbance sites (e.g., Carex spp., Sambucus spp., Rubus ideaus, Pteridium aquilinum, Salvia glutinosa), as well as large populations of regenerating trees (e.g., Castanea sativa, Corylus avellana). This gradual colonization matches with the increased variability of all vegetation indices values except NDREDI, registered in site L1 on the 13rd of September 2019 (Figures 2 and 3, Figures S1-S5). The increase in index variability post-Vaia in C1 and C2 sites suggests that partial damage in these sites may have occurred (Figures 2 and 3, Figures S1-S6). It is probable that both control sites have been subjected to scattered windthrows, considering that their location, vegetation composition, lithology, slope, aspect and elevation is comparable to windthrown sites. These limited damages have been confirmed in the C2 site through the field survey; many scattered trunks and branches were present on the site forest floor. Interestingly, all vegetation indices except NDREDI also caught L1 and L2 batches tree removal; L1 values were almost comparable to L2 in March 2019, lower in April 2019, while L2 had lower values in September 2019 rather than in July 2019, when compared to L1 (Figures 2 and 3, Figures S1-S5).

NDWI8A VID Windthrow Maps
The NDWI8A differencing method successfully detected 2545 ha of forests affected by wind-caused damages, differently from The Regional Forest Service estimation, which indicated damage of 3183 ha. Of the 2545 ha NDWI8A VID delimited areas, 661 ha were included in RFS estimations, corresponding to a 65.60% error (Table S3). It is possible that manually digitalized perimeters have overestimated effectively windthrown areas for practical reasons. First, the Vaia storm did not impact forests uniformly within each RFS polygon. As qualitatively shown during field inspections, RFS included large forest patches that suffered from very scattered damage rather than intense areal damage. Secondly, the proposed method is not suitable for widespread damage as it requires intense and rather uniform damage to satisfy NDWI8A and extension thresholds. These considerations also explain the overall smaller wind-damaged area detected by the NDWI8A vegetation index differencing method than the one offered by the Regional Forest Service ( Figure 5). Windthrown area extension classes differ in the two methods: half of the areas detected by our method are small (0.8-2 ha) whereas the same quantity in RFS is represented by the biggest class (>25 ha). This difference may arise from the difficult conditions of the ground survey in RFS data, due to hard-to-access terrain and remote location of small disturbances. In fact, most areas are hardly accessible and located in remote and steep areas where field inspections and accurate manual perimeters may have been impossible to carry out. In this context, foresters may have preferred a wide perimeter strategy to include all possible damaged areas. On the other hand, according to Figure 5, the damage pattern highlighted following our method was similar to the one described by the RFS, particularly in relation to elevation, aspect and slope characteristics. The elevation, the high slope, and south/east main aspects of the windthrown areas match the Vaia storm evolution and interaction with the mountainous landscape of Friuli Venezia Giulia Alps and pre-Alpes. The Vaia scirocco gusts were channelled in the Alpine valleys causing most damages in the forests along steep mountain slopes where compression forces were amplified [11].
From an ecological perspective, the Vaia storm caused windthrows to many tree species such as Picea abies, Abies alba, Fagus sylvatica and species belonging to Alnus and Salix genera. Interestingly, many species considered as wind-resistant, such as beech, were subjected to windthrow suggesting that gusts' extreme speeds have by far exceeded plant-specific resistance. The vulnerability of Alnus and Salix sp. could have increased due to their exposed position and the presence of soils nearby the river expansion, whose humidity might have been increased by the strong precipitations associated with the Vaia storm.

Study Limitations
In this study, the NDWI8A vegetation index differencing method did not provide a prompt post-storm windthrow map to be used for Vaia impact assessment, as images acquired shortly after this extreme event were having a too low quality to be used as input. This is indeed the limitation of using optical remote sensing in winter storm-associated disturbances studies. Sentinel-2 data use has been strongly limited by shadows and snow presence, as well as cloud cover that occurred in the weeks near to the event [12]. However, as qualitatively shown by our field inspections, the NDWI8A differencing method provides satisfying and reliable windthrow maps. At last, most of the detected areas have intense and uniform damage which was not found in RFS large perimeters. The proposed method may offer a more accurate and easy way to delimit windthrows damages in the forests. Nevertheless, proper validation is still missing. Visual assessment and field surveys offer a broad screening that cannot be used for proper accuracy assessment. In the future, the accuracy assessment of NDWI8A VID-derived windthrow maps should be provided by extensive ground truth data. These data may be accurate forest use notifications data (which report possible occurred wind-damages), forest-stand updated datasets or high-resolution remote sensing data, such as Light detection and ranging (LiDAR) or orthophotos. Moreover, since clear-cutting and phenology effects can bring uncertainties, future studies should pay attention to choosing images for the proposed method. It should be recalled that VID chosen images were taken far away from the date on which the natural phenomena occurred when forests optimal sensing is guaranteed. The early spring period rather than autumn should be favoured when cleaning activities have not started yet.

Conclusions
Vegetation indices that highlight damaged-areas differences in the visible and the SWIR regions are strongly recommended for the detection of windthrows and for the monitoring of both forest die-off and succession dynamics. NDWI8A VID is an easy and cost-effective method able to produce reliable windthrow maps. Our results showed that the proposed method successfully delimited areal windthrows caused by the Vaia storm in forests throughout the Friuli Venezia Giulia Region. This study provides a safer alternative to detect the location of damaged areas, avoiding the risks associated with field inspections, especially in hardly accessible steep areas. Finally, the obtained results could be helpful in identifying sensitive areas to potentially uncontrollable insect/pathogen outbreaks. Such areas are advisable for sanitation as a preventive intervention to remove potential attractive host material for tree bark beetles and other pathogens [64].