## 3 Methods
### Preprocessing  
All the spectral bands contained in the S2 L2A product (Table S1) were first resampled to a 20 m spatial resolution. In order to mask out pixels from further processing and therefore reduce the impact of commission errors in the classification of burned pixels, a raster mask was generated for each single S2 scene using both the S2 L2A MUSCATE Quality Flags and S2 spectral bands (Table S2). Pixels corresponding to water areas were masked out using both the \"water_mask\" quality flag and a modified version of the water pixel formula [2]. Resampled and masked S2 L2A data were later used to compute spectral indices used for the identification of burned areas, specifically normalized burn ratio (NBR) [28], BAIS2 (burned area index for Sentinel-2) [2], and BAIS2\(\alpha\) (albedo-weighted burned area index for Sentinel-2) indices. NBR is a spectral index largely adopted for the identification of burned areas and calculated according to Equation (1). Here a slightly different formulation scaled by the value 1 is used in order to make it comparable to the other considered spectral indices. \[\text{NBR}=1-\frac{\text{B8}-\text{B12}}{\text{B8}+\text{B12}}. \tag{1}\]  
The BAIS2 index, calculated according to Equation (2), has already been demonstrated to be adequate for the detection of fire-affected areas [2]. It combines a band ratio in the red-edge spectral domain and aims to describe vegetation properties, with a band ratio similar to NBR, to detect the radiometric response of the short wave infrared (SWIR) spectral domain. \[\text{BAIS2}=\left(1-\sqrt{\frac{\text{B06}\times\text{B07}\times\text{B8A}}{ \text{B4}}}\right)\times\left(\frac{\text{B12}-\text{B8A}}{\sqrt{\text{B12}+ \text{B8A}}}+1\right). \tag{2}\]  
The BAIS2\(\alpha\) index, presented in this research study and calculated according to Equation (3), is a modification of the BAIS2 index that uses the surface albedo as scaling information. \[\text{BAIS2}\alpha=\left(1-\sqrt{\frac{\text{B06}\times\text{B07}\times\text{ B8A}}{\text{B4}}}\right)\times\left(\frac{\text{B12}-\text{B8A}}{\sqrt{\text{B12}+ \text{B8A}}}+1\right)\times\left(1-\sqrt{\alpha}\right). \tag{3}\]  
Surface albedo \(\alpha\) is computed from bottom-of-atmosphere S2 reflectances using a set of pre-computed band weighting coefficients (reported in Table S1) [52] determined according to the spectral response function of the MSI sensor. Similar to spectral indices, surface albedo has been seen to lead to discrimination capability for burned areas [33]. The biophysical processor [53; 54] was used to compute leaf area index (LAI) from L2A atmospherically corrected data after resampling and masking. Difference in pre-fire and post-fire LAI was later compared to fire severity estimated from spectral indices. Multitemporal observations of the NBR, BAIS2, BAIS2\(\alpha\), and LAI were finally stacked in a multi-dimensional datacube, generated for each granule. Since S2 data acquired during the year 2017 belong to different processing baselines, exhibiting variable geolocation spatial accuracy, each pixel revealed a shaky multitemporal profile as a result of the weak spatial coherence [55]. For this reason, an image co-registration step, performed separately for each S2 granule, was included in the processing chain in order to identify the horizontal linear shift for each observation date and to co-register the multitemporal S2 observations. First, a master S2 image was selected from all the available observations in 2017 according to the following criteria: (i) be acquired by the Sentinel-2A satellite; (ii) belong to the baseline 02.05 (since changes in the processing baseline could modify image geolocation, and baseline 02.05 includes all the S2 data acquired in 2017 from 27 April to 22 October 2017 [56]); (iii) exhibit the highest number of unmasked pixels. Later, for each single S2 observation, nine combinations of horizontal linear shifts, from 0 to 1 pixel (20 m resolution) in both X and Y axes, were applied in order to generate nine slave images.