Fuel and vegetation changes in southwestern, unburned portions of Great Smoky Mountains National Park, USA, 2003–2019

Overstory basal area, ericaceous shrub cover (Kalmia latifolia L. and Rhododendron maximum L.), and fuels (i.e., woody fuel loads and depths and O Horizon thickness) were assessed within Great Smoky Mountains National Park, USA, in 2003 − 2004. Due to recent wildfire activity within the southern Appalachian Mountain region (including Great Smoky Mountains National Park), the potential spread and expansion of ericaceous shrubs, and the impacts of the hemlock woolly adelgid (Adelges tsugae Annand) on eastern hemlock (Tsuga canadensis (L.) Carrière), these same ecosystem components were again assessed in 2019. Elevation and moisture regime (xeric, intermediate, and mesic) were included in this assessment as potential influential factors. An evaluation of repeated measurements from 40 plots suggested that O Horizon thickness did not change significantly over the 16-year period, but increased as elevation increased, and moisture regime (xeric O Horizon thickness > mesic O Horizon thickness) was a significant, related factor. The sum of 1-, 10-, and 100-h fuel loads (fuels less < 7.6 cm diameter) increased, whereas woody fuel depth decreased over the 16-year period. No significant changes in 1000-h fuel loads (> 7.6 cm diameter), total woody fuel loads, ericaceous shrub cover, total basal area, or live T. canadensis basal area were observed. Live T. canadensis basal area decreased with increasing elevation. Dead, standing T. canadensis basal area increased from 2003–2019, and that increase was most pronounced as elevation increased on xeric and intermediate sites. Overall, we found that: 1. hypothesized increases in total woody fuel loads and ericaceous shrub cover were not present; and 2. elevation and moisture regime were most related to observed changes in vegetation and fuel condition.


Introduction
Hemlock wooly adelgid (Adelges tsugae Annand) is a nonnative, invasive insect that has affected eastern hemlock (Tsuga canadensis (L.) Carrière) in many locations of the eastern US since it was first identified in Virginia in 1951 (Limbu et al. 2018). Through 1985, this pest was largely considered a nuisance only to ornamental trees in urban landscapes. Morin et al. (2009) reported that A. tsugae has spread approximately 9 − 20 km a −1 since. A. tsugae was first identified within Great Smoky Mountains National Park (GSMNP) along the North Carolina and Tennessee border in 2002 (Lambdin et al. 2006) and many affected T. canadensis died rapidly, often within only a few years 1 3 following infestation (Ford et al. 2012). As a foundational species, particularly in riparian zones, T. canadensis loss has cascading effects on a myriad of ecosystem properties and species (Abella 2014;Coots et al. 2015). Overstory mortality resulting from the impact of invasive insects can alter dead woody fuel properties, loads, and structure. Although not particularly problematic for increased wildfire hazard on most T. canadensis moist riparian and "cove" sites in the eastern US, increased tree mortality may foster an accelerated expansion of ericaceous shrub species (Krapfl et al. 2012), such as mountain laurel (Kalmia latifolia L.) and great rhododendron (Rhododendron maximum L.), thereby potentially increasing flammability in a given location under dry conditions (Brose 2016(Brose , 2017.
Like much of the public land in the southern Appalachian Mountains, active fire suppression and exclusion guided GSMNP's management policies for much of the twentieth century (Cohen et al. 2007). In GSMNP, drier ridges and slopes at low to mid-elevations (~ 1,200 m), once dominated by American chestnut (Castanea dentata (Marsh.) Borkh.), upland oaks (Quercus spp.), hickories (Carya spp.), and upland pines (shortleaf pine: Pinus echinata Mill.; pitch pine: Pinus rigida Mill.; Table Mountain pine: Pinus pungens Lamb.), have experienced expansion of shade tolerant, mesophytic species such as red maple (Acer rubrum L.), American beech (Fagus grandifolia Ehrh.), blackgum (Nyssa sylvatica Marsh.), and other hardwood species (Jenkins 2007). In fire's absence and the cessation of forest harvesting, the latter species have gradually occupied more space and utilized more resources in historically fire-prone areas (Ryan et al. 2013;Kreye et al. 2018;Alexander et al. 2021). As such, over time, "mesophication" (Nowacki and Abrams 2008) increased, creating conditions less conducive to oak, hickory, and pine regeneration and release (Brose 2016). Ericaceous shrub species have expanded simultaneously and have contributed to these conditions (Royo and Carson 2006;Brose 2016Brose , 2017. The managerial, vegetative, and fuel changes described here have been widespread throughout the eastern US. However, these changes are particularly unique and potentially problematic at GSMNP due to the Park's unique attributes, ecological value, and abundance of sensitive communities and endemic or rare species. Fuel changes after a long history of fire suppression may pose unique challenges and threats within GSMNP, as evidenced by the expansive autumn 2016 wildfires that affected large portions of GSMNP and its surrounding communities (Kuligowski et al. 2020). Recent projections suggest wildfire incidence within the southern Appalachian Mountains will increase for at least the next half century (Prestemon et al. 2016) as constraints to implementing many fuel management strategies have aided fuel load increases and fuel structure changes that may enhance potential wildfire ignitions and spread (Jolly et al. 2015;Kuligowski et al. 2020). Some studies have suggested that fuel depths, thicknesses, and loads may additionally be influenced by factors such as forest class and elevation within this region, predisposing some locations to potentially hazardous wildland fuel conditions (Waldrop et al. 2007;Stottlemyer et al. 2009).
Considering these unique changes and potential implications, we sought to resurvey forest conditions in 40 plots in 2019 that were initially established in 2003-2004 within GSMNP (Waldrop et al. 2007). The original research plots were established to assess coarse woody debris loads and forest composition, and were included in a larger, 1000-plot, comprehensive study of these ecosystem properties throughout the southern US (Waldrop et al. 2007). Generally, onehour fuel loads (0.00 − 0.64 cm diameter) were significantly greater on plots that had been affected by invasive insects than in undisturbed plots. Larger woody fuels (> 7.62 cm diameter) were greater in plots that had been affected by multiple disturbances prior to 2003 − 2004, including insects, fire, and wind. In 2019, plots were revisited to evaluate woody fuel depths and loads, O Horizon thickness, basal area, and ericaceous shrub cover to determine: (1) how forest overstory composition may have changed from 2003 − 2019; (2) how ericaceous shrub cover may have changed from 2003-2019; and, (3) how woody fuel loads and depths and O Horizon thickness may have changed from 2003-2019. We speculated that live basal area increased in all plots and dead basal area increased in plots where T. canadensis occurs or had occurred prior to mortality from A. tsugae. We also predicted that ericaceous shrub cover increased over time. Due to insect-caused mortality and in-growth, we predicted that fuel depths, thicknesses, and loads, except for O Horizon thickness, also increased.

Study Location
Great Smoky Mountains National Park (GSMNP), within the Blue Ridge Physiographic Province portion of the southern Appalachians Mountains, was established in 1934 along the North Carolina and Tennessee border (Fig. 1). Widely recognized for its unique ecological features, GSMNP represents one of the most diverse temperate forests in the world and constitutes one of the largest blocks of core forest in eastern North America (Holzmueller et al. 2006). Most forests and other habitat types, (i.e., high elevation, grassy balds of thick grass/shrub vegetation) in GSMNP were once heavily influenced by both natural and anthropogenic wildland fires (Stambaugh et al. 2015;Lafon et al. 2017). At the height of European settlement through the exploitative logging period prior to Park designation, Harmon (1982) estimated the overall mean fire return interval (MFRI) for GSMNP from 1856 − 1940 as 12.7 years. Fire presence was highly variable depending on forest community and elevation, however, and Flatley et al. (2013) estimated MFRI to be 4 − 7 years in some GSMNP locations.
Our research in GSMNP was conducted in the southwestern portion of the Park near Townsend, Tennessee (Fig. 1). One hundred fifty plots were established in 2003 (Waldrop et al. 2007). Forty plots were identified to potentially revisit in 2019 based upon accessibility, personnel, time, and finances. Specifically, to our knowledge, these plots had not been subjected to any additional anthropogenic disturbances or wildfires from 2003 − 2019, including the 2016 wildfires (R. Klein, National Park Service, personal communication).
The elevation range for the 40 plots was 278-765 m. The area is underlain by sedimentary geology of the Cades, Shields and Wilhite formations (Neuman and Nelson 1965) and soils are Ultisols or Inceptisols (USDA NRCS 2009). Forests were mostly fully mature, second growth, temperate deciduous, pine or mixed types (depending on elevation and aspect) that established following widespread harvesting prior to Park formation in the 1930s (Whittaker 1956). Landscape ecosystem classification (LEC), as measured by Stottlemyer et al. (2009), was used to differentiate potential fuel loads and properties as a function of adjacent topography, slope steepness, and topographic shape (Hutto et al. 1999;Carter et al. 2000;Abella et al. 2003). Four LEC distinctions were noted for the GSMNP sampling locations: mesic, submesic, intermediate, and xeric. For statistical analyses, submesic and intermediate plots were combined into one category (hereafter, intermediate).

− 2004 Sampling
During this time period, overstory trees, ericaceous shrub cover (% cover within a 0.01 ha plot), woody fuel loads (tonnes ha −1 ) and depths (cm), and O Horizon thickness (litter and duff, cm) were measured at each of the original plots. Overstory trees > 2 m tall were identified by species, and diameter at breast height (DBH) was measured for each tree within a 0.02 ha fixed area plot (Fig. 2). Total, live, and dead basal area (m 2 ha −1 ) were calculated based upon these inventories. Ericaceous shrub cover was determined using individual shrub crown measurements assessed within half of the 0.02 ha fixed area plot (or 0.01 ha) (Fig. 2). Shrub crown diameter was measured twice for each ericaceous shrub using a 50-m measuring tape due to the inherent lack of completely circular crowns. These measurements were taken approximately perpendicular to one another. Mean crown diameters were calculated and halved to obtain a mean radius, r. The r measurement was used to determine the area of each individual shrub within the 0.01 ha fixed area plot (i.e., area = r 2 ). The sum of the shrub crown areas in each plot was then divided by the total area of the fixed area plot to determine ericaceous shrub percent cover.
Brown's Planar Intercept Method (Brown 1974), as modified by Coates et al. (2019Coates et al. ( , 2020, was used to estimate woody debris loads. Using this technique, down and dead woody fuels 0 − 0.64 cm, 0.65 − 2.54 cm, 2.55 − 7.62 cm, and > 7.62 cm diameter were tallied by 1-, 10-, 100-, and 1000-h time-lag size classes, respectively, along three, 15.2 m transects established at a 45° angle (Fig. 2). Time-lag Park, North Carolina and Tennessee, USA, and the specific study area located in the western portion of the Park where research plots were established in 2003 (dots). Map created using ESRI ArcGIS refers to how each individual fuel-size class responds to changes in relative humidity (Brown 1974).
O Horizon thickness (cm) and woody fuel depth [(cm; defined as the distance from the top of the litter layer to the top of any woody debris crossing the transect) were measured at three locations along each 15.2 m transect (at 3.7, 7.3 and 12.2 m) (Coates et al. 2019(Coates et al. , 2020. Thus, plot averages for each variable were based upon nine individual measurements.

Sampling
The methods for the 2019 overstory tree, O Horizon thickness, woody fuel load, and woody fuel depth sampling were the same as those used in the 2003 − 2004 sampling. Ericaceous shrub cover estimates differed in 2019, however. Due to time constraints, individual ericaceous shrub crowns were not measured to determine percent cover. Instead, three field technicians visually estimated this percentage within the 0.01 ha fixed area plots and an average of those estimates was derived to represent percent cover. Permanent rebar markers were used at the center point of fuel transects in 2003 − 2004; however no markers were placed at the ends of the sampling transects, therefore the two sampling periods could not be evaluated as a true, paired design.

Statistical Analyses
To assess changes over time in response variables (vegetative: ericaceous shrub percent cover, total basal area, live and dead T. canadensis basal area; fuel: O Horizon thickness, sum of 1-, 10-, 100-h fuel loads, 1000-h fuel loads, total fuel loads, woody fuel depth), a series of generalized linear models was carried out for each variable as y = year + forest class + elevation using the PROC GENMOD procedure in SAS 9.4 (SAS Inc., Cary, NC) whereby year and forest class were considered fixed treatments and elevation a continuous covariate. Using the PROC FFM procedure, it was determined that only O Horizon thickness was normally distributed. For ericaceous shrub percent cover and basal area of dead T. canadensis, a zero-inflated Poisson distribution was used. For all other response variables, a lognormal distribution function was used in the model. Where statistically significant differences in class treatment were observed (α = 0.05), post-hoc differences were examined using LSMEANS separation tests in PROC GENMOD.

Vegetative Changes
As hypothesized, dead T. canadensis basal area increased over the ~ 16-year time period (Table 1) (2012) utilized a 9-to 32-year chronosequence to evaluate "timesince-invasion" of A. tsugae to assess T. canadensis infestations across a range of topographic positions. Contrary to Willeford Bair (2002), Martin and Goebel (2012) found that T. canadensis decline was slowed by increasing elevation (range 345 − 1057 m) and decreased incident radiation. These dynamics were overshadowed by time since infestation, however, as nearly all stem size classes of T. canadensis declined with increasing time since infestation. These authors hypothesized that the broad loss of all T. canadensis stems, regardless of stem crown class, might lead to vegetative increases of non-native species, such as tree-of-heaven (Ailanthus altissima Mill.), and some native species, such as F. grandifolia and R. maximum, due to their ability to persist in semi-shaded conditions for long periods of time as T. canadensis stems declined and eventually contributed to down and dead woody debris pools.
Prior to remeasurement in 2019, it seemed plausible that ericaceous shrub percent cover would increase with increasing T. canadensis mortality. Ford et al. (2012) found that understory shrub growth increased 2.6 times following T. canadensis mortality at the Coweeta  Hydrologic Laboratory in nearby western North Carolina. However, no significant change in ericaceous shrub percent cover was noted at the specific locations measured within GSMNP (Table 1). This result might be related to the specific locations evaluated in this study and may not have encompassed effects present in other portions of GSMNP. It also does not address potential changes that may have occurred at the specific locations measured in 2003-2004 due to our inability to locate some of the rebar markers installed during the initial sampling period and the differences in cover estimation between the two sampling periods (2003 − 2004: measured crown diameters; 2019: visual estimation).
In their examination of multiple locations throughout the southern Appalachian Mountains, including GSMNP, Waldrop et al. (2007) determined that ericaceous shrub percent cover was affected by topographic position: K. latifolia was most common on dry, southwestern slopes and R. maximum on moist, northeastern slopes. However, only 25% of the plots measured in that study contained R. maximum and only 42% had K. latifolia, far less than was hypothesized for both species for the area. Similarly, both the extent and amount of change for the ericaceous shrubs in our study was less than hypothesized.

Fuel Changes
Changes in down and dead woody fuel loads were not uniform at our survey sites: the sum of 1-10-100 h fuel loads increased over the sampling period but 1000-h fuel loads did not increase (Table 1). Added together, the total down and dead woody fuel load did not change significantly over this period (Table 1). Standing and dead, or snag, T. canadensis basal area increased (Table 1) and these dead stems have not been added to the down and dead woody debris pool. Depending on the size of these stems and the overall decay rate of standing and dead T. canadensis, future woody debris loads may increase in some of the locations we resurveyed at GSMNP. In a study of snag longevity in central Maine for 10 native tree species, including T. canadensis, Garber et al. (2005) found that snag longevity increased with diameter and decreased with stand basal area. Other research suggests initial lag times for snag fall range from 1 to 15 years (McArdle 1931;Keen 1955;Lyon 1977), followed by  rison and Raphael 1993;Bull et al. 1997;Garber et al. 2005). It is plausible that the decrease in woody fuel depth from 2004 − 2019 (Table 1) was related to the initial, immediate increase in branch and limb fall resulting from the A. tsugae infestation that was actively occurring within GSMNP in 2003GSMNP in − 2004 Forest class and elevation were not significantly related to the woody fuel load and depth changes within GSMNP (Tables 1, 2). In their study of fuels and vegetative cover within the Chauga Ridges in northwestern South Carolina, a location with similar soils and elevation (305 − 80 m) as this part of GSMNP, Stottlemyer et al. (2009) determined that xeric, intermediate, submesic, and mesic sites could be distinguished by differences in woody fuel depth, duff thickness, 1000-h woody fuel loading, and ground cover composition (including R. maximum). The 1-, 10-, and 100-h woody fuel loads in their study did not differ between site designations, litter depth, or K. latifolia cover. O Horizon thickness was greatest on xeric sites in the Chauga Ridges, like our results (Tables 1, 2). Higher elevation was also associated with greater O Horizon thickness (Table 1) in our study, like the findings of Waldrop et al. (2007) in their more comprehensive study of the southern Appalachian Mountains.

Management Implications
As our results show and those of others within this region (Waldrop et al. 2007;Stottlemyer et al. 2009), vegetation and fuel dynamics within the southern Appalachian Mountains are complex. Interactions of successional dynamics, sites, moisture regimes, and disturbances are among the factors that can impact changes in these variables over time (Stottlemyer et al. 2009;Limbu et al. 2018). Changes in these factors may influence potential wildfire ignitions, including their intensity and severity, a common concern among land managers in the region. This potential was actualized within GSMNP and many other locations throughout the southern Appalachian Mountains in autumn 2016. At that time, a significant drought and an advancing weather system allowed fire behavior and spread to be amplified from an anthropogenic ignition in the Park (Guthrie et al. 2017). Under those conditions, any fuel load increases that had occurred because of natural disturbances, such as mortality from A. tsugae or other wind and ice events, may have potentially amplified wildfire behavior. Additionally, increases in ericaceous shrubs, even in pockets on specific slopes, could potentially affect fuel loading and fuel structure, increasing the likelihood that crown fires could be generated due to the ladder effect of fuel structure (Whelan 1995;Keane 2015). The wildland-urban interface present in the  surrounding communities of Gatlinburg and Pigeon Forge, Tennessee was greatly impacted by this wildfire event (Kuligowski et al. 2020). The prevalence of standing and dead T. canadensis stems may be an important consideration for future land management planning, along with vegetation that might continue to fill any gaps created by the loss of these stems. Additionally, our study was limited to resampling at 40 locations within GSMNP, therefore the results may not fully encapsulate changes that have taken place throughout GSMNP or the Appalachian Mountains as a whole.
In most forest systems, standing snags and down and dead woody fuels are foundational ecosystem components (Yuan et al. 2017;Campbell et al. 2019). These components support: wildlife habitat and refugia; feeding substrate for invertebrates, fungi, and microorganisms; and nutrient cycling and soil development processes, including carbon sequestration (Abella 2014;Magnússon et al. 2016;Marcot 2017). In eastern hardwood or mixed-pine hardwood forests, most research has focused on the value and use of snags by cavity-nesting avifauna (Kilgo and Vukovich 2014) or bats (Silvis et al. 2016), and that of down and dead woody fuels with woodland salamanders (Ford et al. 2002;Moseley et al. 2008) or small mammals (Menzel et al. 1999;Owens et al. 2008). Baird et al. (2007Baird et al. ( , 2009) identified arthropod species within T. canadensis-dominated areas of GSMNP that had not been previously identified in any other ecological system. Rohr et al. (2009) determined that seven taxa associated with T. canadensis could decline following A. tsugae infestation. Therefore, impacts to these taxa could have substantial, broader community integrity impacts. Continued monitoring of T. canadensis changes within the confines of GSMNP and other locations relative to fire frequency, fire effects, and other ecological components are warranted.

Conclusions
In our study, some vegetative and fuel components changed in the measured areas of GSMNP from 2003 − 2019, however, these changes did not fully conform with our predictions. For example, total woody fuel loads, woody fuel depth, and ericaceous shrub percent cover did not increase within the study area. With additional time, standing and dead T. canadensis stems will fall to the forest floor and add to the total down and dead woody fuel load. As with most natural systems, the ecology of GSMNP is not static (Boyle et al. 2014) -the complex topography, community diversity and myriad of natural disturbance agents means the Park has been and will continue to be a dynamic system that is difficult to model. Moreover, embedded in a humanmodified system, other stressors also affect and/or interact with ecological processes in GSMNP, as is the case with other ecosystems globally. It is important to take account of the changes occurring in these systems at both landscape and site-specific levels to provide managers with the information needed to respond accordingly and meet their desired, management objectives.
Acknowledgements The authors wish to acknowledge Rob Klein for his assistance with field operations in 2019 and for a thoughtful review of an earlier draft of this manuscript. Kevin Hiers also reviewed an earlier version of this manuscript and his comments and suggestions greatly improved the quality of this work. The authors wish to acknowledge John Kese, Brandon Hughes, and Nick Boley for their hard work and perseverance in challenging field conditions during the 2019 resampling. We would also like to acknowledge Helen Mohr for sharing the original 2003−2004 data. Many thanks, as well, to Tom Waldrop, Gregg Chapman, Lucy Brudnak, Chuck Flint, Mitch Smith, and others within the USDA Forest Service Southern Station Center for Forested Disturbance who envisioned and conducted the 2003−2004 sampling and analyses in GSMNP. Lastly, we thank the National Park Service for funding this effort (Task Agreement P19AC01059).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.