Damage Analysis of Three Long-Track Tornadoes Using High-Resolution Satellite Imagery

: Remote sensing of tornado damage can provide valuable observations for post-event surveys and reconstructions. The tornadoes of 3 March 2019 in the southeastern United States are an ideal opportunity to relate high-resolution satellite imagery of damage with estimated wind speeds from post-event surveys, as well as with the Rankine vortex tornado wind ﬁeld model. Of the spectral metrics tested, the strongest correlations with survey-estimated wind speeds are found using a Normalized Di ﬀ erence Vegetation Index (NDVI, used as a proxy for vegetation health) di ﬀ erence image and a principal components analysis emphasizing di ﬀ erences in red and blue band reﬂectance. NDVI-di ﬀ erenced values across the width of the EF-4 Beauregard-Smiths Station, Alabama, tornado path resemble the pattern of maximum ground-relative wind speeds across the width of the Rankine vortex model. Maximum damage sampled using these techniques occurred within 130 m of the tornado vortex center. The ﬁndings presented herein establish the utility of widely accessible Sentinel imagery, which is shown to have su ﬃ cient spatial resolution to make inferences about the intensity and dynamics of violent tornadoes occurring in vegetated areas.


Introduction
Short-term meteorological events such as severe thunderstorms and tornadoes can leave lasting effects on the natural landscape. Wind, hail, and tornado damage to forests and agricultural regions can be severe enough to identify with both ground level observations [1][2][3][4][5] and satellite imagery [6,7]. Aerial imagery of tornado damage has been used since the mid-20th century and allows for a "birds-eye" view to identify patterns that may not be discernible at ground-level, or would otherwise be unidentified [8,9]. This approach can be particularly useful in rural or remote regions that are difficult to survey after the event because they are inaccessible at ground level [10][11][12][13]. Additionally, tornado intensity is difficult to infer in areas where damage indicators are lacking [10]. Recent advances in moderate-and high-resolution multispectral imagery have allowed for further study of vegetation damage by tornadoes [12], as well as tornado characteristics that may not be apparent at the ground level, including the extent of tornado path and effects of local topography [14][15][16]. These techniques are most effective when using imagery with a high spatial resolution in regions of strong damage and homogenous land cover [17,18]. Normalized Difference Vegetation Index (NDVI, used as a proxy for vegetation health) differences and principal components techniques can be used to discern damage as weak as (E)F-1 [14,19,20] on the Enhanced Fujita (EF) scale, with more pronounced decreases in NDVI corresponding to greater damage in stronger tornadoes [14,21]. NDVI difference calculated from an image obtained before the tornado event and an image obtained after the event, as is performed in this study, is an effective tool in damage identification [14]. However, debris removal and other non-meteorological processes occurring in the weeks after the event can also affect NDVI in ways that

Data Description
The study region of southeast Alabama and southwest Georgia was well sampled by Sentinel-2 MSI on 24 February 2019, seven days before the tornadoes, and 6 March 2019, three days after the tornadoes, producing > 95% cloud-free images on both dates. The use of imagery from only three days after the event limits the effects of long-term cleanup, post-event vegetation regrowth, etc., on spectral differences between the two images. A majority of the landcover in the region is forested, and the tornadoes occurred early in the growing season. Thus, this event represents ideal conditions for a unique opportunity to study and compare remotely sensed damage patterns between these three strong tornadoes using high-resolution imagery with minimal cloud contamination. Images representing Level-1C geometrically and radiometrically corrected top-of-atmosphere reflectance in the near-infrared and visible red, green, and blue bands were downloaded from the United States Geological Survey (USGS) EarthExplorer website [33]. These bands have central wavelengths at 490 (blue), 560 (green), 665 (red), and 842 (near-infrared) nanometers, each with a spatial resolution of 10 m. Level-1C Sentinel products are registered to subpixel accuracy, and the images used in this study were projected into North American Datum 1983 Universal Transverse Mercator Zone 16N, as appropriate for the longitude of this study region. Using a geographic information system (GIS) platform, each of the images were clipped to a spatial extent representing a three-kilometer buffer around the paths of each tornado, as defined by the surveys conducted by local National Weather Service (NWS) meteorologists. The three-kilometer width of the buffer is about twice the maximum width of the Beauregard-Smiths Station tornado, the widest of the three tornadoes studied here, and allows for spectral patterns of damaged areas to be interpreted in the spatial context of undamaged areas. The landcovers within the tornado damage paths were identified using data from the National Land Cover Database 2016 (NLCD 2016) database. NLCD 2016 is a robust and commonly used data source that identifies land covers across the United States at a 30-m spatial resolution [34].
Data from the NWS Birmingham post-event surveys included 493 damage indicator point locations and estimated wind speeds at those locations. These surveys were performed from 4 March to 7 March 2019, the days immediately following the tornado event. A wide variety of damage indicators were surveyed, including both built and natural features such as trees, agricultural buildings, manufactured homes, and single family residences [31]. Detailed descriptions of the NWS survey process for strong tornadoes can be found in Edwards et al. (2013) [10] and Burgess et al. (2014) [35]. Estimated wind speeds at damage indicators in the NWS Peachtree City survey were not available, so these point locations were not used in this study.

Normalized Difference Vegetation Index
NDVI was calculated with the red and near-infrared reflectance images from both 24 February and 6 March using Equation (1) in which Re f NIR is the reflectivity of a given pixel in the near-infrared band, and Re f Red is the reflectivity of that pixel in the red band. Pixels containing photosynthetically active vegetation are more reflective in near-infrared than in red wavelengths, and thus NDVI is commonly used to measure vegetation health and activity [36].

Principal Components Analysis
The second technique used was a principal component analysis (PCA). PCA is an image enhancement technique that involves an orthogonal transformation of reflectance data into a new set of variables, called principal components, which capture as much variability in the data as possible. While NDVI-differencing uses information from a total of four bands-two bands (red and near-infrared) each on two separate dates (24 February and 6 March)-PCA allows for the addition of more bands into the analysis. The blue and green bands from 24 February and 6 March were included in the PCA, for a total of eight bands. Since PCA is sensitive to extreme values and distributions with uneven scaling, the images were centered and standardized. Pixels with values equal to the mean reflectance in each image were assigned a value of 0.0. Values one standard deviation above and below the mean were assigned the 1.0 and -1.0, respectively, and so on. All eight bands were centered and standardized in this manner before performing the PCA, which produced eight principal component images. This standardization technique has been observed to enhance signal to noise ratio [37] and is used in change detection research to improve results (e.g., [38]).

Correlating Remotely Sensed Products and Estimated Wind Speeds
To determine the correlations between survey-estimated wind speeds and multispectral indices such as NDVI and PCA, the values of the pixels containing damage indicator locations were obtained from the 6 March NDVI image, the NDVI-difference image (6 March minus 24 February), and the eight principal component images. Then, these pixel values were correlated with the estimated wind speed at that location from the surveys using a simple linear regression for each product.
The two products that had the strongest correlations with estimated wind speeds were NDVI-difference and PC7. From these products, new images were created by applying a square-shaped smoothing window to each pixel in both the NDVI-difference and PC7 image, in which the pixel at the center of the window was reassigned a value equal to the mean of the original values of every pixel in the window. This process was performed using a three-pixel-by-three-pixel (30 m by 30 m) window on the NDVI-difference and PC7 images, and then again using a five-pixel-by-five-pixel (50 m by 50 m) window. This procedure enhanced the correlations between pixel value and estimated wind speed in all four of the new images by reducing statistical noise.

Evaluating Remotely Sensed Products Across Tornado Widths
The patterns in NDVI difference and PC7 values across the width of each tornado were evaluated by randomly placing 1000 points within the paths of the three tornadoes listed in Table 1. This number of pixels was chosen to ensure sampling throughout the length and width of each tornado path without necessitating the selection of pixels nearby to one another, whose values may be spatially autocorrelated, especially after smoothing windows were applied. The tornado paths were defined by the polygon shapefiles in the NWS Birmingham and NWS Peachtree City storm survey GIS data. The only alterations made to these polygons were to ensure continuity in tornado path width at the Chattahoochee River, which is the boundary between the NWS Birmingham and NWS Peachtree City county warning areas. Values of pixels containing each of the 2000 points placed within the paths of the Davisville and Weedon Field tornadoes were extracted using the NDVI difference and PC7 images created through application of the three-pixel-by-three-pixel smoothing window. The points placed within these two tornado paths were a minimum of 55 m apart to ensure that no two nearby points fell in pixels with values calculated from overlapping smoothing windows. The ratio of sampled pixels to total pixels within the surveyed tornado path was 1:403 for the Davisville tornado and 1:243 for the Weedon Field tornado. This process was repeated in a similar manner for the Beauregard-Smiths Station tornado. Since the Beauregard-Smiths Station tornado path had a much larger area than the other two tornado paths, a larger minimum-distance threshold of 85 m was used for random point placement. This larger minimum-distance threshold allowed for values to be sampled from the Atmosphere 2020, 11, 613 6 of 20 NDVI difference and PC7 images created with a five-pixel-by-five-pixel smoothing window, while still ensuring that no two sampled values were calculated from overlapping smoothing windows. The ratio of sampled pixels to total pixels in the Beauregard-Smiths Station tornado path was 1:873.
For each of the three tornadoes, any value from a point that fell in a pixel with a negative NDVI value was excluded from the further analysis, since these pixels were likely water bodies and thus would be inappropriate for estimating tornado wind speeds. This totaled five points in the Beauregard-Smiths Station tornado path, six points in the Davisville tornado path, and fifty-one points in the Weedon Field tornado path, which crossed the Chattahoochee at a wide point in the river. The remaining NDVI-difference and PC7 values were plotted against the distance from the corresponding point to the center of the tornado path. The distances of the points to the left of the tornado path center, relative to storm motion, were multiplied by −1, so that when plotted on a traditional cartesian grid, these points would be separated from points on the right of the path center. The relationship between NDVI difference and point distance from tornado path center, as well as the relationship between PC7 and point distance from tornado path center, were evaluated by fitting a smoothing spline to each of the six plots (NDVI difference and PC7 for each of the three tornadoes). Smoothing splines were used for this purpose because of their flexibility in representing noisy data in which the relationship between dependent and independent variables is not consistent across the domain of independent variables. Since the relationship between wind speed and distance from tornado path center is complex even in basic tornadic wind field models, more simple regressions would not be flexible enough for this purpose. The shapes of these smoothing splines were then compared to the pattern of maximum ground-relative wind speeds across the lateral width of a tornado modelled by a Rankine vortex, such as that shown in Figure 2.
the relationship between PC7 and point distance from tornado path center, were evaluated by fitting a smoothing spline to each of the six plots (NDVI difference and PC7 for each of the three tornadoes).
Smoothing splines were used for this purpose because of their flexibility in representing noisy data in which the relationship between dependent and independent variables is not consistent across the domain of independent variables. Since the relationship between wind speed and distance from tornado path center is complex even in basic tornadic wind field models, more simple regressions would not be flexible enough for this purpose. The shapes of these smoothing splines were then compared to the pattern of maximum ground-relative wind speeds across the lateral width of a tornado modelled by a Rankine vortex, such as that shown in Figure 2. Figure 2 approximates how maximum ground-relative wind speeds vary across the lateral width of a tornado, perpendicular to its translational motion, in a simplistic Rankine vortex model. This figure is used as a template for visual comparison to the distribution of remotely sensed values within the paths of the three tornadoes listed in Table 1. The values in Figure 2 were calculated using the following equations, adapted from those used in Holland et al. (2006), Beck and Dotzek (2010), and Karstens et al. (2013) [26][27][28] and translated to cartesian coordinates in which is the maximum tangential wind velocity in the vortex, is the radius from vortex center at which occurs, is the maximum radial wind velocity towards the center of the vortex, is the translational velocity of the tornado vortex, is the distance from tornado center measured on a plane oriented perpendicular to its translational vector, and is the approximate maximum ground-relative wind speed at distance . Equation (2) is used when | | < , and Equation (3)   To examine the effects of landcover, the 1000 pixels randomly selected from each of the three tornado paths were placed into four different landcover groups depending on the landcover in which V t is the maximum tangential wind velocity in the vortex, R max is the radius from vortex center at which V t occurs, V r is the maximum radial wind velocity towards the center of the vortex, V s is the translational velocity of the tornado vortex, x is the distance from tornado center measured on a plane oriented perpendicular to its translational vector, and V gr is the approximate maximum ground-relative wind speed at distance x. Equation (2) is used when |x| < R max , and Equation (3) is To examine the effects of landcover, the 1000 pixels randomly selected from each of the three tornado paths were placed into four different landcover groups depending on the landcover identified at their locations in NLCD 2016. The first group was comprised of forested landcovers, including deciduous forest (class 41), evergreen forest (class 42), mixed forest (class 43), and woody wetlands (class 90). The second group consisted of agricultural landcovers, including pasture/hay (class 81) and cultivated crops (class 82). The third group was comprised of shrubs and grassland, including shrub/scrub (class 52), and grassland/herbaceous (class 71). The fourth group consisted of urban landcovers, including developed open space (class 21), developed low intensity (class 22), developed medium intensity (class 23), and developed high intensity (class 24). Of the 3000 randomly selected pixels, 87 were identified by the NLCD 2016 as other landcover types that were either unlikely to exhibit spectral changes due to tornado damage, including open water (class 11) or barren rock (class 31), or were otherwise not common enough in the study area to draw meaningful conclusions.

Principal Component Analysis of Pre-and Post-Event Imagery
A PCA was performed on eight centered and standardized bands: blue, green, red, and near infrared obtained pre-event on 24 February 2019, and the same four bands obtained post-event on 6 March 2019. The eigenvalues and eigenvectors produced in this analysis are shown in Table 2. Each of the eight principal component images were then correlated with estimated wind speeds at the locations of 493 damage indicators from the NWS Birmingham storm survey. Several of the principal components, such as PC3 and PC6, emphasized differences in reflectance on the two dates, with eigenvalues of opposite signs for the same band on each date.

Correlating Remotely Sensed Products and Estimated Wind Speeds
The 6 March NDVI, NDVI-difference, and eight PCA values from non-water pixels containing damage indicators from the NWS Birmingham storm survey were then correlated to the estimated wind speeds at those locations to determine which two of these images were the best proxies for damaging tornadic wind speeds. The slopes, intercepts, and R 2 values for each of these relationships are listed below in Table 3. This table shows that the NDVI-difference and PC7 images exhibited the strongest correlations with estimated wind speeds. Negative NDVI differences within the tornado paths are likely due to vegetation damage, the degree of which is strongly dependent on wind speed. The negative slope in Table 3 indicates that greater vegetation damage is associated with stronger estimated wind speeds. PC7 values also exhibit a relatively strong correlation to estimated wind speeds in the NWS Birmingham survey compared to all other PC values. Smoothing filters were Atmosphere 2020, 11, 613 8 of 20 applied to the NDVI-difference and PC7 images, as noted in the Methods section, which increased the correlation strengths by eliminating some statistical noise. These correlations are shown in Table 3, as well. Smoothed NDVI-difference and PC7 images show stronger correlations with estimated wind speeds, as shown by the R 2 values that are higher for these images than the others listed in Table 3. Based on these results, the smoothed NDVI-difference and PC7 values were evaluated across the width of the tornado paths. Visually, remotely sensed signatures consistent with swaths of tornado damage were apparent in many locations within the extent of NWS survey-defined tornado paths. Three examples using post-event NDVI, NDVI-difference, and PC7, the unsmoothed images most strongly correlated with estimated wind speed in Table 2, are shown in Figure 3. Figure 3a shows 6 March (post-event) NDVI in portions of both the Davisville and Beauregard-Smiths Station tornado paths. The NDVI imagery exhibits patches of both high NDVI in dark green and low NDVI in light green, indicating high and low photosynthetic activity, respectively. The heterogenous landcover in this portion of the study area, a mix between forest and agricultural fields, gives the image its "patchy" look. A faint, elongated region of low NDVI values is visible within the Beauregard-Smiths Station tornado path, in the upper right quadrant of Figure 3a, evidence of reduced vegetative activity due to tornado damage. Figure 3b shows NDVI difference in a portion of intense damage in the Beauregard-Smiths Station tornado path, near the Chattahoochee River. Decreases in NDVI between 24 February and 6 March are evident as dark brown within the NWS survey-defined path of the tornado. Figure 3c shows PC7 values along the Weedon Field tornado path. High PC7 values, correlated with stronger estimated wind speeds in Table 3, are colored in dark brown and exist both inside and outside of the tornado path. This is likely a region of unidentified tornado damage, as will be discussed further. Elongated regions of lighter colored pixels are regions of riparian vegetation, which are detected by PC7 as well. Atmosphere 2020, 11, x FOR PEER REVIEW 9 of 20

Examination of Damage Signatures Across Tornado Path Width
Figures 4-6 show how NDVI-difference and PC7 values vary across the lateral width of each of the three tornado paths. Since NDVI decreases are evidence of stronger wind speeds, the NDVIdifference axes have been reversed for easier comparison to Figure 2. Figure 4 shows the NDVI-difference and PC7 values of pixels containing points placed randomly within the Beauregard-Smiths Station tornado path, plotted in blue. The red lines are smoothing splines fitted with six degrees of freedom. Splines with degrees of freedom varying from two to ten were also fitted to these plots, but the spline with six degrees of freedom for Figure 4a most strongly resembled the curve shown in Figure 2, hence why it is shown here. The spline in Figure 4a, showing NDVI difference, exhibits an absolute minimum at x = 140 m, as well as a relative minimum at x = −120 m (note that the NDVI-difference values on the y-axis have been reversed). The spline shows rapidly increasing NDVI differences at x-values outside of this range. This suggests that the strongest vegetation damage, and inferred strongest wind speeds, occurred within 130 ± 10 m of the center of the tornado path. Additionally, the shape of the spline in Figure 4a bears a resemblance to the shape of maximum ground-relative wind speeds in Figure 2, which shows quickly decreasing wind speeds outside a radius of (1.0 in Figure 2). The smoothing spline in Figure 4b does not exhibit  this shape. Although there is a relative maximum in PC7 at about x = 140 m, the shape of the smoothing spline is relatively flat, or even increasing, at values of |x| > 200 m, unlike the shape evident in Figure 2. This suggests that PC7, unlike NDVI difference, is ineffective in this tornado path as a proxy for the wind field pattern in the Rankine vortex, despite the relatively strong correlation with estimated wind speeds at damage indicators in the NWS Birmingham survey (Table 3).    Table 4 shows the NDVI difference and PC7 values of the randomly selected pixels in each of the three tornado paths, categorized into four groups of similar landcover classes in NLCD 2016, as  Figure 4 shows the NDVI-difference and PC7 values of pixels containing points placed randomly within the Beauregard-Smiths Station tornado path, plotted in blue. The red lines are smoothing splines fitted with six degrees of freedom. Splines with degrees of freedom varying from two to ten were also fitted to these plots, but the spline with six degrees of freedom for Figure 4a most strongly resembled the curve shown in Figure 2, hence why it is shown here. The spline in Figure 4a, showing NDVI difference, exhibits an absolute minimum at x = 140 m, as well as a relative minimum at x = −120 m (note that the NDVI-difference values on the y-axis have been reversed). The spline shows rapidly increasing NDVI differences at x-values outside of this range. This suggests that the strongest vegetation damage, and inferred strongest wind speeds, occurred within 130 ± 10 m of the center of the tornado path. Additionally, the shape of the spline in Figure 4a bears a resemblance to the shape of maximum ground-relative wind speeds in Figure 2, which shows quickly decreasing wind speeds outside a radius of R max (1.0 in Figure 2). The smoothing spline in Figure 4b does not exhibit this shape. Although there is a relative maximum in PC7 at about x = 140 m, the shape of the smoothing spline is relatively flat, or even increasing, at values of |x| > 200 m, unlike the shape evident in Figure 2. This suggests that PC7, unlike NDVI difference, is ineffective in this tornado path as a proxy for the wind field pattern in the Rankine vortex, despite the relatively strong correlation with estimated wind speeds at damage indicators in the NWS Birmingham survey (Table 3). Figure 5 shows the same NDVI-difference and PC7 plots for randomly placed points within the path of the Davisville tornado. The smoothing splines in Figure 5 show little evidence of strong patterns in NDVI difference and PC7 across the width of the Davisville tornado path. These splines have six degrees of freedom for consistency with the smoothing splines in Figure 4, but splines with varying degrees of freedom were fitted to the data and none of them resembled the curve of Figure 2. There is a weak minimum in the Figure 5a spline at x = 100, as well as a weak maximum in Figure 5b at x = 140. These are shown as dashed vertical lines in Figure 5. While the strongest ground-relative wind speeds, and thus the greatest damages, usually occur on the right side of a cyclonically rotating tornado, the smoothing splines in Figure 5 are relatively flat compared to the curve of maximum wind speeds in Figure 2. The remotely sensed damage signature in the Davisville tornado path does not match well with the Rankine vortex model of maximum ground-relative wind speeds using either NDVI difference or PC7. One explanation for this lack of resemblance is that the Davisville tornado was only rated EF2 with peak wind speeds of 185 km/h, which is weaker than the Beauregard Smiths-Station tornado's peak estimated wind speeds of 284 km/h.

Landcover
The plots for pixels in the Weedon Field tornado path are shown in Figure 6. The smoothing splines in Figure 6 are relatively flat across the width of the Weedon Field tornado path. Once again, these splines have six degrees of freedom for consistency with Figures 4 and 5, but no spline from these data resembled the shape of the curve in Figure 2, even when the degrees of freedom were adjusted. A weak minimum in NDVI-difference values is evident near the center of the tornado path in Figure 6a, as are marginal increases in NDVI difference and decreases in PC7 (Figure 6b) at distances of |x| > 200 m. However, as in Figure 5, no notable variations in damage in Figure 6 are comparable to the wind field in Figure 2. Table 4 shows the NDVI difference and PC7 values of the randomly selected pixels in each of the three tornado paths, categorized into four groups of similar landcover classes in NLCD 2016, as described in the Methods section. Mean NDVI-difference values tend to be most negative in forested and urban landcover groups. PC7 values are highest for these landcover groups as well. The shrub/grassland and agricultural landcover groups exhibit the highest NDVI-difference and lowest PC7 values. However, variance is relatively high, and the differences between the forested and agricultural landcover groups are greater in magnitude than the standard deviations of either of these groups.

A region of Displacement between the Survey-Defined Tornado Path and NDVI Decrease
While the tornado-path polygons produced in the post-event surveys [30,31] matched the remotely sensed damage areas from this study quite well, one region of apparent displacement was located in southern Stewart County and northern Quitman County, Georgia, shown in Figure 7. The brown pattern below (to the south of) the survey-defined tornado path in Figure 7 is an area of NDVI decrease, indicating vegetation damage, between 24 February and 6 March. This region runs parallel one kilometer to the south of the survey-defined tornado path for a distance of about five kilometers along it, or about one-tenth the length of the Weedon Field tornado path. This was the only portion of the three tornado paths that exhibited an elongated spatial displacement between apparent damage and the path polygons themselves, so displacement between the remotely sensed damage signatures and survey-defined paths is unlikely to substantially compromise the results shown in Figures 4-6.   Table 3 show that NDVI differencing and PCA are effective in identifying areas of tornado damage and approximating damage intensity. However, neither product is without its  Table 3 show that NDVI differencing and PCA are effective in identifying areas of tornado damage and approximating damage intensity. However, neither product is without its drawbacks. NDVI differencing is less effective outside the growing season or in regions with little vegetation, as several studies have noted [17][18][19]21]. The relative lack of damage indicators for vegetation [10] may mean that the intensities of these tornadoes were underestimated in forested areas, where NDVI decreases are most pronounced (Table 4). A comparison of NDVI difference and EF-scale damage estimates in forested areas (e.g., [11]) would shed light on this problem. The interpretation of PCA varies by study area, band inclusion, and other factors, and is thus difficult to compare from one event to the other, despite favorable results in this study and others [14,19,20].

Results in
The correlations of both PCA and NDVI-difference values with estimated wind speeds were improved by a smoothing window, which served to dampen noise in the reflectance data. While high-resolution imagery provides more detail on damage patterns than imagery from lower-resolution sensors [14], the positive results obtained through use of the smoothing window suggest that substantial variability in the reflectance of neighboring pixels in high-resolution imagery presents a challenge in estimating wind speeds over small spatial scales. A hypothetical example of this would be in a low-density residential region. A pixel comprised of reflectance from a non-vegetated surface like a house would be unlikely to exhibit a substantial decrease in NDVI due to tornado damage. However, neighboring pixels comprised of vegetated surfaces such as lawns and suburban tree cover would exhibit stronger decreases in NDVI as this vegetation is damaged by the tornado or covered in debris. The pixel containing reflectance from the house would have its NDVI-difference value decreased by a smoothing window that includes NDVI decreases in neighboring vegetated pixels, more accurately reflecting strong wind speeds.
Another notable finding in this study is the high eigenvector magnitude on the blue bands in PC7 (Table 2) and the relatively strong correlation between PC7 and estimated tornadic wind speeds (Table 3) compared to other PCs. For the 24 February image, blue reflectance has a positive PC7 eigenvector, and red reflectance has a negative eigenvector. However, for the 6 March image, blue reflectance has a negative eigenvector, while red reflectance has a positive eigenvector. This means that pixels exhibiting a higher ratio of red reflectance to blue reflectance on 6 March than on 24 February will have high PC7 values, which are correlated to stronger wind speeds in Table 3. For example, a pixel dominated by vegetation in the 24 February image will have a nearly even ratio of red to blue reflectance because vegetation reflects at a similar rate in these wavelengths. If vegetation in that pixel is damaged by a tornado, exposing soil beneath, the ratio of red to blue reflectance will increase because soil is more reflective in red wavelengths than in blue wavelengths. As a result, this pixel would register a high PC7 value given the eigenvectors in Table 2, indicating strong winds at that location (Table 3). Some previous studies that used PCA to delineate tornado paths did not include blue reflectance in the analysis, either because the sensor of choice did not obtain blue reflectance data [19], or for the sake of consistency with previous research [14]. While these are logical reasons for the omission of blue reflectance, including it in either PCA or other multispectral indices may be beneficial to future studies like this one.
The similarity between the smoothing spline of NDVI difference across the width of the Beauregard-Smiths Station tornado path (Figure 4a) to the pattern of maximum ground-relative wind speeds across the width of a tornado in a Rankine vortex model ( Figure 2) is another intriguing finding of this study. NDVI-difference values reached minima, indicating strong vegetation damage, at distances of 130 ± 10 m on both sides of the center of the tornado path in Figure 6a. This distance may serve as a starting point for an estimation of R max , although it should be noted that direct estimates of R max using radar data can differ from the width of maximum on-the-ground damage [39]. Empirical calculations of V t and V r in Equations (2) and (3) are difficult, requiring the use of high-resolution tree fall data [26][27][28]40], and are not likely feasible using NDVI-difference products at a 10-m resolution. Surprisingly, a signature similar to Figure 2 was not apparent in the PC7 values in the Beauregard-Smiths Station tornado path (Figure 4b), despite their relatively strong correlation to estimated wind speeds in the NWS Birmingham survey (Table 2). A maximum in the smoothing spline of PC7 values, indicating stronger wind speeds, was located at 140 m to the right of the tornado path center, but a strong decrease in PC7 values was not evident in the outer flanks of the path. The reasons for this are unclear, but perhaps small-scale variations in tornadic wind fields are best proxied by indices that detect vegetation damage, such as NDVI, rather than PCA techniques that might be affected by other unrelated factors. Yuan et al. (2002) suggested that NDVI differencing outperforms PCA in detecting areas of F2 (the Enhanced Fujita scale was not yet implemented at the time the study was published) and F1 damage [19]. It is possible that the NDVI-differencing technique replicates damage patterns resembling the Rankine vortex wind field in Figure 2 because it is more sensitive specifically to relatively weak damage, which is likely to be found near the exterior of the tornado path. If this is the case, detecting a damage signature consistent with the wind field in Figure 2 requires both intense tornado damage and a spectral index that is sufficiently sensitive to weaker damage as well.
The Davisville and Weedon Field tornadoes did not exhibit patterns in NDVI differences or PC7 values consistent with the Rankine vortex wind field in Figure 2. This is probably attributable to their weaker intensities (Table 1). Even in violent tornadoes, the strongest wind speeds usually encompass only a fraction of the tornado path [35,41,42], and regions of greatest vegetation damage are often confined to small areas [15]. NDVI differencing becomes much less effective at detecting tornado damage at < EF2 intensities [14,19,21]. The Davisville and Weedon Field tornadoes likely did not produce sufficiently intense damage over large enough areas for the randomly sampled pixels in Figures 4 and 5 to show substantial variation in NDVI difference and PC7 across the width of the tornado paths.
Variations in landcover within the tornado paths may also play a role in masking the relationship between multispectral indices and damage. While a detailed examination of damage signatures in various landcovers is outside the scope of this study, Table 4 suggests that damage tends to be more evident in regions of forested and urban landcovers. Forested areas are comprised of dense vegetation, and vegetation indices such as NDVI are particularly effective for identifying vegetation damage [17]. Of the 118 pixels in urban areas (Table 4), 113 (96 percent) were classified in NLCD 2016 as either "Developed, open space" or "Developed, low intensity", which are both classes in which lawn grasses and trees may comprise a substantial portion of cover. As such, NDVI and other vegetation indices can still be effective in detecting damage to this vegetation. The destruction of houses and other buildings found in urban areas also produces large amounts of debris, and regions of extensive debris coverage are easily detectable using multispectral methods [19]. Further analysis of the relationship between landcover and high -resolution damage detection (e.g., [18]) is an opportunity for further study.
Image selection is another notable factor for consideration of the results presented here. Molthan et al. (2014) found that NDVI decreases in the weeks following a tornado became marginally more pronounced with increases in EF-scale damage, but variance was relatively high within each damage level [14]. The authors suggested that non-meteorological processes such as cleanup efforts are likely the cause of these NDVI decreases [14]. Thus, detailed post-event damage analysis should ideally be conducted with imagery obtained as soon as possible after the event, since non-meteorological processes that inevitably occur in the days after the event affect the landscape in complex ways that might mask fine-scale variations in damage. This may especially be true in urban areas, where anthropogenic activities are most concentrated. While some cleanup undoubtedly occurred in the days immediately following the tornadoes examined in this study, the three-day lag between the tornado event  [14,19,20]. Kingfield and de Beurs (2017) provided a detailed exploration of vegetation recovery in the five years following the 27 April 2011 tornado outbreak using Landsat imagery, finding that substantial recovery, relative to nearby undamaged regions, did not begin until six months after the event [18]. Future studies could conduct a similar examination of short-and medium-term reflectance changes in damaged areas over small spatial scales, as this could shed light on the importance of data latency in post-event image selection.
An advantage of this study's methods is in their ease of being replicable. Detailed studies of tornado damage using aerial imagery often involve manual inspection and tree fall digitization [11,27,28], which produce very detailed datasets, but are also very time intensive and require aerial imagery of the affected area, which can be expensive to obtain. NDVI, on the other hand, is a very common metric in the Earth observation literature and is easy to calculate at the expense of some level of detail. Imagery obtained from UAVs and satellites with a very fine (< 10-m) spatial resolution will increasingly be used for future damage assessments of this nature [12,13,43,44].
This approach can also be used to refine ground-based damage surveys. One example of this is the region of displacement between NDVI decreases and the survey-defined path polygon shown in Figure 7. Such displacement between remotely sensed damage and ground-level surveys, especially perpendicular to storm motion, is not uncommon in rural regions [14,16,45]. While Karstens et al. (2013) [28] noted tree falls associated with rear flank downdraft (RFD) surges to the right of the 22 May 2011 Joplin, Missouri, tornado track, these surges caused substantially less damage than was found within the path of the Joplin tornado. Given the severity of the damage shown in Figure 7 outside of the path relative to the damage inside it, Figure 7 likely shows damage from the tornado itself that was not included in the survey-defined path because it would have been nearly impossible for surveyors to notice or access from ground level.

Conclusions
This study examines damage from three tornadoes using imagery from Sentinel-2 MSI with a 10-m spatial resolution. Pixel values from post-event NDVI, NDVI differences, and eight principal component images are correlated with estimated wind speeds at damage indicators noted in post-event surveys. The strongest correlations with these wind speeds are in the NDVI-difference image and one of the principal component images (PC7), after applying a five-pixel-by-five-pixel filter window to each image to diminish statistical noise. Random points are then placed within each of the three tornado paths, and NDVI-difference and PC7 values at the location of these points are plotted against the distance from each respective point to the center of the tornado path. The shape of a smoothing spline fitted to the NDVI-difference values in the Beauregard-Smiths Station tornado path resembles the pattern of ground-relative maximum wind speeds across the width of a tornado in the Rankine vortex model. The most intense damage detected with these methods occurred within 130 m of the tornado path center, which may serve as a rough estimate for R max for this tornado. However, the shapes of the smoothing splines fitted to PC7 values in the Beauregard-Smiths Station tornado, as well as both NDVI-difference and PC7 values in the weaker Davisville and Weedon Field tornadoes, are relatively flat across the width of the tornado paths. It is concluded that, with little cloud contamination, imagery with ≤ 10-m spatial resolution such as that from Sentinel-2 MSI can be used to approximate damage intensity of vegetated areas and fine-tune the extent of tornado paths delineated in on-the-ground surveys. Pixel-level atmospheric corrections may be able to provide even more detail on ground-level damage patterns, as may UAV imagery, which is obtained from a low altitude [13]. The techniques and findings presented in this study may be of value to future research on vegetative response to tornado damage (e.g., [1,11,18]), as well as comparisons between ground surveys, radar observations, and passively sensed imagery (e.g., [12,16,35,39,42,46]). As high-resolution UAV and satellite imagery becomes increasingly available, these multi-platform studies will be able to make detailed observations of tornado dynamics to aid advances in modeling, building safety, and similar efforts.