Predicting Risks of Severe Windstorm Damage to Southeastern U.S. forests


 ContextThe southeastern U.S. experiences tornadoes and severe thunderstorms that can cause economic and ecological damage to forest stands resulting in loss of timber, reduction in short-term carbon sequestration, and increasing forest pests and pathogens. ObjectivesThis project sought to determine landscape-scale patterns of recurring wind damages and their relationships to topographic attributes, overall climatic patterns and soil characteristics in southeastern forests. MethodsWe assembled post-damage assessment data collected since 2012 by the National Oceanic and Atmospheric Administration (NOAA). We utilized a regularized Generalized Additive Model (GAM) framework to identify and select influencing topographic, soil and climate variables and to discriminate between damage levels (broken branches, uprooting, or trunk breakage). Further, we applied a multinomial GAM utilizing the identified variables to generate predictions and interpolated the results to create predictive maps for tree damage. ResultsTerrain characteristics of slope and valley depth, soil characteristics including erodibility factor and bedrock depth, and climatic variables including temperatures and precipitation levels contributed to damage severity for pine trees. In contrast, valley depth and soil pH, along with climactic variables of isothermality and temperature contributed to damage severity for hardwood trees. Areas in the mid-south from Mississippi to Alabama, and portions of central Arkansas and Oklahoma showed increased probabilities of more severe levels of tree damage. ConclusionsOur project identified important soil and climatic predictors of tree damage levels, and areas in the southeastern U.S. that are at greater risk of severe wind damage, with management implications under continuing climate change.


Introduction
The southeastern U.S. represents America's "wood basket", producing 60% of the total timber volume of the U.S., and representing 80 million ha of private forest land and 13 million ha of public forest land (Butler and Wear 2013; Oswalt and Smith 2014). This region is subject to signi cant and increasing tornado activity (Ashley and Strader 2016;Dixon et al. 2011), and exposure to both Atlantic and Gulf hurricanes. While these wind events are part of natural disturbance in southeastern forests, they can cause substantive economic and ecological damage to forest stands resulting in the loss of timber, increased reforestation efforts and costs, reduction in carbon sequestration, rapid expansion of invasive plant populations, and a potential increase in forest pest and pathogen populations (Chapman et al. 2008; Gandhi et al. 2007; Marini et al. 2013;McNulty 2002;Vogt et al. 2020). Recovery periods for damaged forest stands to reach structural similarity to undamaged forest stands can take up to 25 years in U.S. temperate forests (Allen et al. 2012). Frequent wind events have major impacts on forest resilience and economic sustainability in commercial forests, and therefore, understanding patterns of events and risks associated with damage from wind events is critical for risk management and planning.
At the stand scale, damage patterns (i.e., the distribution of damage within the storm path) and damage severity (i.e., total area affected or damage magnitude) may vary with the topographic con guration of the landscape. Surface topography affects tornado direction and intensity, with even small changes in uencing near-surface in ow, structure, and path (Lewellen 2012). Steeper slopes have been shown to increase probability of tornado occurrence, likely due to the role of slope in the creation of updraft, although the relationship between elevation and tornadic activity can vary by region and storm characteristics (Frazier et al. 2019;Houser et al. 2020). Remote sensing studies after tornado damage in Georgia and Tennessee forests found more severe damage as tornados descended ridges and diminished damage on the ascent, and this pattern was more pronounced on shallower than steeper slopes ). Numeric simulations suggest surface roughness and topography can in uence the near surface swirl of tornados, thus modifying the intensity at ground level (Lewellen 2012).
Further, terrain features like slope and aspect can sometimes be more strongly associated with damage than the measures of storm meteorology (Kupfer et al. 2008).
Soils represent another important consideration for tree and stand vulnerability, as substrate conditions, soil moisture, and depth of bedrock can affect tree stability during wind events (Loope et al. 1994;Phillips et al. 2008). Poor soil drainage can lead to shallow root development, potentially increasing risk of windfall, however shallow root systems can also equal a larger proportion of biomass allocated to roots, as well as increased root contact between trees, both of which can increase stability from windthrow (Nicoll and Ray 1996). Different tree species have adaptations for rooting into different types of soil (Wang and Xu 2009). Hence, soils and drainage class can interact with species vulnerability, for example oak (Quercus spp.) trees are more vulnerable to damage in drier soils while pine (Pinus spp.) trees have increased resistance to wind damage in drier soils (Rutledge et al. 2021). Chemical and physical properties of the soil can affect root development and nutrient availability, which in turn in uences tree health and root stability. For example, low soil pH can be related to wind damage vulnerability in forests, possibly due to lower nutrient availability in these soils (Mayer et al. 2005).
Vulnerability at the stand and individual tree scales is one facet of an overall understanding of risk, but can be limited in terms of broader spatial application (Kupfer et al. 2008). Risk must take into account likelihood of an event along with expected damage magnitude (Hanewinkel et al. 2011). This estimation can be accomplished through integration of probability models with spatial data, which can provide highresolution disturbance risk information over broad geographic areas (Suvanto et al. 2019). Landscapescale patterns of damage and the relationships of those patterns to topography and soils are needed to ll gaps in our understanding of wind disturbance in southeastern forests . Understanding these patterns of periodic wind disturbances can allow forest systems to be managed prior to disturbances in ways that alter how the system would respond based on its particular vulnerabilities (Kupfer et al. 2008). In this way, management for wind damage can take place within a larger context of risk management, and integrating understanding of likelihood and severity into the broader aims of stand management (Mitchell 2013).
Our research objectives were to: 1) determine the landscape-scale patterns of event probability and damage across the region; 2) determine the topographic, climactic and soil variables which separate the probabilities of different damage classes at the landscape scale; and 3) develop predictive models of wind damage severity in forests over the southeastern U.S. We used post-damage assessments conducted during 2012-2020 across the southeastern region by the National Oceanic and Atmospheric Administration (NOAA). Topographic, soil and climate data were utilized as predictor variables, and event probability was integrated to develop risk maps. This information can be integrated with on-the-ground vulnerability analyses based on species and stand structure characteristics at the local scale for a more complete understanding of risk for the individual landowner. For the purposes of this paper, we focused our attention on tornados and severe thunderstorms because of the availability of NOAA's post-damage assessment data, which allowed for a regional scale analysis of damage risk related to these speci c events. However similar predictive modelling processes could be applied to hurricanes and other wind events if region-wide post-event damage data were readily available.

Data collection and preparation
To develop a regional scale spatial context of probability and magnitude, wind event history data were obtained from the NOAA, National Weather Service (NWS) Storm Prediction Center (https://www.spc.noaa.gov/gis/svrgis/). Tornado path data for 1950-2019 includes start and end points, width of track, magnitude on the F scale before 2007, and the EF scale after 2007. Line shape les for tornado events were rasterized to 1 km 2 resolution using R statistical software package "raster" (Hijmans and Etten 2012). The probabilities of tornadic events within each pixel were calculated based on numbers of events divided by the number of years of data. Average magnitude per pixel was calculated based on the minimum 3-second gust associated with the EF scale for tornados 2007 onwards. Magnitude was only considered for 2007 onwards due to the change in magnitude scale applied in 2007 by NOAA from the F to EF scale.
Tree damage data for the southeastern U.S. was gathered from NOAA post-damage assessments, obtained via the damage assessment viewer website (https://apps.dat.noaa.gov/StormDamage/DamageViewer/) for 2012-2020. Only tree damage was included for these analyses (damage types 27 and 28, representing "hardwood" and "softwood"). Damage points were rasterized by maximum damage category and maximum damaging wind speed to 1 km 2 pixel resolution. These data yielded a total of 28,449 tree damage points related to tornados and 6,541 tree damage points related to thunderstorm and tropical storm winds, including 19,450 hardwood tree damage observations and 15,562 softwood (or thereafter pine tree) damage observations. Independent predictor variables fell into three main categories: terrain, soil, and climate. Terrain attributes were derived from the Shuttle Radar Topography Mission (SRTM) 30 Arc-Second elevation map, using SAGA GIS terrain modules. The attributes include aspect, channel distance, convergence index, diurnal anisotropic heating, LS factor (slope and length steepness factor), slope, slope position, terrain ruggedness index, topographic wetness index, total catchment area, valley depth, and wind exposition index. In addition, distance from coast was calculated for each 1 km 2 pixel. Soil data were collected from the USDA State Soil Geographic Database (STATSGO) Database including available water storage capacity in the top 250 cm of soil, a surface soil erodibility factor, depth of bedrock, pH of soil top layer, permeability and bulk density, and rasterized to 1 km 2 resolution. Historic climate data values averaged for the period between 1970 and 2000 were obtained from Worldclim.org, and included mean, minimum and maximum annual temperatures, average diurnal range (mean of monthly max temp -min temp), Isothermality [(Diurnal Range / Annual Range) * 100], temperature seasonality (standard deviation x 100), as well as mean temperatures of the wettest, driest, warmest, and coldest quarters each year.
Precipitation variables included annual precipitation, seasonal precipitation (coe cient of variation), as well as precipitation of the wettest and driest month, and precipitation of the wettest, driest, warmest, and coldest quarters each year. Tree damage is categorized by NOAA based on a score of 1-5, with score of 1 representing small limbs broken (up to 1" diameter, or ~2.5cm), 2 indicating large branches broken (1-3" diameter or ~2.5-7.6 cm), 3 representing trees uprooted, 4 indicating trunk breakage and 5 indicating trees debarked with only stubs of largest branches remaining (McDonald and Mehta 2004). For the purposes of our analysis, we merged these ve categories into three categories as follows: 1) branch damage -NOAA category 1 and 2 together; 2) uprooting -NOAA category 3 alone; and 3) trunk breakage -NOAA categories 4 and 5 together. It is notable that category 5 is an extremely rare event, and only eight data points in the data set of 28,449 represented this damage.

Statistical analyses
The analyses were conducted in two steps as follows: 1) we rst selected the most important predictors that could discriminate between consecutive damage classes; and 2) once predictors were found, we tted a multinomial model that separated all classes as a function of the predictors found in the rst step. To select the variables associated with damage categories, we utilized a Generalized Additive Model (GAM), to avoid making parametric assumptions about the form of the response. The GAM framework, as implemented in the 'gamsel' package in R (Chouldechova et al. 2018), sequentially penalizes a likelihood function to automatically select the most important linear and non-linear predictors. This framework allows each variable to be estimated as null, linear, or low-complexity curve in a nonparametric manner as given by the data. GAMSEL utilizes a LASSO (Least Absolute Shrinkage and Selection Operator) type of regularization method that penalizes regression coe cients with poor predictive powers, shrinking them towards zero while nding the strength of the regularization term sequentially using 10-fold cross validation against a training set of 80% of the data to select the optimal model. This process reduces over tting and removes extraneous and highly correlated variables. The response variable was assumed to be binomially distributed for each damage level, beginning with branch damage set to zero and both uprooting and trunk breakage levels set to one, followed by the branch and uprooting damage categories set to zero and the trunk breakage category set to one. The GAMSEL process was applied to in this manner to determine which variables were separating the lower vs. higher levels of damage. A tuning parameter (gamma) was set to 0.5 to penalize linear and non-linear components equally, and each variable were allowed 10 basis functions and a maximum of seven degrees of freedom to allow some exibility while limiting over tting of non-linear forms. Variables selected via this process were then utilized in a second step using a multinomial Generalized Additive Model (GAM) as implemented in the "mgcv" package in R (Wood 2017). Smoothing splines were applied to estimate the non-parametric functions for all selected variables except for windspeed, for which a tensor product smooth was applied. For a multinomial GAM, categories must be coded as 0 to K, so that there are K+1 categories and K linear predictors. In our model, branch damage was coded 0, uprooting damage was coded 1, and trunk breakage was coded 2, thus K = 2 with three categories and two linear predictors. Each predictor depends on smooth functions of predictor variables, in this case the variables selected via the GAMSEL as separating low damage category from the higher damage categories (uprooting or trunk breakage), and the variables separating the uprooting category from the trunk breakage category. Predictions were then generated across the region using the selected predictor variables for terrain, soil, and climate, and using average windspeeds from historical NOAA tornado path data. Errors were calculated from the logit values, and predictions on the response scale (probability) were generated. Response scale predictions were then multiplied by event probability to give an overall probability of damage at each level (low, uprooting, or trunk breakage). Finally, future predictions were calculated by multiplying the probabilities across 10, 20, and 30 year intervals, to capture the increasing probability of events over time as follows: Where p = the immediate event probability and t = the amount of time in years. This process was applied for pine and hardwood tree damage separately (Fig. 1).

Results
Landscape-scale patterns of wind event probability

{ [ ]}
The probability of tornadic events based on the 70-year history of data revealed hotspots of tornadic activity in the southeastern U.S. Parts of Oklahoma and Arkansas, along with the mid-south from southern Mississippi through northern Alabama showed a higher probability of tornadoes (Fig. 2a). When examining the spatial pattern of tornado magnitude (Fig. 2b), the mid-south was prominent again as an area of not only high probability, but also higher magnitude tornadoes, particularly in Alabama. Areas closer to the coast, even those which experience frequent tornados such as parts of Florida and coastal Texas, had on average much lower magnitudes than inland areas.

Variables discriminating between damage classes
For pine trees, seven variables were identi ed as predictors of damage beyond broken branches (including uprooting or trunk breakage  (Fig. 3). Windspeed was the primary variable separating trunk breakage and uprooting as higher windspeeds increased the chance of trunk breakage, while a high soil erodibility factor decreased the chance of trunk breakage, making uprooting a more likely scenario in highly eroded soils (Fig. 4). For hardwood trees, seven variables were identi ed as predictors of damage beyond broken branches (including uprooting or trunk breakage). Those variables included windspeed, valley depth, soil pH, isothermality, average temperature of the driest quarter, precipitation in the coldest quarter, and longitude.
The variables which separated trunk breakage from uprooting included windspeed, distance from the coast, soil erodibility factor, and longitude. Precipitation in the coldest quarter and isothermality were not signi cant predictors of damage level in the nal model (Table 1), all other variables were signi cant predictors. The model explained 60.9% of the deviance. The likelihood of hardwoods experiencing damage categories of uprooting and/or trunk breakage increased with increasing windspeed, and also when average temperatures in the driest quarter were greater than 10°C. Damage category tended to decrease with increasing valley depth, and decreased when soil pH was >6.5 (Fig. 5). Windspeed was the primary variable separating trunk breakage and uprooting, with higher windspeeds increasing the chance of trunk breakage, while soil erodibility was not signi cant in separating trunk breakage and uprooting for hardwoods as it was for pine trees (Fig. 6).

Predictive models of wind damage in forests
Predictions of different levels of damage showed broad spatial patterns, with very similar patterns for both pine and hardwood trees (Fig. 7, S1, S2). Branch damage was more common in northwest Texas and in areas closer to the coast, where tornadoes tend to be lower in magnitude. Uprooting damage showed hotspots in southern Mississippi and mid Arkansas, with scattered patches of higher probability through much of the mid-south and Oklahoma. The areas with highest probabilities of trunk breakage aligned with the areas of areas associated with highest event probabilities and highest magnitude tornadic activity, primarily the mid-south through Mississippi and Alabama, and parts of Arkansas and Oklahoma. When projected forward to 10, 20, and 30-year probabilities, parts of the mid-south, mid-Arkansas and Oklahoma exceeded a 60% probability of either uprooting or trunk breakage for pine trees in within 30 years (Fig. 8). For hardwood trees, a portion of southern Mississippi shows >50% probability of uprooting within 30 years, along with scattered patches in Alabama, Arkansas, and Mississippi, while hardwood trunk breakage shows >50% probability of within 30 years concentrated in a portion of northern Alabama and north-central Arkansas, areas where tornadoes have also been historically higher in magnitude (Fig. 9). As a general trend, trunk breakage appears more likely for pine trees and uprooting more likely for hardwoods, but this general trend depended on location and storm intensity.

Discussion
Our study provides the rst regional assessment of tree damage risk in the southeastern U.S. based on tornadic history of the region, damage assessments, and landscape attributes. We found that terrain characteristics of slope and valley depth can contribute to damage severity, along with soil characteristics such as erodibility factors, soil pH, and bedrock depth. In addition, climatic variables including average wind gust speed, average temperatures, and precipitation can also in uence tree damage levels. We found that areas in the mid-south through Mississippi and Alabama, and portions of central Arkansas and Oklahoma had increased probabilities of more severe levels of tree damage (uprooting and trunk breakage) which has implications for forest restoration and resilience activities in these states.
Climatic, landscape, and stand level predictor variables were variously important and different for pine and hardwood trees. Precipitation, in particular in the coldest quarter, was a selected predictor for both pine and hardwood tree damage, but interestingly showed opposite patterns for tree-types with increased risk for pine and decreased risk (though non-signi cant) for hardwood trees. This is consistent with observations that soil characteristics can interact with tree species vulnerability, as hardwood trees are often more vulnerable to wind damage in drier soils, while pine trees are more vulnerable to wind damage in wetter soils (Rutledge et al. 2021). For hardwood tree damage, higher temperatures during the dry quarter were associated with increased damage severity, again suggesting that dry soils increased damage risk. Soil pH was related to higher probabilities of uprooting for hardwoods in soils with higher acidity (or lower pH). Similarly, studies in Europe have found soil pH to be a signi cant factor related to damage after storm events, with greater damage in areas of higher soil acidity, although the reason for this relationship is poorly understood (Mayer et al. 2005). Bedrock has also been noted as a factor in tree damage, as shallower bedrock often facilitates stability for trees which are able to penetrate the bedrock, reducing uprooting potential, though sometimes increasing trunk breakage potential (Foster and Boose 1992;Loope et al. 1994). In another instance, 93% of uprooted trees (mostly pines) after a tornado in Arkansas were found to have penetrated the bedrock (Phillips et al. 2008). Hence, this effect can be variable based on tree species and the type of bedrock, and the capacity for very deep rooting in deep soils can have the same effects (Peterson 2000a). In our study, deeper bedrock increased damage risk for pine and not for hardwood trees, and it could be due to different rooting strategies for these tree-types.
Slope had a very moderate in uence on the damage category and for pine trees only, there was increased damage on moderate slopes of 10° to 30°. Previous research has indicated that greater damage occurred as tornadoes descended slopes and less damage as they ascended slopes ). This could be part of the reason slope does not stand out in this analysis, as the damage level would depend on the slope's relation to the path of the tornado. While most tornados in this region have a general southwest to northeast direction (Thom 1963), aspect was also not an important factor in predicting tree damage. Other studies have noted variability between these factors and storm intensity based on location and storm characteristics (Frazier et al. 2019;Houser et al. 2020). It is possible that the slope/aspect relationship deserves more focused attention in localized risk assessments, and a broad regional approach such as our study may be insu cient for elucidating a clear relationship between those variables and damage risk.
Our models explained 53% and 61% of the deviance for pine and hardwood trees, respectively. This likely re ects the fact that many local variables come into play when explaining damage variation, and thus landscape level variables will generally be limited in their ability to account for substantial deviance. Tree characteristics including species, age, diameter, height, crown size, and crown density often explain ) and stand structure. For example, mechanistic models have shown that irregular spaced stands tend to experience scattered damage, where regular stands can be resilient to lower wind speeds and experienced collapse above a critical wind speed (Ancelin et al. 2004). Nonetheless, we observed some broader spatial patterns which can be useful to land managers, and they can be integrated into their local risk analysis considering tree and stand structure characteristics.
The frequency and strength of tornadoes through the mid-south area from Mississippi through Alabama has been noted in meteorological literature, and the area has been dubbed "Dixie Alley" ( , it stands to reason that this mid-south region would experience higher risks of tree damage given both higher frequency and magnitude of tornadic events along with signi cant area of forest land. There were only minor differences between the spatial patterns of predicted damage to pines vs. hardwoods, but the differences between spatial patterns of the type of damage were striking. Higher damage categories for both hardwood and pine trees were observed through the mid-south, with uprooting hotspots appearing in southern Mississippi and central Arkansas, and trunk breakage hotspots extending into northern Alabama and central Oklahoma. Hence, differential management of stands under wind disturbances as based on location may be a future consideration. For the most part, the immediate damage probability did not exceed 5% with the exception of uprooting probability for hardwood trees (highest overall probability 6%), and broken trunks for pine trees (highest overall probability 7%). However, a typical stand rotation of 25-30 years would experience successively increasing risk over time, and the longer a stand has gone without a disturbance event, the more likely a disturbance event will occur. It is hence, useful for land managers to know their relative risk with time, which will not be equal in all locations. Areas which show a high relative risk within a 30 year time frame may consider damage mitigation within their risk management plan for each rotation. Understanding the likely types of damage can also be helpful for risk mitigation, as salvage operations for uprooted trees vs. broken trunks may have differential costs associated with them.
Our study's projections are only relative to tornado or severe thunderstorm damage, and do not include hurricanes and tropical storms. Areas like south Florida, costal Texas, and coastal Louisiana may be at low risk for tornado damage, however these areas are likely to have higher risk from hurricane and tropical storm damage. Future risk assessments will incorporate hurricane damage risk for a more complete picture of windthrow tree damage throughout the southeastern U.S. We chose to focus rst on tornados whose patterns are less predictable, and tend to do more damage to inland forests. In addition, it is important to note that these damage probabilities indicate the likely level of damage (broken branches, uprooting, or trunk breakage) but do not predict the percent of treefall within a stand. This is because the current EF scale post-damage assessments do not account for percent of treefall in an area, only the type of damage. Satellite and remote sensing methods are now making it possible for more precise assessments of treefall damage for EF scale estimates, even in remote areas, rather than on the ground observations (Godfrey and Peterson 2017). As this more precise and inclusive data are being gathered over time, future models will increase the precision of the current estimates.
A limitation of the current study is the lack of detail relating to tree characteristics associated with the NOAA damage data set beyond "hardwood" or "pine" trees. A more precise model would include tree attributes such as age class, height, diameter, and species. Future work may re ne these models and improve the deviance explained by incorporating these variables if su cient data are collected or becomes available. The variable selection process for determining predictors of damage utilized in these analyses can also be applied to other types of damaging events when post-storm damage data are available for such predictive purposes. The purpose of these analyses was to provide landowners and managers some insight as to the relative risk of tree damage given their location, recognizing that local factors and stand characteristics will also in uence their risk pro le. Tornadoes are di cult events to predict with any precision, and the occurrence of a tornado or severe thunderstorm event does not inherently mean damage will be severe. Yet, we noted some general spatial patterns in the type and severity of damage that are not random and give some insight into relative risk associated with these catastrophic events in southeastern U.S. forests.

Statements and Declarations
Declarations Figure 1 Model process work ow for our study on assessing wind risks in the southeastern U.S. forests. GAMSEL = Generalized additive model selection GAM = generalized additive model, EF = measure of tornado magnitude.  Variables separating branch damage from higher damage categories (uprooting or trunk breakage) for pine trees in the southeastern U.S., including the average 3-second gust windspeed associated with tornado or severe thunderstorm (a), slope (b), valley depth (c), bedrock depth (d), average temperatures in the warmest quarter of the year (e), precipitation in the coldest quarter of the year (f), and longitude indicating the partial effects of along the longitudinal axis (g). Grey area represents standard errors. Y axis indicates the partial effect of the predictor variable on damage level. Y axis label indicates: smoothing function (predictor variable, effective degrees of freedom). Effective degrees of freedom = 1 indicates a linear relationship, higher effective degrees of freedom represent a non-linear relationship.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.