North Alabama is among the most tornado-prone regions in the United States and is composed of more spatially variable terrain and land cover than the frequently studied North American Great Plains region. Because of the high tornado frequency observed across north Alabama, there is a need to understand how land surface roughness heterogeneity influences tornadogenesis, particularly for weak-intensity tornadoes. This study investigates whether horizontal gradients in land surface roughness exist surrounding locations of tornadogenesis for weak (EF0–EF1) tornadoes. The existence of the horizontal gradients could lead to the generation of positive values of the vertical components of the 3D vorticity vector near the surface that may aid in the tornadogenesis process. In this study, surface roughness was estimated using parameterizations from the Noah land surface model with inputs from MODIS 500-m and Landsat 30-m data. Spatial variations in the parameterized roughness lengths were assessed using GIS-based grid and quadrant pattern analyses to quantify observed variation of land surface features surrounding tornadogenesis locations across spatial scales. This analysis determined that statistically significant horizontal gradients in surface roughness exist surrounding tornadogenesis locations.
The United States leads the world in tornado occurrence, averaging around 1000 events recorded annually (NCEI 2016). Within the United States, some locations experience a disproportionately higher frequency of tornadoes, which poses a higher risk to inhabiting citizens who suffer from the economic and societal impacts. North Alabama is one of these locations found by Coleman and Dixon (2014) to be among the most at risk areas for significant tornadoes [enhanced Fujita scale 2 (EF2) and greater] in the country. In addition, north Alabama has suffered from the most prolific tornado outbreaks in recorded history including the 3 April 1974 and 25–27 April 2011 outbreaks (Knupp 2014).
Within north Alabama, the NOAA National Weather Service Weather Forecast Office (WFO) operational forecasters, broadcast meteorologists, private industry weather companies, and the public have observed spatial patterns in tornado distribution and frequency. A spatial density map generated using the NOAA Severe Weather GIS (SVRGIS) dataset EF0–EF5 tornado tracks observed within the University of Alabama in Huntsville (UAH) operated Advanced Radar for Meteorological and Operational Research (ARMOR) coverage area is shown in Figure 1 (NOAA 2015). For tornado events captured prior to the 1973 enactment of the Fujita scale, tornado intensity was determined according to a property loss threshold if monetary damage information was available (NOAA 2015). In 2007 when the EF scale began being used, the Fujita scale rankings were converted to the EF scale and catalogued within the dataset. The generated map indicates two areas with high densities, or high frequencies of tornadoes for a given area. One area is located over north-central Alabama into southern Tennessee, while a second is located over the Cumberland Plateau extending into Marshall, Jackson, and DeKalb Counties. Compared to the frequently studied Great Plains region, north Alabama has greater variability in topography and vegetation, posing questions on how the lowest several hundred meters of the atmosphere, that is, the surface layer to lower boundary layer are influenced by horizontal gradients in the land surface roughness and how the tornadogenesis process may be affected.
Tornadogenesis occurs within both supercell storms and quasi-linear convective systems (QLCS). Since a significant majority of high-intensity (EF2–EF5) tornadoes occurs within supercell storms, previous research has focused on the development of low-level vertical vorticity within supercell storms via horizontal gradients in buoyancy at low levels (Davies-Jones 2015; Markowski et al. 2018, and references cited therein). The development of mesovortices within QLCSs has also been related to the development of horizontal gradients in buoyancy (Atkins and Laurent 2009), development of rotors by surface friction (Schenkman et al. 2012), and horizontal shearing instability (Conrad and Knupp 2019). This study focuses on weak-intensity tornadoes (EF0–EF1) since they are more numerous and are potentially influenced by horizontal variations in land surface roughness, according to the primary hypothesis defined in the following. For both storm types, Markowski et al. (2018, p. 3623) has stated, “A better understanding is needed of what can trigger a storm in a favorable environment to suddenly make a tornado at a particular stage in its evolution.”
Decades of land surface remote sensing, tornado event observations, and geospatial tools enable variations in land surface features and tornado events to be assessed geostatistically as a first step toward characterizing how they affect tornadogenesis. This study integrates land surface remote sensing parameters, tornadogenesis locations, geospatial analysis, and the vertical component of the 3D vorticity equation to propose a novel, interdisciplinary approach to studying spatial relationships between land surface roughness variation and tornadogenesis. The goal of this study was to determine whether horizontal gradients in land surface roughness exist in the surrounding locations of tornadogenesis for weak (EF0–EF1) tornadoes. The existence of the horizontal gradients could lead to the generation of positive values of the vertical components of the 3D vorticity vector near the surface in the lowest several hundred meters of the atmosphere (surface layer to lower boundary layer) that may aid in the tornadogenesis process. We accomplish this goal by 1) parameterizing surface roughness using satellite remote sensing and land surface model schemes, 2) running a spatial pattern analysis to assess spatial roughness variation surrounding locations of tornadogenesis, 3) conducting a quadrant pattern analysis to statistically identify significant relationships between horizontal roughness gradients and tornadogenesis locations, and 4) conduct this study using 30- and 500-m land surface data to assess if the same patterns are observed at varying spatial scales. The developed methods present a comprehensive and novel GIS-based approach for quantifying variations in surface friction using land surface roughness parameterization schemes and spatially assessing its variation in north Alabama surrounding EF0–EF1 tornadogenesis locations.
2.1. Defining roughness length
Land surface roughness is a measure of the dynamic interaction between the land surface and the wind profile near the surface used to quantify the reduction in wind speed by drag force and friction (Raupach 1994; Holton 2004). The roughness of the land surface depends on the vertical deviation from an ideal, homogeneous surface at neutral stability created by variable land cover and topography. Within the logarithmic wind profile equation, the roughness length term z0 is quantified in units of length. Equation (1) shows the terms of the logarithmic wind profile equation rearranged to equate to z0 (m):
where u(z) is the mean flow velocity profile (m s−1) over a surface, u* is the friction velocity (m s−1), k is the von Kármán constant (0.4), and z is the height above ground level (m). The friction velocity u* shown in Equation (2) accounts for shearing stress near the ground:
where τ is the shear stress (kg m−1 s−2), and ρ is the density of a fluid (kg m−3). The larger the value of z0, the greater amount of frictional drag exerted on air moving over the land surface, thus corresponding to reduced wind shear.
2.2. Vorticity and tornadogenesis
Between areas with different magnitudes of z0 there exist horizontal gradients in low-level flow moving over these surfaces that are most notable between areas of high (forest) and low (water, smooth grass surfaces) z0. Coleman and Knupp (2009) proposed that when surface winds move over heterogeneous land surfaces where the gradient in z0 is located normal to the direction of low-level flow (lowest several hundred meters of the atmosphere), then positive vertical vorticity could be produced above the surface. This process produces nonzero values of the vertical component of the 3D vorticity vector at the surface that may enhance tornadogenesis. Similar interactions create turbulent eddies and mixing when water or air moves over obstacles such as trees, buildings, and rocks that induce frictional drag on the fluid (air, water) creating small vortices (Stull 1988; Bou-Zeid et al. 2004). This can be explained mathematically by the friction term in the simplified vertical component of the 3D vorticity equation ζ (Bluestein 1992):
where all other terms are assumed to be zero, is the vertical unit vector, and is the curl of the friction force (Coleman and Knupp 2009). Figure 2 illustrates this process, where represents the horizontal gradient in z0 along a rough forested area (z0 = 0.5 m) adjacent to a smoother grassland (z0 = 0.1 m) area. The arrows represent differing velocity vectors and the resultant curl of the friction force as low-level flow moves over the area with a horizontal gradient in z0. The resultant curl, shown by the curved arrow, indicates a location where tornadogenesis could be aided.
Although both the horizontal and vertical components of the vorticity vector require consideration, we focus herein on the generation of the vertical component of vorticity ζ within the lowest several hundred meters of the atmosphere, that is, the surface layer to lower boundary layer. It is recognized that the magnitude of horizontal vorticity associated with vertical wind shear, and ζ within supercell mesocyclones or QLCS mesovortices, is 10−2 s−1. The magnitude of ζ produced by horizontal gradients in z0 associated with a water–land interface, documented by Asefi-Najafabady et al. (2010) for a moderate wind case (e.g., 5 m s−1 wind speed), can exceed 10−3 s−1 over a significant depth of the ABL (Figure 3). This value could conceivably approach 5 × 10−3 s−1 for stronger surface-layer flow on much smaller scales on the order of 100 m. For example, a surface-layer wind differential of 5 m s−1 over a 500-m horizontal scale yields a ζ magnitude of 10−2 s−1. Values of around this magnitude are common along the leading edge of QLCSs, and within some well-defined cold fronts and outflow boundaries associated with mesocyclones (Marquis et al. 2007).
Laboratory simulations and numerical modeling studies have also investigated the influence that land surface roughness heterogeneities have exerted on surface flow, and have demonstrated that changes in the land surface contribute to production of vorticity and generation of tornado-like flows, illustrating a need for further investigation (Dessens 1972; Avissar and Pielke 1989; Natarajan and Hangan 2009; Kellner and Niyogi 2014; Liu and Ishihara 2016; Wang et al. 2017).
2.3. Methods to derive z0
Approaches to derive and represent z0 for research and modeling purposes include using lookup tables, satellite-derived land surface and vegetation parameters, and hybrid approaches. Traditional atmospheric modeling techniques implement lookup tables that equate a land-cover class to a representative z value derived from field observations (Borak et al. 2005; Jasinski et al. 2006; Zheng et al. 2014). These approaches fail to capture variation within land-cover classes and account for seasonal vegetation changes, resulting in errors in heat fluxes and skin temperature values used in weather forecasting models (Zheng et al. 2014). Alternatively, remote sensing methods have been used to relate satellite observations and ancillary data to better capture the spatial and temporal variability of z0. Hybrid approaches using the Moderate Resolution Imaging Spectroradiometer (MODIS) vegetation indices and lookup tables have been proposed by Wu et al. (2015), Borak et al. (2005), Jasinski et al. (2006), and Jasinski et al. (2005). Backscatter caused by variation in the land surface measured by spaceborne radiometers, such as QuikSCAT, have been used to quantify roughness (Goodberlet and Mead 2014; Bergen et al. 2009). Zhang et al. (2004) implemented a multisensor approach using RADARSAT synthetic aperture radar (SAR), Landsat 5 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper Plus (ETM+) data. The SAR and lidar data approaches have been developed to capture the roughness of land-cover types such as bare soil, ice, and urban areas, which are difficult to model (Chen et al. 2009; Rees and Arnold 2006; Rivas et al. 2006).
Atmospheric models implement z0 parameterization schemes to account for energy transfer and interaction between the land surface and atmosphere using lookup table values and satellite retrievals. The Noah land surface model (LSM) (Chen et al. 1996; Ek et al. 2003) is widely used as the land surface component for regional and global weather forecasting models at the National Centers for Environmental Prediction and by the Weather Research and Forecasting (WRF) Model at the National Center for Atmospheric Research (Zheng et al. 2014). The Noah LSM parameterization scheme may fail to accurately capture z0 ranges given its reliance on land-class-based lookup table values and time-series-derived z0 ranges, which in some cases offer constant values of land-cover classes. Nevertheless, the Noah LSM scheme was used to quantify z0 across the study area because it offers a hybrid approach using land surface remote sensing data and lookup table values that are implemented widely across modeling applications, and accounts for seasonal vegetation variation (Zheng et al. 2014; Chen and Zhang 2009; Ek et al. 2003; Niu et al. 2011).
2.4. Study area and period
The ARMOR is a scanning dual-polarimetric Doppler radar operated by the UAH Department of Atmospheric Science (University of Alabama in Huntsville 2004). ARMOR is located at the Huntsville International Airport and typically scans a 120-km range across north Alabama and portions of Tennessee and Georgia. The study area (Figure 4) is located west of the southern Appalachian Mountains covering portions of the Cumberland Plateau and Tennessee River valley with topography ranging approximately 80–700 m in elevation.
The study area is classified as humid subtropical under the Köppen climate classification scheme (Kottek et al. 2006). During the March–May severe weather season, the climate is warm and humid with average 3-month temperatures near 60°F and a relative humidity of 70% (NCDC 2018). In the fall and, primarily, the spring, severe weather and tornadoes are most prevalent; however, thunderstorms occur throughout the summer with 50–60 thunderstorm days observed each year (Burt and Stroud 2007). During the spring tornado months, conditions are very dynamic with respect to the emergence of vegetation following the winter months. The tornado events used in this study (shown in Figure 4 as inverted triangles) included all tornadoes that occurred within 120 km of the ARMOR radar over the 2005–14 period during March–May, the climatological peak in tornado occurrence within the state of Alabama (NWS 2018).
3. Data and methods
This study combines geospatial data, satellite remote sensing retrievals, and the Noah LSM parameterization scheme to capture z0 variation within the ARMOR coverage area and assesses its relationship to tornadogenesis locations of weak-intensity tornadoes (EF0–EF1). The selected datasets represent locations of tornadogenesis, land cover, and vegetation cover. The Noah LSM was used to parameterize z0, followed by a geostatistical analysis to assess the spatial relationships between tornadogenesis locations and z0 variability. A z0 layer was created at two spatial resolutions to assess relationships between larger-scale (500 m) and smaller-scale features (30 m) surrounding tornadogenesis locations. An overview of the data and methods used in this study is shown in Figure 5 and outlined in the proceeding section.
3.1. Tornado database
The NOAA Storm Prediction Center (SPC) SVRGIS tornado-track dataset (NOAA 2015) and the unpublished UAH Department of Atmospheric Science tornado event catalog (University of Alabama in Huntsville 2014, unpublished material, accessed 1 June 2014) were used to identify tornadogenesis locations and information on maximum tornado intensity, wind speed, fatalities, and injuries. Each dataset is compiled from National Weather Service storm reports; however, the UAH catalog is a spatial and temporal subset of the SVRGIS dataset containing only tornado events that initiated within the ARMOR radar coverage area from 2005 to 2014. These two datasets were quality controlled by cross comparing event information to identify errors before plotting the tornadogenesis point locations. The points were then filtered to identify EF0–EF1 tornadogenesis events during the study period that initiated within the study area, leaving 131 tornadogenesis points.
Several tornadoes produced during the 25–27 April 2011 tornado outbreak and other events left behind tornado scars (long-lasting damage to vegetation, especially trees) that greatly altered the land cover and were captured by the MODIS and Landsat retrievals. The z0 parameterization scheme is dependent on representative land cover and vegetation for the spring tornado seasons of 2006 and 2011, thus tornadogenesis locations occurring between 2005 and 2008 correspond to 2006 land cover and those occurring between 2009 and 2014 correspond to 2011 land cover. Changes to the underlying land features affect the accuracy of the z0 representation particularly for tornadogenesis events that occurred before the tornado scars were made and captured within the 2011 land cover (from 2009 to 27 April 2011), thus, it was necessary to correct for these changes in z0. The NOAA geospatial tornado-track survey data for 25–27 April 2011 produced by NWS Birmingham (2011) offers detailed tornado swath information. These data were used to identify significant (EF4–EF5) tornado events that occurred within 0–4 km of a tornadogenesis location. Tornadogenesis events within the 2011 sample group occurring from 2009 to 27 April 2011 within a 4-km area of a significant tornado event were omitted from this study. A total of 6 tornadoes from the 2011 sample group were removed reducing the final sample size to 125 tornadogenesis points.
3.2. Land cover and reclassification
Data acquired from the Multi-Resolution Land Characteristics Consortium (MRLC) National Land Cover Database (NLCD) served as the high-resolution 30-m land-cover layer (Homer et al. 2011; Fry et al. 2006). NLCD data are derived from Landsat imagery and use a 16-class land-cover classification scheme based on the Anderson et al. (1976) scheme containing an accuracy of 84% (Wickham et al. 2013). Spatial variations in land-cover type are minimal seasonally and annually, thus NLCD land cover is produced in 5-yr intervals. The 2006 and 2011 datasets provided representative land cover for two time spans of tornadogenesis points: 2005–08 for the 2006 dataset and 2009–14 for the 2011 NLCD dataset.
The MODIS collection 5 level 3 annual 500-m land-cover composite product, MCD12Q1, is derived from the MODIS sensor aboard the Aqua and Terra satellites, with class accuracy between 72% and 75%, provided a medium-resolution land-cover representation (Friedl et al. 2010; Strahler et al. 1999; Friedl and Sulla-Menashe 2015). MODIS land-cover data for 2006 and 2011 were acquired from the Land Processes Distributed Active Archive Center (LP DAAC) to align with the NLCD land-cover data. The IGBP classification scheme was selected because it was developed to represent land classes most useful across modeling disciplines Strahler et al. (1999).
The Noah LSM combines remote sensing measurements and z0 values associated with the 17-class MODIS IGBP or 24-class USGS land-cover classification schemes. The NLCD land-cover data were reclassified to match the 24 USGS land-cover classes, with the exception of the urban areas. The NLCD urban land-cover classes capture four development classes that vary widely in z0: open space, low intensity, medium intensity, and high intensity. Because z0 values for these classes are not part of the Noah LSM USGS classification scheme, associated z0 value ranges were determined from the Vickery et al. (2006) and Simpson et al. (2012) predefined parameters.
3.3. NDVI and GVF data, analysis, and corrections
The vegetation of the study area was represented using normalized difference vegetation index (NDVI) datasets derived from two satellite instruments. The Terra MODIS, version 3, 500-m MOD13A1 product (Didan 2015) was acquired from the LP DAAC for the full years of 2006 and 2011. This NDVI data layer contains an accuracy within ±0.025 and was quality controlled to remove pixels under clouds and cloud shadow that cause brightening and darken effects, respectively, reducing the accuracy of atmospheric correction and computed NDVI values (Huete et al. 1999; MODIS Land Team 2016; Zhu and Woodcock 2012).
The U.S. Geological Survey (USGS) Landsat Surface Reflectance Level-2 Science Product data (USGS 2011, 2006) were acquired from the Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) interface (USGS 2019a,b). For this data product, the 6S radiative transfer model (Vermote 1997) was applied to Landsat 4–7 retrievals, from which NDVI and quality assessment (QA) bands are derived. For this study, 286 Landsat 7 ETM+ and Landsat 5 TM images with 60% cloud clover or less were acquired for 2006 and 2011. The NDVI data were quality controlled using the QA bands to remove contaminated pixels caused by oversaturation, cloud cover, and cloud shadow established by Zhu and Woodcock (2012), Arvidson et al. (2001), Arvidson et al. (2000), and Simpson and Stitt (1998).
The Noah LSM uses green vegetation fraction (GVF) to quantify and control water and energy transfer between the land surface and the atmosphere and capture seasonal variations in z0 over vegetated areas (Zeng et al. 2000). To determine z0 of a particular area, the GVF is determined using
where NDVI is the NDVI value of a specific pixel at a point in time, NDVImin is the constant minimum NDVI value for bare soil (0.01), and NDVImax is the maximum NDVI for a land-cover class over a period of time.
A temporal analysis of the NDVI values was performed to determine the NDVImax for USGS and IGBP land-cover types in the study area. Two full years of MODIS and Landsat NDVI data were analyzed to ensure values captured vegetative ranges across multiple years. NDVI pixel values for each land-cover class were extracted for the study area to which a histogram thresholding method by Zeng et al. (2000) and Gutman and Ignatov (1998) was applied to account for pixel contamination near land-cover -class boundaries. The NDVImax values from the sample were selected as follows: 90th-percentile value for sparsely vegetated areas (e.g., urban and bare soil) and the 75th-percentile value for vegetated areas. The highest NDVImax values across 2006 and 2011 are shown in Table 1.
The NDVI variable in Equation (4) was sampled from NDVI values for the climatological peak spring tornado season for north Alabama (March–May). Representative spring layers for the Landsat and MODIS NDVI were generated for 2006 and 2011 by determining the maximum NDVI for each pixel across the spring months. The maximum spring NDVI provided the best NDVI representation for the study area that corrected striping from the 2003 Landsat 7 ETM+ sensor scan line corrector failure without interpolating Landsat retrievals. Once the NDVI was determined, GVF was calculated for 2006 and 2011 shown in Figure 6.
3.4. Noah LSM parameterization and z0 calculations
The Noah LSM uses two types of roughness values to account for the heat transfer Z0h between the land surface and low-tropospheric levels, and momentum transfer using a momentum roughness length Z0m. This study uses the generalized version of the Monin–Obukhov similarity theory that assumes Z0h = Z0m and implements the Z0m parameterization scheme, from here on out we refer to roughness as Z0m (Brutsaert 1982; Niu et al. 2011). The Chen and Zhang (2009) implementation of the Noah LSM Z0m parameterizations scheme was selected since it was found to consistently perform better over changing surface conditions and seasons (Zheng et al. 2014). This approach uses the aforementioned Z0m lookup table values (Table 1) and the GVF to solve for Z0m:
where GVF is the spring layers calculated in the previous section, Z0m,min is the minimum seasonal roughness length value in meters, and Z0m,max is the maximum seasonal roughness length value in meters. This equation was applied to generate four corresponding Z0m layers based on each pixel’s land-cover class using the GVF and the variables outlined in Table 1. The final Z0m layers for the 2006 and 2011 Landsat and MODIS data are shown in Figure 7.
3.5. Spatial analysis architecture
A two-part methodology was developed to assess spatial relationships between locations of tornadogenesis and surrounding Z0m variation and patterns. A gridcell analysis was conducted to acquire a qualitative overview of the variation in Z0m observed by the MODIS and Landsat instruments to evaluate patterns captured across different satellite instruments on equal-area grids. To quantify the observed patterns, a quadrant pattern analysis was used to determine the significance of the observed patterns. ArcGIS was used in both these analyses to construct the grids and quadrants and perform the geospatial statistics.
3.5.1. Gridcell analysis
An analytical constraint was set to focus this investigation on local mesoscale features surrounding each tornadogenesis, or tornado-track initiation point location. Mesoscale interactions occur within a 2–2000-km horizontal scale; therefore, anything below this threshold is considered microscale (Orlanski 1975). A 4-km threshold surrounding each tornadogenesis point was selected based on 4 km being the coarsest spatial resolution used in convection allowing mesoscale models such as the WRF Nonhydrostatic Mesoscale Model and Advanced Research WRF, and this threshold being close to the approximate microscale where the proposed tornadogenesis processes occur (Markowski and Richardson 2009). The minimum gridcell dimension was set to 500 m, the spatial resolution of the MODIS data, enabling the comparison of the different resolution data across equal spatial scales. Grids with 500 m × 500 m cells within an 8 km × 8 km area surrounding each tornadogenesis point were created and oriented along each tornadogenesis point associated tornado track bearing direction θ, which is used as a proxy for surface-level flow recognizing there is a 45° offset between surface flow and tornado-track orientation. Figure 8 provides a schematic of one of the 125 resultant grids, where the tornadogenesis location is shown as the black dot, surrounded by 256 grid cells within a 64-km2 area. For each tornadogenesis point grid, the Z0m pixel values falling within each grid cell were averaged for the Landsat and MODIS data. These 125 MODIS and 125 Landsat grids were then averaged across the sample for all corresponding gridcell locations to generate comparable results for the two Z0m datasets outlined in the proceeding section.
3.5.2. Quadrant pattern analysis
Square quadrants surrounding each tornadogenesis point were separately created (in addition to the grids) using the same 500-m intervals (Δr and Δd shown in Figure 9) oriented along θ from 500 to 4000 m in size. To generate the quadrants, the coordinates of two points along the θ axis were determined at varying along-track distances from the tornadogenesis location at Δd ahead of and before the location. A schematic of this process is shown in Figure 9: the point of tornadogenesis is centered, surrounded by a point forward (after) the tornadogenesis location and a point backward (before) tornadogenesis oriented along θ. From these lines and using the same 500-m intervals, square quadrants were created along each axis to a cross-track distance r, two before the tornadogenesis point (backward-left and backward-right quadrants), and two after tornadogenesis (forward-left and forward-right quadrant). A UTM projection was critical for determining the coordinate locations used to construct the quadrants because it offers a higher precision than degree-based projections, which reduced the resultant offset between the orientation of the constructed quadrant and the original tornado track. Similar to the grid analysis, for each quadrant the mean Z0m values of pixels falling within each of the areas were sampled across all tornadogenesis locations and repeated across the four Z0m layers for each variable quadrant size.
3.5.3. Quadrant pattern groups
Separately, a pattern analysis was used to quantify spatial patterns and variation observed in the Z0m grid analysis to determine relative differences in Z0m between six predefined pattern groups, whether these differences were significant, and if significant, the relevant spatial scales they differ. The mean Z0m values collected for each quadrant were combined to test six possible quadrant patterns combinations surrounding the tornadogenesis locations shown in Figure 10. For each pattern combination, the average quadrant values of Z0m were separated into two groups to compare the relative difference in the mean quadrant Z0m for each pattern combination where gray represents group A and black represents group B. The quadrants were combined into two pattern groups by averaging the per-quadrant mean Z0m within each pattern group. The resultant mean values of Z0m for each pattern group were differenced across the 125 tornadogenesis points and spatial scales. Eight iterations of pattern differences were run for each of the six patterns for each tornadogenesis point, generating 1000 data values for each pattern group.
where K represents the kernel, xi is the univariate independent variable, h is the smoothing parameter called the bandwidth, n represents the sample size, and Kh is the scaled kernel. Section 4 will discuss the results of both the grid and the quadrant pattern analyses.
4. Results and discussion
4.1. Grid analysis results
The grid analysis enabled a qualitative assessment of spatial variations in Z0m across two spatial scales. For each grid cell, the average Z0m pixel value was extracted and averaged across the sample for MODIS and Landsat independently, shown in Figure 11, where the location of tornadogenesis (black inverted triangle), θ (solid arrow), and surface flow (dashed arrow) are depicted. The grid analysis revealed the mean Z0m for Landsat to be higher than MODIS with value ranges from 0.180 to 0.228 m and 0.128 to 0.176 m, respectively; however, the range for each dataset is similar, falling near 0.05 m. The spatial variation in Z0m appears similar for each sample with lower Z0m located in the forward-right quadrants and higher Z0m located in the other three quadrants. Most notable are the highest Z0m located within the backward quadrants, primarily the backward right and along the bearing axis. The differences between the MODIS and Landsat representations of Z0m are resultant from the spatial resolutions, where the MODIS representation shows larger-scale features, and more variation in Z0m is observed for the Landsat data across the same general pattern. A gradient in Z0m is shown to exist in the Landsat and MODIS layers perpendicular to the surface-layer wind. The tornadogenesis points are also very close to the minimum in Z0m supporting the hypothesis in that a horizontal gradient exists with the potential of generating positive values of the vertical components of the 3D vorticity vector that may aid in tornadogenesis. This is consistent with the conceptual model shown in Figure 2.
4.2. Quadrant pattern analysis results
The KDF distributions of each pattern group difference for each tornadogenesis location were computed independently for the Landsat and MODIS data to quantify observations made from the grid analysis. Plots were generated that depict a series of distribution lines, each representative of the pattern group difference values calculated across the 125 tornadogenesis locations, eight spatial scales, and the six pattern groups.
4.2.1. Landsat results
The Landsat KDF plots (Figure 12) reveal the distribution of pattern group difference values to become more widely distributed as the quadrant size increases. This variation likely results from the larger area sampled from which more Z0m pixel values are averaged yielding greater variability and wider distributions. In addition, noticeable skewness in both larger and smaller quadrant scales are observed across several of the pattern groups, most notable in 1, 2, 3, and 6. This relative skewness implies larger differences in Z0m between pattern groups across certain spatial scales further indicating horizontal gradients in Z0m.
The mean, standard deviation, skewness, and kurtosis were computed to further investigate the statistical significance of the KDF pattern group difference observations. A one-sample two-tailed t test was performed to identify pattern groups meeting a 95% significance level with a null hypothesis that there is no difference between the two pattern groups, or that the resulting Z0m difference is zero (implying a homogeneous flat surface). Otherwise, those that test as significant imply that a notable horizontal gradient exists between the two pattern groups. The statistical results for Landsat are provided in Table 2 where boldface values indicate statistical significance.
The two-tailed t test revealed two pattern groups for the Landsat data with p values meeting the 95% significance level; pattern group 1 at 500 m had a p value of 0.047 and Pattern group 6 at 4000 m the p value was 0.045. Further examination reveals that the distribution of values is skewed toward the positive end, meaning quadrant group A has a higher mean Z0m value than quadrant group B. In each of these pattern groups, group A is located to the left of the tornado track bearing direction θ, which was observed to be significantly higher for these specific spatial scales and pattern groups, hence rejecting that null hypothesis that a homogeneous surface surrounds the tornadogenesis locations. This result supports the assessed hypothesis revealing Z0m gradients conducive for the formation of positive vertical vorticity to form near the surface, and consistent with the findings from the grid analysis (Figure 11).
4.2.2. MODIS results
The MODIS KDF plots in Figure 13 reveal wider distributions than the Landsat plots and less variation across spatial scales. The tail end of most distributions reveal a smaller secondary peak for all patterns, however, the KDF distributions maximums remain centered around zero. Less variability between spatial scale and pattern groups can be attributed to the lower resolution of the MODIS data, which captured larger-scale features in the grid analysis compared to Landsat.
The two-tailed t test performed on the MODIS data did not reveal any pattern combinations that met the 95% significance level in Table 3. These results were expected given the lack of qualitative skewness in the KDF plots and reveal the 500-m MODIS data likely have a spatial resolution too coarse to study Z0m variation within an 8 km × 8 km area, failing to reject the null hypothesis.
4.3. Discussion of results
The grid and quadrant pattern analyses were performed to qualitatively and quantitatively evaluate horizontal variation in Z0m, and statistically identify the spatial patterns and scales with significant differences between pattern groups. The grid analysis revealed a decrease in Z0m within the forward-right quadrant in the MODIS and Landsat data, with higher Z0m in the backward quadrants. These results suggest that across the sample of EF0–EF1 tornadogenesis locations, there may be gradients in Z0m conducive for the formation of positive vertical vorticity near the surface. Building on the grid analysis, the quadrant pattern analysis was developed to evaluate six quadrant pattern group differences across the MODIS and Landsat datasets and within varying areas surrounding each tornadogenesis location. To accomplish this, KDF plots for each pattern group and spatial scale were generated, followed by a series of statistics to test for significance. The resultant Landsat KDF plots revealed relative skewness between pattern difference groups; however, as the spatial scale increased, less skewness and wider distributions were observed. In comparison, the MODIS KDF plots revealed wider distributions than Landsat with little change in the distribution width as the quadrant size increased. In addition, secondary peaks in MODIS KDF values were observed; however, the maximum of each curve was centered on zero with little relative skewness observed. The two-tailed t test was performed using a 95% significance level to identify significant differences between the pattern groups tested. The test revealed two pattern groups and spatial scales with significant differences in Z0m for the Landsat data; however, the MODIS statistics were unable to replicate these findings. The Landsat test identified pattern group 1 at 500 m and pattern group 6 at 4000 m to have value distributions skewed toward the positive end, meaning quadrant group A, on the left side of the flow, has a higher mean Z0m value than quadrant group B. These Landsat results found statistically significant differences in Z0m, supporting the grid analysis, and finding a horizontal gradient in surface roughness. This finding supports the theory that horizontal surface roughness gradients coinciding with weak-intensity tornadoes exist that may influence the vertical component of the 3D vorticity vector near the surface. More work is needed to investigate at what magnitude these horizontal gradients influence the 3D vorticity vector and tornadogenesis. Results indicate that within an 8 km × 8 km area surrounding each tornadogenesis location, Landsat was able to detect more variability and smaller-scale detail because of its higher spatial resolution. The MODIS data proved useful for identifying larger-scale variation in Z0m; however, the variation was not statistically significant at smaller scales where resultant KDF distributions were found to be similar across the spatial scales and patterns tested. For future work, the area surrounding each tornadogenesis location should be expanded to assess results across larger grid and quadrant areas that fall within the mesoscale.
Errors in the outlined findings may result from inaccuracies in the data and methodologies implemented. First, the Noah LSM parameterization scheme may result in errors in the resultant Z0m layers given its dependence on a land-class-based lookup table. In addition, the accuracy of the NDVI and land-cover classes could affect the Z0m layer derived because of errors in the land-cover classification schemes, and oversaturation of NDVI values in highly vegetated areas. In final, it must be noted that the tornado dataset used does not represent all events that occurred within the study area given is it reliance on human reports, thus these data are skewed toward more populated areas with human observations and access, keeping in mind that errors in the location of the actual tornadogenesis point exist, but are unknown. Future efforts should test and compare results from alternative Z0m parameterization and remote sensing approaches. In addition, it must be noted that tornadoes are highly dynamic events that rely on atmospheric as well as other environmental and land surface parameters. This first-order analysis was conducted to identify horizontal gradients on the land surface without considering the other atmospheric dynamics at hand. Further research is needed to incorporate more detailed modeling efforts and case study analysis to expand on this work.
The goal of this study was to determine whether horizontal gradients in land surface roughness exist surrounding locations of tornadogenesis for weak (EF0–EF1) tornadoes. This study found evidence of statistically significant horizontal gradients in Z0m surrounding EF0–EF1 tornadogenesis locations. Observations acquired from Landsat at 30-m spatial resolution revealed statistically significant horizontal gradients Z0m that may be favorable for supporting tornadogenesis. It is possible that the weaker tornadoes occur in environments that are marginally supportive of tornadogenesis and that the roughness gradient provides the forcing required for tornadogenesis in such cases. Future studies should address the physical processes using idealized numerical simulations. Additional research is needed to test the developed methodology across other geographic regions to validate this study’s findings in modeling efforts. Other remote sensing and modeling schemes should be incorporated into further analysis to ensure land surface variation is captured effectively. Finally, a primarily outcome of this work was the need to expand the grid and quadrant area sampled around each tornadogenesis location in order to investigate larger-scale features captured by coarser-resolution data. These results present a case to further investigate the effect variation in land surface features has on tornado dynamics. Increased understanding of the interactions between the atmosphere and land surface will improve weather forecasting and modeling and may, in the future, lead to tornado hazard mapping techniques to aid in disaster management efforts.
We thank the University of Alabama in Huntsville Department of Atmospheric Science and Earth System Science Center for supporting this research. Additionally, thanks to Kel Markert for his technical programming support in executing this study.