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

T.Adam Coates·W.Mark Ford

Abstract 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.

Keywords Eastern hemlock·Hemlock woolly adelgid·Mountain laurel·Rhododendron

Abbreviations

GSMNP Great Smoky Mountains National Park

DBH Diameter at breast height

Introduction

Hemlock wooly adelgid (Adelges tsugaeAnnand) 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 thatA.tsugaehas spread approximately 9-20 km a-1since.A.tsugaewas first identified within Great Smoky Mountains National Park (GSMNP) along the North Carolina and Tennessee border in 2002 (Lambdin et al.2006) and many affectedT.canadensisdied rapidly，often within only a few years following infestation (Ford et al.2012).As a foundational species，particularly in riparian zones，T.canadensisloss 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 mostT.canadensismoist 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 latifoliaL.) and great rhododendron (Rhododendron maximumL.)，thereby potentially increasing flammability in a given location under dry conditions (Brose 2016，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 (Quercusspp.)，hickories (Caryaspp.)，and upland pines (shortleaf pine:Pinus echinataMill.;pitch pine:Pinus rigidaMill.;Table Mountain pine:Pinus pungensLamb.)，have experienced expansion of shade tolerant，mesophytic species such as red maple (Acer rubrumL.)，American beech (Fagus grandifoliaEhrh.)，blackgum (Nyssa sylvaticaMarsh.)，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 2016，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 whereT.canadensisoccurs or had occurred prior to mortality fromA.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.

Materials and Methods

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.2 015;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.

Fig.1 Location of Great Smoky Mountains National 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

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).

2003-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 (m2ha-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.Thermeasurement was used to determine the area of each individual shrub within the 0.01 ha fixed area plot (i.e.，area=πr2).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.(2019，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 refers to how each individual fuel-size class responds to changes in relative humidity (Brown 1974).

Fig.2 Orientation of overstory tree，ericaceous shrub，and coarse woody debris transects at each sampling location within Great Smoky Mountains National Park 2019

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，2020).Thus，plot averages for each variable were based upon nine individual measurements.

2019 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 deadT.canadensisbasal 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 asy=year+forest class+elevationusing 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 deadT.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-hocdifferences were examined using LSMEANS separation tests in PROC GENMOD.

Results

Year

DeadT.canadensisbasal area increased from 2004-2019(p＜0.0001) (Table 1，Fig.3).The sum of 1-，10-，and 100-h woody fuel loads also increased (p=0.0494) (Table 1，Fig.4).Woody fuel depth decreased from 2004-2019(p＜0.0001) (Table 2，Fig.5).Total overstory basal area(Fig.3)，liveT.canadensisbasal area (Fig.3)，1000-h fuel loads and total woody fuel loads (Fig.4)，O Horizon thickness (Fig.5)，and ericaceous shrub percent cover (Fig.6) did not change significantly (Table 1).

Forest Class

DeadT.canadensisbasal area was significantly less on mesic sites than xeric or intermediate sites (p＜0.0001)(Tables 1，2，Fig.3).O Horizon thickness was significantly less on mesic sites than xeric sites but did not differ from intermediate sites (p=0.0274) (Tables 1，2，Fig.5).Total overstory basal area (Fig.3)，liveT.canadensisbasal area (Fig.3)，woody fuel loads (Fig.4)，woody fuel depth(Fig.5)，and ericaceous shrub percent cover (Fig.6) did not differ by forest class (Tables 1，2).

Table 1 Parameter estimates for vegetative and fuel response variables，Great Smoky Mountains National Park，Tennessee，2003-2004 and 2019

Table 1 (continued)

Elevation

LiveT.canadensisbasal area decreased (p=0.0361) and deadT.canadensisbasal area increased with elevation(p＜0.0001) (Table 1，Fig.3).O Horizon thickness increased with elevation (p＜0.0001) (Table 1，Fig.5).Total basal area (Fig.3)，woody fuel loads (Fig.4)，woody fuel depth(Fig.5)，and ericaceous shrub cover (Fig.6) did not change significantly with changes in elevation (Table 1).

Fig.3 (A) Mean total，(B) live T.canadensis and (C) dead T.canadensis basal area (m2 ha-1) from 2003-2004 to 2019 across site conditions and elevation gradients;colored areas represent 95% confidence interval for each site class

Fig.4 (A) Mean 1-，10-，and 100-h fuel loads，(B) 1000-h fuel loads，and (C) total down and dead woody fuel loads (tonnes ha-1) from 2003-2004 to 2019 across site conditions and elevation gradients;colored areas represent the 95% confidence interval for each site class

Fig.5 (A) Mean woody fuel depth and (B) O Horizon thickness(cm) from 2003-2004 to 2019 across site conditions and elevation gradients;colored areas represent 95% confidence interval for each site class

Fig.6 Mean ericaceous shrub cover (%) from 2003-2004 to 2019 across site conditions and elevation gradients;colored areas represent 95% confidence interval for each site class

Discussion

Vegetative Changes

As hypothesized，deadT.canadensisbasal area increased over the～ 16-year time period (Table 1).This change was most pronounced on high elevation (up to 765 m)，xeric and intermediate locations.LiveT.canadensisdecreased with increasing elevation，as well，but forest class was not significantly related to this decrease (Tables 1，2).In New England，Orwig et al.(2002) observedT.canadensisdecline and found that xeric locations succumbed most rapidly toA.tsugaeinfestation，but slope and elevation were not significantly related to eventualT.canadensisstem mortality.Willeford Bair (2002) determined thatT.canadensisstems with an excellent crown condition declined from 80% in 1990 to near 0% in 2000 in Shenandoah National Park，Virginia.The most significant declines in crown condition occurred above elevations of 750 m in their study.Few plots in our study of GSMNP were present at this elevation.In the central Appalachian Mountains of Virginia and West Virginia，Martin and Goebel (2012)utilized a 9-to 32-year chronosequence to evaluate “timesince-invasion” ofA.tsugaeto assessT.canadensisinfestations across a range of topographic positions.Contrary to Willeford Bair (2002)，Martin and Goebel (2012) found thatT.canadensisdecline 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 ofT.canadensisdeclined with increasing time since infestation.These authors hypothesized that the broad loss of allT.canadensisstems，regardless of stem crown class，might lead to vegetative increases of non-native species，such as tree-of-heaven (Ailanthus altissimaMill.)，and some native species，such asF.grandifoliaandR.maximum，due to their ability to persist in semi-shaded conditions for long periods of time asT.canadensisstems 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 increasingT.canadensismortality.Ford et al.(2012) found that understory shrub growth increased 2.6 times followingT.canadensismortality at the CoweetaHydrologic 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.latifoliawas most common on dry，southwestern slopes andR.maximumon moist，northeastern slopes.However，only 25% of the plots measured in that study containedR.maximumand only 42% hadK.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.canadensisbasal 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 deadT.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，includingT.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 prolonged fall rates affected by cause of mortality，climate，disturbance regime，soil type，species，and stem size (Morrison 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 theA.tsugaeinfestation that was actively occurring within GSMNP in 2003-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，duffthickness，1000-h woody fuel loading，and ground cover composition (includingR.maximum).The 1-，10-，and 100-h woody fuel loads in their study did not differ between site designations，litter depth，orK.latifoliacover.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 fromA.tsugaeor 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 wildf rie event (Kuligowski et al.2020).The prevalence of standing and deadT.canadensisstems 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.(2007，2009) identified arthropod species withinT.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 withT.canadensiscould decline followingA.tsugaeinfestation.Therefore，impacts to these taxa could have substantial，broader community integrity impacts.Continued monitoring ofT.canadensischanges 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 deadT.canadensisstems 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.

AcknowledgementsThe 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 AccessThis 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://creativecommons.org/licenses/by/4.0/.