2.1. Study Area and Sampling Design
The research was conducted in an old-growth forest in Central Amazon, located near the city of Manaus, Brazil (2°53'41"S, 60°16'26"W). This forest was impacted by a convective storm occurred in November 2015, which propagated destructive wind gusts and caused widespread tree-mortality across an area of ~70 hectares (
Figure 1a, b, and c).
The local relief in our study region is undulating, including plateaus, slopes and small valleys associated with perennial drainages [
39,
40]. The windthrown area is covered by a forest transitioning from
terra-firme dominating plateaus to
campinarana on the slopes and valley bottoms, such as described in adjacent areas [
41,
42,
43,
44]. In this type of forest, a high diversity of tree species can be found, for instance, from 63 to 137 species in a single hectare [
41]. The average annual temperature and annual precipitation in the region of Manaus (40 km from our study site) is 26.9 ± 0.17º C e 2.231 ± 118 mm (mean ± 95% confidence interval for period of 1970-2016, [
2]). This region has an evident dry season from July to September, with monthly precipitation less than 100 mm [
6,
9]. Soils on the plateaus have generally high clay content transitioning to sandy soils in the lower portions that are seasonally flooded [
40,
45]. In general, the soils have low fertility, low pH, high aluminum concentration and low organic carbon content [
39].
To quantify windthrow tree-mortality in the field we carried out a detailed forest inventory following a protocol built up from previous studies [
1,
2,
6,
7,
16]. The forest inventory was completed in two campaigns, conducted in December 2016 (~ one year after the windthrow) and April 2017 (~ one and a half years after the windthrow). We established three 20 m x 125 m transects, subdivided into 10 subplots of 10 m x 25 m (250 m
2;) each (total of 30 subplots, referred here as
field subplots).
Our transects crossed the entire disturbance gradient, ranging from areas with little or no disturbance (i.e., unimpacted old-growth forest) to severely impacted forest with few or no survivor trees. In the field subplots, we recorded, identified and measured the diameter at breast height (DBH) all living trees with DBH ≥10 cm. Due to logistics and safety of our field crew, and for reducing random errors due to possible missing trees hidden by the large amount of coarse wood debris typical of windthrows (
Figures S1–S4) [
1,
7], we did not count dead trees in severely impacted areas. Instead, we estimated the number of dead trees by subtracting the number of remaining living trees recorded in each subplot from the mean tree density of living trees in old-growth forests in our study region (i.e., ~ 600 trees.ha
-1, or 15 trees in 250 m
2,
Table S1, [
41,
42,
43,
46,
47,
48]). The subplots with more than 15 recorded living-trees were considered as undisturbed (i.e., no associated windthrow tree-mortality).
2.2. Spectral Mixture Analysis and Remote Sensing Estimates of Windthrow Tree-Mortality
We used Landsat 8, Sentinel 2 satellite imagery acquired before and after the studied windthrow event. For the WorldView 2 satellite, we used only data acquired after the event. Landsat 8 and Sentinel 2 were downloaded from the
Google Earth Engine platform (
https://earthengine.google.com/) [
49]. WorldView 2 was purchased from Digital Globe. We selected images with the minimum percent cloud-cover, the closest pass dates to each other and the greatest proximity between the acquisition date and the occurrence of the studied windthrow (
Table S2).
We used Spectral Mixture Analysis (SMA) [
50,
51] to quantify the fractions of endmembers [
20] following the routine of previous studies [
6,
7,
9,
16,
36]. We used the endmembers
photosynthetic vegetation (GV),
non-photosynthetic vegetation (NPV), and
shade (SHD) [
20]. The endmembers contain specific spectral signatures of multiple elements that make up the forest surface, and can be used to compute fraction-images on different targets of interest [
52,
53].
We used an analytical routine on the ENVI 5.3 software to extract the endmembers within each image in a standardized way [
20,
54,
55,
56]. We corrected the scales of Top of the Atmosphere (TOA) reflectance at the pixel level to allow a direct comparison between satellites. Sentinel 2 images are delivered in TOA values divided by 10
-4. Worldview 2 requires a band-by-band radiometric correction factor [
57] to bring pixel values to the same range as Landsat 8 and Sentinel 2.
As the evaluated satellites have different number of bands and respective spectral ranges, we used the bands common to all of them to obtain the spectral signatures of the endmembers. For Landsat 8 and Sentinel 2 we use the blue, green, red, NIR and SWIR bands. For Worldview 2, due to the smaller spectral range, we use the blue, green, red, red edge and NIR bands. The selected bands are sensitive to the physical, chemical and anatomical characteristics of the leaves and trunks, maximizing the distinction between GV and NPV [
50,
52,
53].
We used a spectral toolkit ([
58],
Figure S8) to select the purest pixels in all images (e.g. areas where all trees were downed versus unaffected old growth forest canopy) to acquire the most accurate/representative endmembers of interest [
20]. Finally, we conducted a SMA to compute the GV, NPV and SHD fraction-images [
51] covering our study area.
Windthrows have a specific spectral signature for the NPV fraction, which is originated from the deposition of large amounts of leaves, branches and trunks on the forest floor [
59]. To ensure that the GV fraction represented forest patches not affected by the studied windthrow, we selected endmembers within the adjacent old growth forest. The SHD endmember was selected from rivers in the vicinity of the impacted area. To generate final images, the NPV fractions were normalized without SHD as:
where
NPVnorm is the normalized values of the endmembers ranging from 0 to 1 [
20]. The best fraction images were selected according to [
60], but also based on field observations and the spatial distribution of
NPVnorm values with respect to windthrow patches visible in a RGB composition. We checked the histograms of the
NPVnorm fraction-images and we selected images that had more than 98% of the pixels with values between 0 and 1 [
61]. We also used the residual error (RMS) from the SMA as a selection criterium (i.e., the lower the better) (
Figure 2).
The differences between the
NPVnorm before (i.e., old-growth) and after disturbance (ΔNPV) provide a quantitative measure of tree-mortality due the windthrow [
36,
62]. The larger the ΔNPV values, the greater is the windthrow tree-mortality [
7,
16,
36]. To assess whether the NPV signal was influenced by previous disturbances such as selective logging [
61] and human-caused deforestation [
63,
64], for Landsat 8 and Sentinel 2 we calculated the ΔNPV fraction using a pair of images acquired before and after the studied windthrow date (
Table S6). This approach was not used for WorldView 2 as we did not have access to an image before the windthrow. Since
NPVnorm and ΔNPV for Landsat 8 and Sentinel 2 showed similar amplitude and range of values (
Figure S5,
Figure S6,
Figure S7,
Table S4,
Table S5,
Table S6), and a stable agreement, subsequent analysis with the three sensors were conducted based on the
NPVnorm.
We calculate a
NPVnorm value to each field subplot (
Table S3). For this, we converted
NPVnorm fraction images in raster format to
NPVnorm polygons in shapefile format using QGIS [
65]. After that, we overlayed the subplot shapefile with the
NPVnorm polygons for isolating all segments inside respective subplots. Thus, for each field subplot we obtained one or more segments of
NPVnorm. We calculated the area of each
NPVnorm polygon to obtain a weighted NPV value for each subplot [
1,
2,
4] as:
where
NPV is the weighted value within each subplot,
NPVnorm is the normalized value of each pixel and
A is the area of each pixel (m
2;) that is fully or partially included in the respective subplot; 250 is equivalent to the area (m
2) of each subplot.
Apart from the 30 field subplots, we estimated windthrow tree-mortality remotely for 100 virtual subplots (also 10 m x 25m; hereafter referred as
virtual subplots) randomly distributed across the disturbed forest (
Figure 1d and e). These virtual subplots were used to evaluate the robustness of tree-mortality estimates by the three satellites in adjacent areas containing a greater variation windthrow severity (
Figure 1e). The
NPV values for virtual subplots were obtained using the same method as described for field subplots.
2.3. Remote Sensing Estimates of Windthrow Tree-Mortality
We fitted a Generalized Linear Model (GLM) [
66] relating
NPV to the percentage of windthrow tree-mortality recorded in the forest inventory. The GLMs were fitted for each satellite using the binomial family for describing the distribution of residuals. We used the
logit link function, in which linear predictors must be reversed to the scale of the observations by means of an inverse function [
66].
We used the field subplots only for model fitting, thus avoiding overfitting problems for the cases in which the validation was carried using the same subplots [
67]. Virtual subplots were not used for validation/training. We rather focused on evaluating the performance of the models by partitioning the
Residual Deviance using Analysis of Variance (ANOVA). This allowed us to check which model best captured field measurements of windthrow tree-mortality (i.e., ground truth). We also used the Akaike Information Criteria (AIC, [
68]), Standard Error of the Estimates (S
yx), Root Mean Square Error (RMSE), Residual Standard Deviation (Sigma), and the Adjusted Coefficient of Determination (R
2;) computed with the Kullback-Leibler divergence formula (R
2;
KL, [
69]) as parameters for model selection. The R
2;
KL is suitable for exponential family models (e.g. logit), which retains the informative properties of the fit due to the inclusion of regressors [
69,
70]. The best model was the one with the lower Residual Deviance, AIC, S
yx, RMSE, Sigma, and higher R
2;
KL.
We evaluated the quality of remote sensing estimates of windthrow tree-mortality as metrics of precision/accuracy (i.e., minimum [Min], maximum [Max], median and mean) and dispersion/uncertainty (i.e., standard deviation [SD], standard error [SE] and 95% confidence interval of the mean [hereafter referred as 95% CI]). For comparing our results with those reported for other Amazon regions [
2,
4], we analyzed tree-mortality within categories describing windthrow severities: old-growth/undisturbed forest [≤4% of windthrow tree-mortality]; low windthrow severity [4% < windthrow tree-mortality ≤ 20%]; moderate windthrow severity [20% < windthrow tree-mortality ≤ 40%]; high windthrow severity [40% < windthrow tree-mortality ≤ 60%]; and extreme windthrow severity [> 60% of windthrow tree-mortality]).