Assessing Typhoon-Induced Canopy Damage Using Vegetation Indices in the Fushan Experimental Forest, Taiwan

Cyclonic windstorms profoundly affect forest structure and function throughout the tropics and subtropics. Remote sensing techniques and vegetation indices (VIs) have improved our ability to characterize cyclone impacts over broad spatial scales. Although VIs are useful for understanding changes in forest cover, their consistency on detecting changes in vegetation cover is not well understood. A better understanding of the similarities and differences in commonly used VIs across disturbance events and forest types is needed to reconcile the results from different studies. Using Landsat imagery, we analyzed the change between preand post-typhoon VI values (ΔVIs) of four VIs for five typhoons (local name of cyclones in the North Pacific) that affected the Fushan Experimental Forest of Taiwan. We found that typhoons varied in their effect on forest canopy cover even when they had comparable trajectories, wind speeds, and rainfall. Most VIs measured a decrease in forest cover following typhoons, ranging from -1.18% to -19.87%; however, the direction of ΔVI–topography relationships varied among events. All typhoons significantly increased vegetation heterogeneity, and ΔVI was negatively related to pre-typhoon VI across all typhoons. Four of the five typhoons showed that more frequently affected sites had greater VI decreases. VIs ranged in their sensitivity to detect typhoon-induced changes in canopy coverage, and no single VI was most sensitive across all typhoons. Therefore, we recommend using VIs in combination—for example Normalized Difference Infrared Index (NDII) and Enhanced Vegetation Index (EVI), when comparing cyclone-disturbance-induced changes in vegetation cover among disturbances and across forests.


Introduction
Tropical forests harbor a large portion of global biodiversity [1] and play a key role in global carbon cycling [2][3][4]. An estimated 52% of the terrestrial carbon pool is contained by tropical and subtropical forests, although they cover less than 15% of the land surfaces excluding Antarctica [5]. Thus, changes in the geographic extent and quality (i.e., structure) of tropical forests have major effects on biodiversity conservation efforts [6,7] and important ecological processes such as carbon sequestration [8,9]. Disturbances such as fires, pathogens, and tropical cyclones are important agents affecting the dynamics of many tropical and subtropical forests [10][11][12]. Tropical cyclones are among the most common natural disturbances to tropical and subtropical forests [13]. Cyclones vary in their effects on forest structure and damage to individual trees, with damage ranging from defoliation and branch-stripping to bole-snapping and the uprooting of trees [14]. These events have consequences for forest functioning through the changing of canopy cover, gap formation, and tree population demography [12,15].
Spatial patterns of damage are key to understanding forest ecosystem dynamics in relation to cyclone disturbance at the landscape scale. In addition to vegetation characteristics such as forest type, species, phenological and successional state [16][17][18][19][20][21][22], and the distance to cyclone eye [23][24][25][26], landscape topography is also a main factor which modulates cyclone disturbance effects on vegetation [27,28]. Cyclones tend to have greater effects in exposed areas, such as ridges or windward slopes [25,26,[29][30][31][32][33]. In addition, cyclone damages may vary with elevation [31,34], slope steepness [20,32] and topographic position [18,20,29,32]. The interactions among the topographical variables, biotic characteristics, and wind properties lead to complex patterns of cyclone effects across the landscape [25,32,35,36]. Thus, the spatial patterns of cyclone effects are likely to vary among different cyclone events. In fact, dissimilarities in ecological effects among multiple cyclones have been reported for various sites, including when several typhoons (local name of tropical cyclones in the Northwest Pacific) have occurred in a short time period with similar tracks [21,25,[37][38][39]. Because topographical characteristics are key to cyclone effects, to better predict cyclone effects on forest landscapes it is important to examine which topographic properties have consistent relationships with cyclone disturbance and which ones are more variable. Few studies have compared landscape damage patterns of multiple cyclones [36,40], possibly due to the scarcity of cyclones over a period of time during which the changes in vegetation patterns are minimal except for those caused by cyclones in most regions.
Cyclone effects on vegetation can be assessed through direct field observations (e.g., [16,41,42]) or the use of remote sensing techniques (e.g., [23,32,43]). Remote sensing offers the possibility of assessment of cyclone effects over broad areas [23,32,44], and the results have been used to formulate conservation recommendations [45][46][47]. Used in combination with other data, remote sensing can effectively identify broad spatial patterns of cyclone damage in relation to vegetation types and topographical properties [20,32,43,47] that could be otherwise difficult to assess with ground surveys. Some of the most commonly used remote sensing data are vegetation indices (VIs), which are derived from satellite or aircraft-collected spectral reflectance measurements.
Many VIs have been developed to assess landscape vegetation patterns and dynamics [48], with each having its strengths and weaknesses [49]. VIs based on measurements using near-infrared (NIR) spectral bands, such as the Enhanced Vegetation Index (EVI), the Normalized Difference Vegetation Index (NDVI), and the Soil-Adjusted Vegetation Index (SAVI), have been widely used in vegetation assessments [48,[50][51][52] although the shortwave infrared band of Landsat has also been shown to be effective for evaluating photosynthesis and forest canopy cover [53,54]. Relationships between VIs and forest characteristics are diverse and can be context-dependent; however, generally, the EVI is more closely related to canopy structure [55] and has a better relationship with photosynthesis than NDVI under high leaf area index (LAI, [56]) because EVI is less likely to saturate than NDVI at high vegetation cover [57][58][59]. The Normalized Difference Infrared Index (NDII), another NIR-based index, has been shown to be more sensitive than the NDVI to changes in vegetation water content [60], particularly at high vegetation cover [61,62]. The SAVI, also an NIR-based index, could be an option when the influence from exposed soil is a concern, such as in disturbed areas with low vegetation cover [63]. Disparate uses of VIs in studies of forest disturbances make cross-study comparisons difficult. Hence, a better comprehension of the comparability among frequently used VIs (i.e., NDVI, EVI, NDII) in relation to forest disturbance and recovery is needed for studies of vegetation disturbance ecology using VIs. Given the widespread occurrence and frequency of cyclones, they provide an opportunity to comparatively evaluate the effectiveness and consistency of VIs to capture vegetation disturbance and recovery within and across cyclones and forests.
The Fushan Experimental Forest (FEF) in north-eastern Taiwan provides a unique opportunity  to study and compare the effects of multiple typhoons, as it was hit by nearly one typhoon per year  between 1951 and 2005, including many years with multiple typhoons of different trajectories [42,64,65]. The frequent typhoon occurrence in this site has permitted the study of disturbance representativeness of large forest plots over a wider landscape [66], but understanding of damage variation and the comparability of different VIs across typhoons is lacking. Here, taking advantage of the FEF location, we use remotely sensed data to examine: 1) the consistency of patterns of typhoon damages among different typhoons; 2) the spatial patterns of typhoon damage in relation to cyclone properties (i.e., wind direction) and site topographical features (e.g., elevation, slope, and aspect); 3) the relationships between typhoon frequency and severity of typhoon damage; and 4) whether patterns of typhoon damage are consistent among different vegetation indices.

The Fushan Experimental Forest and Typhoons
The 1000-ha FEF is a nature reserve located in northeastern Taiwan, approximately 25 km from the east coast of the island ( Figure 1). Elevation ranges from 400 to 1400 m above sea level (asl). From 1993 to 2004, mean annual temperature measured 18.2 °C, mean annual precipitation measured 4271 mm, and mean relative humidity was 95% [67,68]. The forest is classified as an old-growth submontane evergreen broadleaf forest dominated by trees species of Lauraceae and Fagaceae [67][68][69]. The FEF is regularly affected by typhoons between June and October, with an average of 0.74 typhoons per year between 1951 and 2005 [64]. Five recent typhoons (Herb, Nari, Aere, Soudelor, and Dujuan) were selected based on the availability of Landsat images with less than 50% cloud cover of the FEF within seven weeks of the typhoon passage ( [23,72] but see [44,73]). We chose the temporal limit of seven weeks in order to minimize any changes in leaf phenology, which can occur within a matter of weeks. Although 50% cloud cover is substantial, it is difficult to find images with low cloud cover in cyclone-impacted regions [44], especially for the FEF, where it rains more than 220 days per year on average [74]. The five typhoons passed within 100 km of the FEF, and were all category 2 or 3 on the Saffir-Simpson scale ( [75], Figure 1). Given that each of the typhoon air masses were within 100 km of the FEF, a distance envelope within which winds are strongest, the FEF was affected by all five typhoons.
According to the IBTrACS cyclone tracks dataset [70,71], Typhoon Herb (1996, category 3) was the first typhoon of its season to cross the 100-km range envelope of the FEF. In contrast, Typhoon Nari (September 16th 2001) was preceded by Typhoon Toraji, which passed in proximity to the FEF on July 29th-30th (categories 3 to 1), and Tropical Storm Mindulle (July 4th to 5th 2004) preceded Typhoon Aere (August 8th 2004). No typhoon made landfall before Soudelor and Dujuan (2005) near the FEF, and the two typhoons were seven weeks apart ( Figure 1). Among the studied typhoons, Typhoon Aere did not make landfall, but its center was 80 km from the FEF at its closest point. All typhoons were associated with high rainfall (690-1300 mm), and winds over 20 m s -1 at the FEF (Table  1). Nari was the wettest typhoon, whereas Herb had the highest wind speed (Table 1). The winds of the typhoons had different trajectories ( Figure 2). The winds of Soudelor and Dujuan were mostly out of the south and west, whereas the winds of Nari and Aere were of northern origin ( Figure 2). Although within each typhoon, at the hourly timescale, wind direction and speed were largely consistent, and daily data showed that the strongest wind gusts could come from different directions. , and direction of daily strongest wind during five typhoons that affected the Fushan Experimental Forest. Hourly mean wind direction data was not available for Typhoon Herb. Percentages do not necessarily sum up to 100% as a wind speed of 0 m s -1 is not shown. Data are from the Central Weather Bureau of Taiwan and were analyzed using the 'openair' R package [76].

Satellite Images
Landsat images were used to analyze vegetation cover change and recovery within the FEF in relation to typhoon disturbances. For each typhoon, pre-and post-disturbance images were taken within seven weeks of the typhoon passage to minimize phenological change in vegetation cover [77,78]. Images used to study vegetation recovery were taken during the growth season before typhoon passage, and one year later during at a comparable time with minimal cloud cover. However, the recovery associated with Typhoon Herb was not studied because of the lack of images with low cloud cover. Basic information on the images used in this study is given in Table 2, images used to study pre-and post-typhoon states are shown in Figure S1 in the Supplementary Material. Data from satellites Landsat 5, 7, and 8 were downloaded from the EarthExplorer website [79] as surface reflectance at 30 m resolution, with atmospheric and radiometric corrections performed with the LEDAPS (Typhoons Aere, Herb and Nari) or LaSRC (Typhoons Dujuan and Soudelor) algorithm. A digital elevation model (DEM) at a 30-m spatial resolution was downloaded from the ALOS World 3D dataset from JAXA. Non-forested surfaces were identified with the Global Forest Cover dataset version 1.5 (GFC, [80]). Table 2. Basic information on Landsat images used in this study. Null cells were either covered by clouds or defined as non-forested using the Global Forest Cover dataset [80]. Pre-and postdisturbance images are shown in Figure S1.

Pre-Processing
Pre-processing of spectral data followed [81]. Rasters were topographically corrected using the C correction method with the topcor function of the 'RStoolbox' package [82] in R 3.6.1 [83]. Clouds and shadows were identified based on the Pixel QA band provided by the USGS, and derived from [84] with CFmask detection algorithm and visual inspection of true color composite (i.e., for Nari recovery) before being removed from all rasters. In addition, non-forested surfaces were identified and removed using the GFC dataset with a threshold of 75% following studies of other tropical moist forests [85,86].

Processing
Three topographical variables were derived from the DEM using the terrain function from the 'raster' R package [87]: slope steepness, slope aspect, and topographic position index (TPI). The TPI is derived from the type of terrain surrounding a particular cell, with values ranging between -10 and 10. A TPI value of 0 indicates flat terrain surrounds the cell, whereas a cell surrounded by lower terrain has a positive TPI and a cell surrounded by more elevated terrain has a negative TPI. Aspect numerical values were converted into eight categorical variables (e.g., north, northeast, etc.) covering 45° each, centered on the cardinal direction (e.g., from 22.5° to 67.5° for northeast). TPI values were converted into six slope positions following [88] as described in Table 3. Slope steepness was used to differentiate middle slope from flat slope areas, because the two had the same TPI value range (Table  3). Table 3. Conversion thresholds of topographic position index (TPI) to slope positions and based on standard deviation (SD) of TPI and slope steepness following [88].  [90]). Although NDVI has been widely used in ecological studies (see reviews by [48,50,51]), it saturates faster than EVI and NDII at high biomass densities [55]. In contrast, EVI is more sensitive than NDVI to vegetation changes in areas with high biomass [91]. NDII has been shown to provide better monitoring of canopy defoliation and damage to forest structure [92][93][94], while SAVI reduces the soil effect [90]. The four indices are the products of four bands from the Landsat sensors: blue (B), red (R), NIR, and the first short-wave infrared band (SWIR1).
To measure the effects of typhoon disturbance, the variation in VI was calculated using Equation 5. A negative ΔVI indicates a decrease in vegetation cover post-disturbance.

Analysis of Disturbances among Vegetation Indices
Typhoon effects on VIs were first analyzed separately for each of the five typhoons, in which only cloudy cells in the images related to the given event were excluded from the analysis. In contrast, only the shared non-clouded cells were used for cross-typhoon comparisons. For each typhoon, the ΔVI was compared to 0 with a one-sample Wilcoxon signed-rank test. Then, Spearman's ρ was used to explore the relationships among the four ΔVIs (very weak < 0.2 < weak < 0.4 < moderate < 0.6 < strong < 0.8 < very strong) with the corr.test function from the 'psych' R package and p adjustment of Bonferroni for multiple comparisons [96]. In addition, the same test was used to study the correlations between pre-disturbance state ( ) and the ΔVIs to examine if the changes in VIs are related to pre-disturbance vegetation conditions. Across all typhoons, correlations between individual ΔVI (e.g., ΔNDVI among the five typhoons) were tested with Spearman's ρ to examine the consistency of vegetation damage detection among different VIs. Furthermore, to assess the effects of each typhoon disturbance on vegetation, coefficients of variation (CV) of each VI before and after disturbances were compared through a bootstrapped comparison of means (5000 iterations) using Equation 6.

Typhoon Damages and Topography
The relationships between ΔVIs and topographical variables were analyzed with ordinary leastsquares (OLS) regression models, which are commonly used in forest ecology [97][98][99], including elevation and slope as continuous parameters, and TPI and aspect classes as ordinal variables. The number of cloud free pixels across all the images used for model construction was 4593.

Disturbance Frequencies and Intensity
Damage frequency was defined as the number of typhoons (out of 5) affecting a particular cell. Small variations in VI unrelated to wind damage may cause negative values; therefore, a more conservative threshold than < 0 was used to identify typhoon-damaged cells (Equation 7). Only ΔNDII and ΔEVI were used in this analysis, whereas ΔNDVI and ΔSAVI were not included because of their moderate-to-strong correlations with ΔNDII and ΔEVI (ρ values ranging from 0.48 to 0.97). To assess if more frequently disturbed cells (i.e., those with frequencies of 4 or 5) were also more severely disturbed (had more negative ΔVI), mean ΔVIs among the six frequency classes were compared using the multiple comparisons procedure described by [100] for unbalanced designs with unequal variances with the 'multcomp' [101] and 'sandwich' [102,103] packages. Similarly, TPIderived classes as well as aspect classes were compared to examine if cells containing certain topographical features were more frequently affected by typhoons. Relationships between elevation and slope with damage frequencies were tested using OLS regressions.

Vegetation Recovery
Forest canopy recovery following Typhoons Aere, Nari, and Soudelor-Dujuan was analyzed by comparing an image of the FEF before disturbance to an image taken a year later during the same season (i.e., after disturbance, Table 2) through bootstrapped comparisons on means (5000 iterations). Typhoons Soudelor and Dujuan, which were only seven weeks apart, were studied together, as it was not possible to separate their recovery. Typhoon Herb could not be included in the recovery analysis due to the lack of images with low cloud cover following the typhoon.

Vegetation indices, typhoons, and the effect of prior vegetation cover
The four VIs all significantly decreased after the five typhoons (one sample Wilcoxon, V = 63407 to 42679021, all p-values < 0.001, sample size = 6361 to 10128) except for the significant increase of NDII, EVI, and SAVI associated with Typhoon Herb (all p-values < 0.001, Table 4). Compared to the pre-typhoon VI values, the decrease was greatest for Dujuan and smallest for Aere and Soudelor (Table 4).  5). Averages are calculated on variable portions of the Fushan Experimental Forest as the cloud cover varies among the images used to study typhoon events (see Table 2). Correlations among ΔVIs were significant for all typhoons, although they varied in strength ( Table  5). All the correlations were positive except for those between ΔEVI-ΔNDVI and ΔEVI-ΔNDII. Most of the correlations were weak to moderate (ρ = 0.05-0.50) but the correlations between ΔEVI and ΔSAVI were generally strong (ρ up to 0.90) ( Table 5). Pre-disturbance VIs and ΔVIs had significant negative correlations in all cases (p < 0.001, Table S1). These relationships were in general weak for Typhoons Dujuan and Nari and considerably stronger for Typhoons Aere, Herb and Soudelor (Table S1). Table 5. Spearman's correlations among the four ΔVIs measured over five typhoons for the Fushan Experimental Forest. All correlations are statistically significant (p-values < 0.01). ΔVIs are calculated as the difference between post-and pre-typhoon values (Equation 5). The relationships are studied on variable portions of the Fushan Experimental Forest as the cloud cover varied among the images (see Table 2).

Variation of ΔVIs among Typhoons
The correlations of ΔVIs among typhoons were mostly negative between Typhoon Dujuan and other typhoons and mostly positive among other typhoons (Table S2). Most of the correlations were either weak, very weak or not significant. However, moderate positive correlations were detected for the ΔNDVI of Typhoon Herb-Soudelor, Herb-Nari and Soudelor-Nari (Table S2).

Effects on Vegetation Heterogeneity
Except for SAVI in relation to Typhoon Herb, typhoons led to a higher heterogeneity in VI values, as indicated by the significantly greater post-typhoon CV of all VIs except for NDVI for Typhoons Aere and Dujuan, wherein the CV did not change significantly (Figure 3).

Topography and Disturbance Severity
Topographical variables differed in their ability to explain the variation of ΔVIs among the typhoons (Table 6 for Nari and Herb, Table S3 for other typhoons). Topography was a better predictor of ΔNDVI for typhoons Herb and Nari (adjusted R 2 = 0.52 and 0.47, respectively) than it was for other ΔVIs (adjusted R 2 ≤ 0.34). For typhoons Aere, Dujuan, and Soudelor that passed further from the FEF (Figure 1), topographical variables were poor predictors of ΔVI values (adjusted R 2 ≤ 0.16). Elevation and slope were significant in explaining ΔVI variation in all cases except slope for ΔEVI and ΔSAVI for typhoon Dujuan (p > 0.05, Table 6 and Table S3); however, coefficients were small (β < 0.0002 for elevation, β < 0.008 for slope). The magnitude of ΔVI decreased with increasing topographic slope for all VI-typhoon combinations, except for ΔEVI with typhoons Herb and Nari, and ΔSAVI with typhoon Herb. However, increasing elevation led to either higher or lower damage (i.e., change in VI) depending on the typhoon in question, but the change in direction was consistent among VIs for each typhoon. Except for flat slopes, all TPI positions had positive and statistically significant (p < 0.05) regression slope coefficients for Soudelor, Nari, Herb and Aere. On the other hand, Typhoon Dujuan showed negative relationships for all TPI positions except for flat slopes. Nevertheless, the sign of regression slope coefficients remained consistent among ΔVIs for each typhoon. The relationships between aspect and ΔVI values changed among typhoons, showing no clear pattern. Table 6. β Coefficients of ordinary least squares linear models between topographical variables and changes in vegetation index values (ΔVIs) associated with Typhoons Nari and Herb. Sample size of 4593 (equal to the number of pixels from analyzed Landsat images). Significance levels are shown with † (p < 0.05), ‡ (p < 0.01), and * (p < 0.001). Only the first non-zero value is shown, complete results for the five typhoons are in Table S3.

Disturbance Frequency and Severity
For all typhoons but Dujuan, the more frequently affected cells had greater VI losses (Figure 4). Typhoon damage frequency varied with elevation and slope steepness. Elevation was negatively related to EVI-based damage frequencies (β = -6.71, p < 0.001, adjusted R 2 = 0.016) but not to NDIIbased frequencies (p = 0.47). On the other hand, slope was positively related to EVI-and NDII-based frequencies (β = 0.32, adjusted R 2 = 0.001, p = 0.01 for EVI; β = 2.16, adjusted R 2 = 0.07, p < 0.001 for NDII). The proportion of damage frequency varied among aspects ( Figure 5, Table S4). For NDII, all aspects were significantly different in damage frequency (p < 0.001), except for south-east, northeastnorth, southwest-northeast, southwest-north, and west-northwest (p > 0.05, Table S4). For NDII, southwestern to northwestern aspects had the highest damage frequencies and southeastern aspects had the lowest disturbance frequency ( Figure 5A). Considering EVI, northern aspects had more frequent typhoon damage than western and southern aspects, whereas there were no significant differences between other aspects (Table S4, Figure 5B). Among slope positions, according to EVI, the flat slope, lower slope, and valley positions had similar typhoon damage frequencies (Table S4); these three topographic positions were more frequently damaged than other positions. For EVI, ridge areas had lower typhoon damage frequency than all other slope positions (Table S4). With NDII-based frequencies, no significant differences were observed between rides and slope or among other topographical positions (p > 0.05 with adjustment, Table S4).

Recovery
Vegetation cover had recovered in less than a year after Nari, as average VI values measured in June 2002 were equal to or greater than their values in June 2001 (positive 95% CIs, Figure 6). Similarly, all Vis, except NDII, registered vegetation recovery in less than one year after typhoons Soudelor-Dujuan ( Figure 6). In contrast, EVI was the only VI showing recovery after Aere ( Figure 6).

Figure 6.
One-year regeneration of the Fushan Experimental Forest as shown by the vegetation indices measured before typhoon passages and a year later at the same season. Regeneration of Typhoons Soudelor and Dujuan were merged as they passed during the same season; no satisfying images were available to study Typhoon Herb regeneration. All the differences between pre-typhoon and posttyphoon images are significant based on 95% confidence intervals measured through bootstrapped comparisons on means.

Consistency in the Damage Effects among Typhoons
The consistent decrease in VI values following all five typhoons (with very few exceptions; see Table 4) suggest that all the VIs can generally capture typhoon-induced losses in vegetation cover [23,24,43,[104][105][106]. However, we find weak correlations in ΔVIs between different typhoon events (ρ < 0.4, Table S2), indicating that typhoons range in their effects on vegetation cover. Although this is not surprising because typhoons differ in intensity, duration, trajectory and occurrence time relative to plant phenology, the inconsistency among vegetation indices highlights that results derived from one or a few disturbance events are unlikely to represent general trends in disturbance effects. Indeed, our results show that the successive Typhoons Soudelor and Dujuan did not have consistent effects on vegetation although they had comparable paths, wind speeds, and directions ( Figure 1-Figure 2, Table 1). Additionally, the very weak negative correlations of the ΔVIs between typhoons Soudelor and Dujuan (Table S2) indicate that the two typhoons had different effects on the vegetation cover, as observed in other sites subject to successive cyclones [21,22,107].
Powerful cyclones, such as Hurricanes Hugo and Maria in Puerto Rico, Typhoon Herb in Taiwan, and Cyclone Larry in northeastern Australia, attract scientific study of their ecological effects [25,26,35,38,40,43,[106][107][108][109][110][111]. This high prevalence of studies has advanced our understanding of cyclone ecology, especially in relation to the most powerful and damaging storms (reviewed by [12,14,112]). However, in this study, Typhoon Herb was considered to be the most powerful typhoon in several decades [113]. Typhoon Herb had the greatest wind speed among the five typhoons studied, while Nari was most intense in terms of precipitation. Considering these two storms, the resultant changes in vegetation index values (ΔVIs) were only weakly correlated. Future increases in the frequency of the most intense cyclones are predicted [114][115][116]. However, our results suggest that conclusions drawn from studies that document the effects of a single or a few intense cyclones are insufficient for predicting the effects of future cyclones.
The fact that there were small or no changes in VIs associated with Typhoon Herb is somewhat surprising. In addition to the potential influence of image quality on our analyses, the timing of Typhoon Herb is probably the most important factor. Typhoon Herb occurred in the summer (late July) when the forest was in the middle of the main growing season, and the two images were approximately two months apart. Thus, plant growth could be substantial during the period, potentially confounding the detection of typhoon-caused decreases in vegetation cover by VIs. This contrasts with Typhoon Dujuan, which caused the largest decreases in VIs. Typhoon Dujuan occurred near the end of the main growing season (April to September) so that although the two images were also approximately two months apart, there was little vegetation growth during the period. As a result, typhoon-induced vegetation loss can be better detected with the VIs. Thus, image timing in relation to plant growth phenology should be considered when examining disturbanceinduced changes in vegetation cover using VIs.
Not only were the overall effects inconsistent among typhoons, there was large variation in linear regression results, which showed that typhoon damage-topography relationships were inconsistent among typhoons (Table 6 and Table S3). Topography was better at explaining disturbance distribution across the landscape for Typhoons Herb and Nari, the typhoons that passed the closest to the FEF, but they were not the storms that caused the greatest degree of vegetation damage. This result suggests that topography-vegetation damage relationships vary with cyclone distance and that topography is a key determinant of vegetation damage only when typhoons are very close to the study site, despite the magnitude of the damage. It also suggests that factors other than cyclone distance determine the severity of typhoon-induced vegetation cover damage.
Nevertheless, there were some consistencies across typhoons. First, the canopy generally recovered quickly as many VIs returned to their pre-disturbance values within a year ( Figure 6). This result is consistent with the report of the rapid recovery of the FEF observed following Typhoon Bilis using NDVI [106]. However, recovery of a VI does not imply total canopy recovery as LAI and litterfall typically do not recover within a year of damage in the FEF [73]. Second, the relationship between high pre-disturbance NDVI and strong NDVI loss observed by [106] was also detected for all five typhoons in this study despite their differences in paths and intensities (Table S1). This pattern may be explained by the higher aerodynamic drag of dense canopies as suggested by [16], who reported a similar relationship between pre-disturbance LAI and LAI loss (see also [17]). Third, cyclones are disturbance agents which induce heterogeneity in forest landscapes [14,64]. Secondary tree falls and defoliation may have led to the increased heterogeneity observed following almost all typhoons examined here. Finally, most ΔVIs were positively correlated between typhoons (Table S2), except for Soudelor-Dujuan, suggesting that different typhoons could have similar effects although the strength of the relationships varied greatly among typhoons and VIs.

ΔVIs in Relation to Topography
The topography only explained a small proportion of the variation in vegetation cover change. As observed by [32], steeper slopes were associated with greater decreases of all VIs across all typhoons (Table 6 and Table S3), perhaps because of different wind exposures and soil stability (e.g., landslides, [117]). However, our results show that the relationships between vegetation damage and other topographical variables changed among typhoons and VIs, even when typhoons had similar tracks and wind directions. Thus, observations from a single event should be generalized with caution, as they do not necessarily remain true for other cyclones. Indeed, complex interactions between wind and topography have occurred in the FEF, as leeward positions were more affected than eastern positions for Typhoon Nari but not for Dujuan ( Figure 2, Table 6 and Table S3). Less exposed slopes may also offer little protection if the cyclone is particularly strong [40]. The overall low adjusted-R 2 values indicate that factors other than topography (e.g., vegetation conditions and timing of typhoon disturbance) likely play a more important role in determining typhoon-induced changes in vegetation cover.

Typhoon Disturbance Frequency
The greater damage for more frequently than less frequently affected cells (Figure 4) is consistent with a study from a moist forest of Puerto Rico, in which trees sustaining heavy damages from Hurricane Hugo (1989) were more likely to be damaged again nine years later by Hurricane George [118]. The result is also consistent with the mostly positive correlations between ΔVI values among the different typhoons (Table S2). Although it is possible that sites with certain topographical features within the FEF are more prone to typhoon damage, this claim is not supported by the lack of consistent topographical damage patterns among typhoons ( Table 6 and Table S3). Possibly, susceptibility to damage is more related to weakened vegetation in these sites (i.e., repeated disturbances). Further field-based research may help to separate the effect of forest characteristics, topography, and typhoon frequency on canopy damage at the FEF.
The greater elevational signal in typhoon damage frequencies at lower than at higher elevations based on EVI fits the pattern of greater damage frequency at lower topographic positions (lowerslope, valley) than other positions (e.g., middle-and upper-slope). It also supports results from field observations of greater typhoon effects at low than at high elevations along a 2300-m elevational gradient in central Taiwan [65]. It appears that the elevational pattern is consistent across spatial scales including variation in slope across the landscape, possibly because the intensity of typhoons declines quickly in rough topography as they move upward [65]. The differences in NDII and EVI sensitivity to vegetation characteristics may explain the different relationships of their damage frequencies with slope. As observed in Puerto Rico [36], different aspects had different disturbance frequencies, with northern aspects having significantly higher values. However, the relationship with the windward-leeward direction was less clear here, as western aspects (leeward in the FEF) also had high damage frequencies. The very rough topography at the FEF, with a mean slope of 38%, probably obscures relationships between slope, aspect and typhoon damage.

Consistency among Vegetation Indices
We did not detect any VI that was consistently the most sensitive across the five typhoons. EVI and NDII were considered particularly functional to measuring vegetation characteristics such as canopy structure [55] and water content [61,62]. However, in our study, they did not detect strong decreases in vegetation cover relative to pre-typhoon values for Typhoon Herb. Such differences may be the product of different damage variations across disturbance events, as reported for other sites [21,22,107]. It may also be due to the phenological change in leaf properties between pre-and postdisturbance images [77,78], although effort was made to minimize such variation by selecting images from within seven weeks of typhoon events.
Surprisingly, NDVI was the only VI that showed a decreasing canopy cover following Typhoon Herb, whereas it was the least practical index for Dujuan and Nari, despite its lack of sensitivity under high LAI [57][58][59]. NDVI was also the only VI that did not detect change in vegetation cover heterogeneity following Typhoons Aere and Dujuan ( Figure 3). Overall, NDVI is less sensitive to vegetation cover change in the FEF than the other Vis, as reported for other sites [92][93][94]119], and its sensitivity varied with typhoon events or images. Large variability between VIs in their responses to disturbance also took the form of varying correlation strength and direction between and within typhoon events (Table 5). However, ΔVIs based on the same spectral bands, such as ΔEVI-ΔSAVI, had stronger correlations than in other cases (e.g., ΔEVI-ΔNDII). The disparity between ΔVIs for a given typhoon suggests that cross-study comparisons of disturbance effects based on several or different VIs could be problematic.
ΔVIs shared some consistencies, nonetheless, as most VIs indicated increasing heterogeneity, and varied somehow similarly with elevation, slope, and TPI-derived classes for each typhoon ( Table  6, Table S3). In addition, both ΔNDII and ΔEVI detected increasing severity for more frequent canopy damage.
Our study has several constrains. First, although we compared the consistency among different VIs, we could not identify which VI was more accurate in detecting vegetation changes caused by typhoons due to the lack of ground truthing. Conducting ground truthing is difficult, because although the beginning and the end of typhoon seasons are relatively well understood, typhoon events themselves are unpredictable, hence it is difficult to coordinate ground surveys across the studied landscape without knowing whether a disturbance will occur and where clouds would be present on the satellite images. Although it was possible to conduct ground truthing for past disturbance events, we strongly recommend routine ground surveys on several widely spread plots for studies that aim to assess disturbance effects on vegetation whenever possible. Such ground truthing could help to identify which VI would provide most accurate detection of disturbanceinduced vegetation changes.
Second, some of our results may have been affected by cloud obstruction in analyzed images, which is substantial for all images used in this study (17.5%-48.8%). A previous study showed that cloud contamination does not distribute across the FEF, with more cloud cover at higher elevations [66]. Because typhoon-induced vegetation change varied with elevation (Table 6), the observed topographical patterns of vegetation change are likely affected to some degree by the non-random cloud contamination. Unfortunately, cloud contamination is common in the FEF and most humid forests. In fact, cloud cover prevented the inclusion of Typhoon Herb in our analyses of vegetation recovery. Third, ideally, images should be derived from the same sensor as different sensors vary in the width of their spectral bands [120] and radiometric calibrations [121,122]. However, we were constrained by the availability of high-quality (low cloud contamination) images because we studied five typhoons spanning over two decades. With the rapid advancement of remote sensing, highquality images from the same sensor should be increasingly available.

Conclusions
Comparison of the effects of five major typhoons affecting the Fushan Experimental Forest (FEF) showed substantial differences as well as some consistency in their effects. First, while typhoons all led to decreases in vegetation index (VI) values, the magnitude of change (ΔVIs) differed among events. The variability of ΔVIs among typhoons may be the product of complex interactions of their characteristics with landscape topography and the biotic conditions when the forest was disturbed (e.g., recovering from previous disturbance, soil moisture). Indeed, topography alone did not explain variation in ΔVIs among typhoons, and its explanatory power varied among the different indices. However, all typhoons shared the same positive relationship between damage severity and predisturbance vegetation condition, and all typhoons resulted in increased vegetation heterogeneity. Hence, conclusions drawn from the remote sensing studies of one typhoon may or may not stand for other typhoons in the same landscape, depending on the aspects of the effects under concern. Second, observation of greater disturbance severity for more frequently damaged cells of the FEF shows that some sites are more prone (i.e., have decreased resistance) to disturbances than others. On the other hand, the landscape generally has high resilience in order to maintain its forest cover over the many centuries of typhoon disturbance. Third, the four vegetation indices had different relationships with canopy cover damage, probably because of their different sensitivities across the light reflectance spectrum. The Normalized Difference Vegetation Index (NDVI) was, overall, less sensitive to change, which supports the findings of previous studies pointing to its saturation and overall limited sensitivity under high-biomass forest conditions. The Soil-Adjusted Vegetation Index (SAVI) and the Enhanced Vegetation Index (EVI) were highly correlated, but EVI was less related to the other two indices. Although the Normalized Difference Infrared Index (NDII) and EVI had the same change in direction after the five typhoons, they differed in the magnitude of change. Hence, we suggest using them together as complementarity Vis, because NDII and EVI are, respectively, based on short-wave infrared and near-infrared, making them sensitive to different characteristics of vegetation.
Supplementary Materials: The following are available online at www.mdpi.com/2072-4292/12/10/1654/s1. Figure S1: true-color composites of the forest cover and clouds (white) of the Fushan Experimental Forest (FEF, delimited in red) before and after each of the five studied typhoons. Table S1: Spearman's ρ for the correlation between pre-disturbance VIs (VIt0) and ΔVI associated with the typhoons. Table S2: Spearman's ρ (p-value) between corresponding ΔVIs of the five typhoons (p adjustment of Bonferroni). Table S3: β Coefficients of ordinary least squares linear models between topographical variables and changes in vegetation index values caused by Typhoons Aere, Herb, Nari, and Soudelor in relation to topographical variables: elevation, terrain slope, TPI-derived position, and aspects (converted into eight cardinal directions). Table S4: multiple comparisons on mean disturbances for different aspect-and topographical position index (TPI)-derived categories.