Monitoring Forest Infestation and Fire Disturbance in the Southern Appalachian Using a Time Series Analysis of Landsat Imagery

The southern Appalachian forests have been threatened by several large-scale disturbances, such as wildfire and infestation, which alter the forest ecosystem structures and functions. Hemlock Woolly Adelgid (Adelges tsugae Annand, HWA) is a non-native pest that causes widespread foliar damage and eventual mortality, resulting in irreversible tree decline in eastern (Tsuga canadensis) and Carolina (T. caroliniana) hemlocks throughout the eastern United States. It is important to monitor the extent and severity of these disturbances over space and time to better understand their implications in the biogeochemical cycles of forest landscapes. Using all available Landsat images, we investigate and compare the performance of Tasseled Cap Transformation (TCT)-based indices, Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Disturbance Index (DI) in capturing the spectral-temporal trajectory of both abrupt and gradual forest disturbances (e.g., fire and hemlock decline). For each Landsat pixel, the temporal trajectories of these indices were fitted into a time series model, separating the inter-annual disturbance patterns (low frequency) and seasonal phenology (high frequency) signals. We estimated the temporal dynamics of disturbances based on the residuals between the observed and predicted values of the model, investigated the performance of all the indices in capturing the hemlock decline intensity, and further validated the results with the number of individual dead hemlocks identified from high-resolution aerial images. Our results suggested that the overall performance of NDVI, followed by TCT wetness, was most accurate in detecting both the disturbance timing and hemlock decline intensity, explaining over 90% of the variability in the number of dead hemlocks. Despite the overall good performance of TCT wetness in characterizing the disturbance regime, our analysis showed that this index has some limitations in characterizing disturbances due to its recovery patterns following infestation.


Introduction
Over the last century, due to climate change and human activities the eastern North American forests have been affected by a number of large-scale disturbances, including wildfires, insect infestation, frequent droughts, and human-induced deforestation [1][2][3]. Among these, non-native pests and wildfires have been reported as two major forest disturbance agents in the southern Appalachian There are still some limitations in using a small number of before-and after-disturbance images: First, this method can limit the understanding of the disturbance progress and subsequent recovery patterns over time [56]. Second, although this approach can be a good indicator of abrupt changes such as wildfire, it may have limited capability in detecting gradual disturbances, such as gradual defoliation or tree mortality [36]. Third, finding a pair of cloud-free near anniversary images to minimize the incidental changes is often challenging [39]. With the availability of all Landsat archived images [57], many studies used multi-temporal Landsat data to detect long-term changes in the land cover and forest structures [36,[58][59][60]. These studies showed that using dense time series Landsat images is a powerful approach to characterize both abrupt and subtle disturbances [36][37][38]56,[61][62][63][64].
However, it is not clear how well the dense time series of different spectral measurements work for monitoring abrupt and gradual damage in forest ecosystems. In this study, we aim to quantify the performance of different spectral indices derived from Landsat imagery in capturing the spectral-temporal characteristics of both gradual and abrupt forest disturbances by HWA infestation and fire. Our findings provide valuable guidance for future efforts monitoring forest ecosystem disturbance using remote sensing, and the algorithm we developed here can be used as an efficient tool for forest management, planning, and restoration.

Study Site
Linville Gorge Wilderness area is located on the eastern edge of the southern Appalachian Mountains in North Carolina, USA ( Figure 1). The study area has a rugged terrain with an elevation ranging from 280 to 1300 m. This study focuses on the southern part of this wilderness area of approximately 135 km 2 . Although the majority of wilderness area remains unlogged [65], two major forest disturbances have occurred since the early 1990s. The first one is HWA infestation in the lower Linville Gorge watershed ( [66]; Figure 1), and the second is recurrent wildfires in the southern part of the lower watershed [5,65,66]. Hemlock defoliation in this study site was first reported in the early 2000s [67], and prior to this the defoliation riparian area was dominated by hemlocks in the presence of dense Rhododendron maximum shrubs in the understory [5]. Since the 1950s, this region had not been subjected to large-scale fires until November 2000 [68], when about 4000 ha of the wilderness area was burned by the Brushy Ridge human-initiated fire [5]. Two recent human-initiated fires, "Shortoff and Linville Complex" and "Pinnacle", were ignited in May 2007 and through June burned about 970 ha and 2000 ha within the gorge area, respectively [5,66]. A large portion of the Pinnacle fire overlapped with the 2000 Brushy Ridge fire [5,66]. Most recently, the White Creek fire burned the same area in March 2017. However, visual interpretation of 10 June 2016, and 3 July 2018, National Agriculture Imagery Program (NAIP) images (with a 1m resolution) showed that this fire was not as extensive as the 2000 and 2007 wildfires.

Digitizing Dead Hemlock Trees
To evaluate and compare the performance of spectral indices in capturing the disturbance intensity, the spatial distribution of the dead hemlocks was used. To delineate the dead hemlock trees, classification methods such as Maximum Likelihood were tested; however, the accuracy of these methods were not as high as expected. Therefore, to build a more accurate map for the validation process of this study, all the dead hemlock trees were visually identified from aerial images. We digitized all the dead or dying conifers exhibiting the white crown and branching structure of hemlocks (Figures 2 and 3) using a high-resolution NAIP image taken on 5 July 2010 ( Figure 1). Where there was ambiguity, we also relied on a leaf-off Color-Infrared Digital Orthophoto Quarter Quads (DOQQ) aerial image from 1998 to distinguish the evergreen hemlocks from deciduous trees in these mixed forests. A total of 16,425 dead or dying trees were manually identified and digitized in our study area (Figure 1). Note that this only captured dead overstory trees, and not those in the subcanopy or understory. While all the digitized trees are not likely to be hemlocks, the dead trees conform to known hemlocks' habitats and observed declines during this period. Ground-based observations give us confidence that the vast majority of the dead trees mapped in this area are hemlocks. We then aggregated these digitized points to the Landsat resolution (30 m) by counting all the digitized dead trees within a single pixel. In this study, the number of dead trees was used as a proxy of hemlock decline intensity at each pixel.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 20 all the digitized dead trees within a single pixel. In this study, the number of dead trees was used as a proxy of hemlock decline intensity at each pixel.   Early after mortality, the lichen-covered white branches show unique looks, often called "grey ghosts". The recovery pattern of deciduous trees can also be seen in this photo taken by Steven P. Norman, US Forest Service.  Early after mortality, the lichen-covered white branches show unique looks, often called "grey ghosts". The recovery pattern of deciduous trees can also be seen in this photo taken by Steven P. Norman, US Forest Service.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 20 Early after mortality, the lichen-covered white branches show unique looks, often called "grey ghosts". The recovery pattern of deciduous trees can also be seen in this photo taken by Steven P. Norman, US Forest Service.

Landsat Imagery
We downloaded all the available collection-1 level-2 pre-processed Landsat Thematic Mapper (TM) (n = 822) and Operational Land Imager (OLI) (n = 210) scenes from 1990 to 2018 from the USGS Earth Explorer website (https://earthexplorer.usgs.gov/). The surface reflectance products provided by level-2 Landsat TM and Landsat OLI are automatically atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [69,70] and the Land Surface Reflectance Code (LaSRC) [71] algorithms, respectively. In this study, the Landsat Enhanced Thematic Mapper Plus (ETM+) images were not used due to the gaps caused by Scan Line Corrector (SLC) failure since 2003. We masked out the cloud, cloud shadow, and snow pixels using Matlab scripts and the quality assessment (QA) bands provided by the CFMask algorithm [72,73].

Temporal Trajectory of Spectral Indices
We integrated frequent Landsat images in statistical analysis and visual interpretation to compare the performance of different spectral indices in (1) determining the temporal dynamics of major disturbances processes, including timing and duration; and (2) capturing the hemlock defoliation intensity, further cross-validated by fine-resolution imagery. We investigated the performance of eight spectral indices-including three TCT indices, two TCT ratio indices, NDVI, Enhanced Vegetation Index (EVI) [74], and DI trajectories (Table 1)-in monitoring the spectral-temporal patterns of forest disturbances. The TCTs were developed by Kauth and Thomas [49] to reduce the dimensionality of spectral data and compress the spectral variation from the original bands into a few principal components [75,76]. We used the first three principal components: brightness (TCB), greenness (TCG), and wetness (TCW). TCB was designed to capture the overall brightness, TCG represents the degree of green vegetation vigor, and TCW shows the variation in the vegetation/soil moisture [52,[75][76][77]. The TCT indices of Landsat TM and OLI scenes were calculated using the coefficient provided by Crist et al. [75] and Baig et al. [77], respectively. To reduce the topographic illumination effects of the TCT indices, we also tested the performance of "ratio" indices, such as TCB/TCG (B/G) and TCB/TCW (B/W). The primary topographic effects, changes in solar illumination angle and subsequent surface reflectance, are consistent for all the bands [55,78]; hence, we expected the TCT ratios, NDVI, and EVI to cancel out the topographic effects to some degree [55,[79][80][81]. We examined the performance of these eight indices (summarized in Table 1) in highlighting the potential temporal information of disturbance, such as complete tree morbidity, as well as in capturing the hemlock decline intensity. These indices values were rescaled to the [0, 1] range to facilitate the comparison between them.

Disturbance Detection Algorithm
To estimate the time series model of indices under stable conditions before the disturbances, we separated the study period into pre-disturbance (1990 to 2000) and post-disturbance periods (2001 to 2018). To reduce the influence of outliers, mostly caused by missed clouds, cloud shadows, or snow, the values lower than the 5th or higher than the 95th percentiles were removed from the further analyses. A time series model (Equation (1)), used by Zhu and Woodcock [39] to develop the continuous change detection classification (CCDC) algorithm, was used to separate low-frequency patterns (e.g., long-term patterns or disturbances) from high-frequency signals (e.g., seasonal phenology) at each pixel [39].
where i is the Landsat pixel number, p(t) i is the predicted index value for the Julian date t of pixel i, t is the Julian date of clear observations in pixel i, and T is the number of days in a year (365). This model has four fitted parameters: a-intercept; b and c-coefficients to model high-frequency signals; and d-first-order coefficient to model low-frequency signal. High-frequency coefficients are designed to capture seasonal phenology or intra-annual changes, while low-frequency coefficients capture long-term linear trends. An iterative nonlinear least squares method was used to fit the time series of the eight spectral indices at each pixel during the pre-infestation period using the Matlab nlinfit function (MathWorks, Natick, MA, USA). The fitted models were then extended to the post-infestation period to calculate the residuals between the observed and modeled values, being used to capture the characteristics of disturbance and recovery at each pixel.

Maximum Disturbance Timing and Magnitude
For the time series of residual values at each pixel, we applied a moving average method with one-year windows at every three-month time step. This moving average approach was applied (1) to minimize the remnant effects of seasonal variation (e.g., inter-annual variation in phenology) and (2) reduce the effect of outliers and possible short-term changes (e.g., snow or clouds) not filtered by quality control. Note that there is an observational gap between 2012 to 2013 in our Landsat data due to not using the Landsat ETM+ images, so the moving windows overlapping with this gap were not considered in our analysis. From the temporal trajectory of the mean residuals, two metrics were calculated at each pixel: (1) the year with the maximum absolute mean residual as the "maximum disturbance timing", and (2) the maximum absolute mean residual as the "maximum disturbance magnitude" as an indicator of disturbance intensity. To evaluate if the post-disturbance changes in a pixel are significantly different from the pre-disturbance conditions, a two-sample t-test was used. The t-test was applied to the residuals of the pre-disturbance period and the residuals within the 1-year window with the maximum absolute mean residual. Pixels with no significant changes from the pre-disturbance conditions were excluded for the further analysis.
Two sub-regions were selected for further evaluation. Sub-region 1, with a 5.7 km 2 forested area (80 × 80 Landsat pixels) including both HWA defoliated pixels and non-defoliated pixels, was selected to investigate the capability of the indices in capturing the infestation timing and intensity ( Figure 1). No wildfire has been reported in this sub-region. Sub-region 2 with a 47.5 km 2 area was located in the southern Linville Gorge Wilderness area where most wildfires occurred ( Figure 1). To compare the performance of our metrics in capturing the fire disturbance timing, we also used the Linville Gorge fire history map between 1980 and 2017, which has been mapped by the Wildland Fire Decision Support System (WFDSS) (https://wfdss.usgs.gov). We compare the maximum disturbance timing maps from the eight spectral indices in terms of the perimeters of fire and occurrence year of this map.
To evaluate the hemlock decline intensity estimated from the eight spectral indices, we compared the maximum disturbance magnitude maps with the number of dead hemlocks per pixel (as a representative of decline intensity) and their spatial distribution. To compare the indices' performances, the maximum disturbance magnitude maps of sub-region 1 with different dynamic ranges were rescaled using the maximum and minimum values. To evaluate the performance of the indices, we first produced boxplots of the scaled maximum disturbance magnitudes grouped by the number of dead hemlocks. Then, a linear regression analysis was applied between the number of dead hemlocks and the median maximum disturbance magnitude values of pixels with the same number of dead hemlocks.
Due to the small sample size, we grouped the pixels with 10 or more dead hemlocks into one group.

Temporal Dynamics of Indices Time Series
First, we show examples that give the spectral trajectories of TCW and NDVI in two Landsat pixels: (1) a pixel with a fire occurrence in 2007 ( Figure 4) and (2) an HWA-defoliated pixel with 12 dead hemlocks without fire ( Figure 5). The fire in 2007 was characterized as an abrupt change in the TCW (Figure 4a), NDVI (Figure 4b), EVI ( Figure S1c), and B/G ( Figure S1d). The hemlock decline, however, was captured as more a gradual change in all the indices for the selected pixels ( Figure 5 and Figure S2). Compared to the other indices, the TCB time-series is noisy with more inconsistent seasonal patterns (Figures S1a and S2a). TCB slightly increased during the disturbances (Figures S1a and S2a) because disturbed pixels might have a higher background reflectance than undisturbed ones. TCG, TCW, NDVI, and EVI values, on the other hand, decreased due to defoliation, corresponding dryness, and background color changes (Figures 4 and 5, Figures S1 and S2). Subsequently, the B/G and DI values increased with forest disturbances ( Figures S1 and S2

Temporal Dynamics of Indices Time Series
First, we show examples that give the spectral trajectories of TCW and NDVI in two Landsat pixels: (1) a pixel with a fire occurrence in 2007 ( Figure 4) and (2) an HWA-defoliated pixel with 12 dead hemlocks without fire ( Figure 5). The fire in 2007 was characterized as an abrupt change in the TCW (Figure 4a), NDVI (Figure 4b), EVI ( Figure S1c), and B/G ( Figure S1d). The hemlock decline, however, was captured as more a gradual change in all the indices for the selected pixels ( Figure 5 and Figure S2). Compared to the other indices, the TCB time-series is noisy with more inconsistent seasonal patterns (Figures S1a and S2a). TCB slightly increased during the disturbances (Figures S1a and S2a) because disturbed pixels might have a higher background reflectance than undisturbed ones. TCG, TCW, NDVI, and EVI values, on the other hand, decreased due to defoliation, corresponding dryness, and background color changes (Figures 4 and 5, Figures S1 and S2). Subsequently, the B/G and DI values increased with forest disturbances ( Figures S1 and S2), while the B/W values decreased (Figures S1e and S2e). TCW (Figures 4a and 5a) had a smaller intra-annual variability during the pre-disturbance period than the other indices.   Pixel-based observations indicate that the spectral trajectories of these indices have the potential to provide temporal information regarding disturbance onset, maximum disturbance timing, and subsequent sub-canopy recovery (Figures 4 and 5, Figures S1 and S2). The temporal patterns of the HWA-defoliated pixel showed that the gradual declines of TCG, TCW, NDVI, and EVI, as well as gradual increases of TCB, B/G, and DI began around 2003, and maximum forest disturbance occurred between 2007 and 2010 ( Figure 5 and Figure S2). This suggests that it took about 4 to 7 years for the initial HWA outbreak to complete hemlock mortality.
The maximum disturbance timing maps are shown for sub-region 1 ( Figure 6) and sub-region 2 ( Figure 7). NDVI, B/G, and TCG detected the maximum disturbance timing more accurately and consistently than the other indices, while TCB and B/W showed the weakest performance. Based on TCB and B/W, the maximum disturbance for most of the pixels in sub-region 1 occurred before 2006, which is unlikely to be based on the hemlock decline initial timing and duration of HWA defoliation [12,18,67]. Our result showed that in the maximum disturbance timing maps of NDVI and B/G, 66  Pixel-based observations indicate that the spectral trajectories of these indices have the potential to provide temporal information regarding disturbance onset, maximum disturbance timing, and subsequent sub-canopy recovery (Figures 4 and 5, Figures S1 and S2). The temporal patterns of the HWA-defoliated pixel showed that the gradual declines of TCG, TCW, NDVI, and EVI, as well as gradual increases of TCB, B/G, and DI began around 2003, and maximum forest disturbance occurred between 2007 and 2010 ( Figure 5 and Figure S2). This suggests that it took about 4 to 7 years for the initial HWA outbreak to complete hemlock mortality.
The maximum disturbance timing maps are shown for sub-region 1 ( Figure 6) and sub-region 2 ( Figure 7). NDVI, B/G, and TCG detected the maximum disturbance timing more accurately and consistently than the other indices, while TCB and B/W showed the weakest performance. Based on TCB and B/W, the maximum disturbance for most of the pixels in sub-region 1 occurred before 2006, which is unlikely to be based on the hemlock decline initial timing and duration of HWA defoliation [12,18,67]. Our result showed that in the maximum disturbance timing maps of NDVI and B/G, 66 (Table 1) within sub-region 1 with 80 × 80 Landsat pixels (Figure 1). TCW and NDVI timing maps show the best agreements with the distribution of the infested pixels; however, the TCB and B/W timing maps show the weakest performance. The legend in map (b) applies to all maps.
The disturbance timing maps of sub-region 2 from TCG, TCW, NDVI, EVI, B/G, and DI generally agree well with the WFDSS fire history map (Figure 7). The NDVI-derived disturbance timing was slightly more accurate than the other indices (Figure 7e (Figure 7). The overall performance of the brightness-related indices (TCB, B/W, and DI) in detecting this timing for pixels with a wildfire history was lower than that of the other metrics (Figure 7).  (Table 1) within sub-region 1 with 80 × 80 Landsat pixels (Figure 1). TCW and NDVI timing maps show the best agreements with the distribution of the infested pixels; however, the TCB and B/W timing maps show the weakest performance. The legend in map (b) applies to all maps.
The disturbance timing maps of sub-region 2 from TCG, TCW, NDVI, EVI, B/G, and DI generally agree well with the WFDSS fire history map (Figure 7). The NDVI-derived disturbance timing was slightly more accurate than the other indices (Figure 7e Table 1). The legend in map (b) applies to all the maps.

Hemlocks Decline Intensity and Maximum Disturbance Magnitude
Comparisons of the maximum disturbance magnitude maps with the number of dead hemlocks within sub-region 1 also demonstrate that TCW and NDVI were more in agreement with the dead hemlock distribution than the other metrics ( Figure 8). The TCB, TCG, EVI, B/G, B/W, and DI maximum disturbance magnitude maps, on the other hand, did not match with the infestation intensity as well as TCW and NDVI (Figure 8). For each index, the median maximum disturbance magnitudes for pixels with the same number of dead hemlocks were extracted from the metrics' boxplots ( Figure S3) and entered into the linear regression analysis. The relationships were statistically significant for all the indices (p < 0.01), although the accuracy varied (Figure 9). TCW and NDVI were the best predictors of the hemlock decline intensity (R 2 = 0.98 and 0.924, respectively; Figure 9), while TCB, EVI, and B/W showed the weakest performances.

Hemlocks Decline Intensity and Maximum Disturbance Magnitude
Comparisons of the maximum disturbance magnitude maps with the number of dead hemlocks within sub-region 1 also demonstrate that TCW and NDVI were more in agreement with the dead hemlock distribution than the other metrics ( Figure 8). The TCB, TCG, EVI, B/G, B/W, and DI maximum disturbance magnitude maps, on the other hand, did not match with the infestation intensity as well as TCW and NDVI (Figure 8). For each index, the median maximum disturbance magnitudes for pixels with the same number of dead hemlocks were extracted from the metrics' boxplots ( Figure S3) and entered into the linear regression analysis. The relationships were statistically significant for all the indices (p < 0.01), although the accuracy varied (Figure 9). TCW and NDVI were the best predictors of the hemlock decline intensity (R 2 = 0.98 and 0.924, respectively; Figure 9), while TCB, EVI, and B/W showed the weakest performances.  (Table 1). Grey pixels are the pixels with no significant changes from the pre-disturbance condition.  (Table 1). Grey pixels are the pixels with no significant changes from the pre-disturbance condition.  (Table 1). All the regression lines are statistically significant (p < 0.01).

Discussion
In this study, all the available Landsat images were used to explore the capability of eight spectral indices to characterize the spectral-temporal dynamics of the forest fire and infestation disturbances. Our results indicate that dense time series Landsat images enabled us to detect both the timing and intensity of the disturbances accurately. Pixel-based interpretation showed that all eight indices were able to capture the changes associated with the forest disturbances; however, their accuracy in highlighting the characteristics varied. TCB had an inconsistent seasonal cycle, and the spectral differences represented some extreme unexpected values (Figures S1a and S2a), which produced the weakest overall performance in characterizing disturbances. The inconsistent patterns of TCB and TCG are likely attributed to the sensitivity of these indices to the complex terrain [55] and the different sun angles of images with different acquisition dates [82]. Based on pixel-based visualization, the detected timing for TCB and B/G can be also explained by its postdisturbance pattern, where the TCB values during recovery stage are as high as the disturbance period ( Figures S1a and S2a). TCW also showed a smaller seasonal variability compared to TCB during the pre-disturbance period, which is also reported by another study [56]. Increasing or decreasing the pattern of the spectral indices during the disturbances corresponded well with other studies [15,42,50,75].

Maximum Disturbance Timing
The initial timing of hemlock declines from the visual interpretation of trajectories aligns with other studies, where the outbreak of HWA in North Carolina was reported in early 2000s [10,67]. TCW and NDVI, followed by TCG, B/G, EVI, and DI, have the most consistent performance in identifying the maximum disturbance timing for infestation ( Figure 6). The maximum disturbance timing estimated by the spectral indices suggests that HWA takes about 4 to 7 years to complete stand mortality. This is in agreement with other studies in this region, which documented 4-6 years for complete hemlocks mortality [12,18]. For example, Ford et al. [10] reported that 6 years after the initial outbreak of HWA in the southern Appalachian, more than 80% of the hemlocks were completely  (Table 1). All the regression lines are statistically significant (p < 0.01).

Discussion
In this study, all the available Landsat images were used to explore the capability of eight spectral indices to characterize the spectral-temporal dynamics of the forest fire and infestation disturbances. Our results indicate that dense time series Landsat images enabled us to detect both the timing and intensity of the disturbances accurately. Pixel-based interpretation showed that all eight indices were able to capture the changes associated with the forest disturbances; however, their accuracy in highlighting the characteristics varied. TCB had an inconsistent seasonal cycle, and the spectral differences represented some extreme unexpected values (Figures S1a and S2a), which produced the weakest overall performance in characterizing disturbances. The inconsistent patterns of TCB and TCG are likely attributed to the sensitivity of these indices to the complex terrain [55] and the different sun angles of images with different acquisition dates [82]. Based on pixel-based visualization, the detected timing for TCB and B/G can be also explained by its post-disturbance pattern, where the TCB values during recovery stage are as high as the disturbance period ( Figures S1a and S2a). TCW also showed a smaller seasonal variability compared to TCB during the pre-disturbance period, which is also reported by another study [56]. Increasing or decreasing the pattern of the spectral indices during the disturbances corresponded well with other studies [15,42,50,75].

Maximum Disturbance Timing
The initial timing of hemlock declines from the visual interpretation of trajectories aligns with other studies, where the outbreak of HWA in North Carolina was reported in early 2000s [10,67]. TCW and NDVI, followed by TCG, B/G, EVI, and DI, have the most consistent performance in identifying the maximum disturbance timing for infestation ( Figure 6). The maximum disturbance timing estimated by the spectral indices suggests that HWA takes about 4 to 7 years to complete stand mortality. This is in agreement with other studies in this region, which documented 4-6 years for complete hemlocks mortality [12,18]. For example, Ford et al. [10] reported that 6 years after the initial outbreak of HWA in the southern Appalachian, more than 80% of the hemlocks were completely dead. B/G, TCG, and EVI performed well in detecting the maximum disturbance timing for fire (Figure 7), which is in agreement with other studies [83,84] and may suggest that these indices perform better in characterizing abrupt disturbances.
TCW identified the maximum disturbance timing with 1-year or more lags consistently for both the infestation and fire, compared to NDVI (Figures 6 and 7). NDVI captured both the hemlocks' decline and fire disturbance timing more accurately than the other indices (Figures 6 and 7). A possible explanation could be the ability of NDVI to detect changes in vegetation density during the recovery period. Based on the pixel-level observations of the TCW trajectory, no clear recovery was observed (Figures 4 and 5). The lagged detection of the maximum disturbance timing was likely the result of TCW producing large negative values (Figures 4 and 5) even after the subcanopy recovery verified in field observations (Figure 2). The negative residuals and increased seasonal amplitudes in TCW during the recovery influenced the mean residuals and produced the delayed maximum disturbance timing, while NDVI showed steep increasing patterns at the initial stage of recovery. Another explanation for the delayed responses of all indices for the 2000 fire event (Figure 7) may be because this fire occurred late in the year (November); hence, in our moving average approach the maximum absolute mean residuals could happen with at least a one-year delay (Figure 7). In addition, our analysis did not detect the 2000 timing for many pixels within the perimeter, which might be explained by the fact that the WFDSS map includes many lightly burned or even unburned areas. For example, the 2010 NAIP image still showed that some hemlocks stands within the 2000 fire perimeter from the WFDSS map, and it was very likely that the fire in 2000 lightly burned the study area.

Maximum Disturbance Magnitude
Although TCG was a good indicator of disturbance timing for both the infestation and fire (Figures 6c and 7c), it was not a good predictor of the hemlock decline intensity (Figure 8c). It was likely that the biased values in the TCG trajectory affected the maximum disturbance magnitudes. Both the B/G and B/W ratios, designed to minimize the possible topographic effects, performed better than TCB in capturing the disturbance timing and magnitude. However, these ratio indices were sensitive to extremely small values in the denominator, and therefore produced quite unstable time series. Although the maximum disturbance timing maps derived from B/G were reliable (Figures 6g and 7g), the maximum disturbance magnitude did not match well with the number of dead hemlocks (Figure 8g). The Disturbance Index (DI), introduced to accentuate the difference between the disturbed and undisturbed forest conditions [42,85], did not perform better than TCW and NDVI in detecting the maximum disturbance timing and hemlock decline intensity (Figures 6i, 7i and 8i). The better performance of DI in detecting infestation intensity, compared to TCB and TCG, might be explained by the usage of the normalized TCT values in its calculation, which would minimize the effects of seasonality and variability due to bidirectional reflectance [85,86].
The maximum disturbance magnitudes of TCW and NDVI were relatively comparable with the spatial distribution of the dead hemlock trees (Figure 8). In all the disturbance magnitude maps (Figure 8), some pixels with no digitized dead hemlocks had large absolute disturbance magnitudes, while some severely infested pixels had small absolute disturbance magnitudes. These unexpected results might be explained in a few ways. First, the complex rugged terrain of the study area can limit the performance of spectral indices, especially the ones that are more sensitive to topography variation, such as TCB and TCG [55,87], leading to biased high or low relative residuals and inaccurate estimates. Song and Woodcock [55] also reported that while TCB and TCG were sensitive to topographic effects especially when the solar zenith angle was large, TCW and NDVI are less sensitive. While DI could minimize the topographic effect by using normalized TCT, it can still be affected by topography, inherited from the TCB and TCG variability. Studies also reported that although EVI can reduce the effects of atmospheric conditions and vegetation saturation, compared to NDVI it is sensitive to topographic and viewing illumination effects due to the presence of the soil adjustment factor in its equation [80,[88][89][90]. Second, there could be uncertainty in identifying the hemlock decline intensity at each pixel based on the number of dead hemlocks. HWA kills hemlocks across all ages and sizes eventually [7], but their progresses are different based on their sizes [25]. Therefore, the number of overstory dead hemlocks, as a proxy of hemlock decline intensity, without considering the projected crown sizes, overlapping, and possible understory dead hemlocks, may either underestimate or overestimate the infestation intensity within a single pixel, not even mentioning geometric errors. Third, requiring clear remote sensing observations can produce small sample sizes in pixels with persistent clouds or cloud shadows, thereby undermining the reliability and robustness of parameter estimates in our regression models [39].

Post-Infestation Recovery
Another important information that can be retrieved is the recovery progress after the disturbances. For example, the spectral trajectories of two example pixels (Figures 4 and 5, Figures S1 and S2) indicated that all indices, except TCB and B/W, indicated onsets of recovery after 7-8 years of the initial infestation ( Figure 5), which also showed the stable seasonal patterns around 2014. While TCW reliably identified the timing and magnitude of the disturbance, it has limitations in determining the potential post-infestation recovery. The TCW residuals decreased toward negative values during the entire period of infestation and recovery, and they never recovered, especially during the dormancy seasons, to the pre-infestation level ( Figure 5). The current evergreen hemlocks stands are being replaced by deciduous trees, such as black birch (Betula lenta L.) and maple (Acer spp.) (Figure 2) [91]. Emerging leaf-off periods during this transition from evergreen to deciduous trees might indicate that TCW might have limitations in detecting forest decline and following recoveries using image difference [15]. Overall, our time series analysis also suggested that the seasonal dynamic of some indices could be appropriate to characterize forest recovery following the disturbances. However, such indices should be chosen more carefully, and the interpretation of the recovery patterns should depend on subsequent species replacement. Future studies might be needed to validate the results with ground-based observations in terms of recovery patterns.

Limitations of the Developed Method
Our method has some limitations. First, to fit the time series model during an undisturbed period, it is necessary to have detailed historical information of the study area, such as the fire history and perimeter, as well as high-resolution aerial photos. Therefore, in a region with frequent disturbances, such as wildfires, defining a pre-disturbance period can be challenging. Second, in this study to evaluate the performance of the indices we developed our algorithm to focus only on the maximum disturbance events over the period of recording. However, it should be noted that this method could be adapted to identify (local) maxima when the analysis of multiple disturbances is of interest. Our analysis on pixels with both the 2000 and 2007 fires in the study site ( Figure S4) showed that our approach could be modified to detect multiple disturbances if we were to dig into several disturbance timings.

Conclusions
In this study, we examined the capability of eight spectral indices to capture the spectral and temporal dynamics of fire and infestation disturbances in the southern Appalachian forest using all available Landsat images. We estimated these disturbance characteristics in space and time based on the residuals between the observed and predicted values of multi-frequency remote sensing data, separating the inter-annual disturbance patterns (low frequency) and seasonal phenology (high frequency) signals. In addition, we investigated the performance of these indices in capturing the hemlock decline intensity, validated with the number of dead trees identified at each pixel. Our results suggested that NDVI, followed by TCW, showed the most reliable performance in detecting the temporal-spectral properties of forest disturbances in complex terrain, while TCB, B/W, and EVI have the poorest performance in characterizing disturbances. Although TCW showed relatively consistent results in capturing the disturbance regimes, it might have limitations in detecting recovery patterns.
Our study highlights the value of using long-term Landsat imagery to monitor forest conditions over space and time for forest disturbances and recovery.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/15/2412/s1, Figure S1: The spectral trajectory of TCB, TCG, EVI, B/G, B/W, and DI; their residuals; and the mean values within 1-year moving windows with 3-month steps at a Landsat pixel with a fire in 2007. Figure S2: The spectral trajectory of TCB, TCG, EVI, B/G, B/W, and DI; their residuals; and the mean values within 1-year moving windows with 3-month steps at a Landsat pixel with 12 dead hemlocks. Figure S3: Boxplots of indices' maximum disturbance magnitude as a function of the number of dead hemlocks in a pixel. Figure S4: