Spring and Autumn Phenological Variability across Environmental Gradients of Great Smoky Mountains National Park, USA

Abstract

Mountainous regions experience complex phenological behavior along climatic, vegetational and topographic gradients. In this paper, we use a MODIS time series of the Normalized Difference Vegetation Index (NDVI) to understand the causes of variations in spring and autumn timing from 2000 to 2015, for a landscape renowned for its biological diversity. By filtering for cover type, topography and disturbance history, we achieved an improved understanding of the effects of seasonal weather variation on land surface phenology (LSP). Elevational effects were greatest in spring and were more important than site moisture effects. The spring and autumn NDVI of deciduous forests were found to increase in response to antecedent warm temperatures, with evidence of possible cross-seasonal lag effects, including possible accelerated green-up after cold Januarys and early brown-down following warm springs. Areas that were disturbed by the hemlock woolly adelgid and a severe tornado showed a weaker sensitivity to cross-year temperature and precipitation variation, while low severity wildland fire had no discernable effect. Use of ancillary datasets to filter for disturbance and vegetation type improves our understanding of vegetation’s phenological responsiveness to climate dynamics across complex environmental gradients.

1. Introduction

Remote sensing excels at efficiently tracking environmental change over broad areas, yet its success at downscaling observations made at coarse resolution has been mixed [1,2]. This challenge of scale is greatest when the phenomenon to be monitored is highly complex in type or composition or where the dynamic to be captured changes rapidly, such as with near-real-time disturbance detection or seasonal phenological transitions [3]. Our ability to efficiently and accurately track growing season dynamics is particularly critical as climate and disturbance regimes are changing. The onset of spring and duration of the photosynthetically active period link climate dynamics to annual net primary productivity and carbon flux [4,5], nutrient dynamics [6], and water use [7], while providing an indicator of disturbances [8,9,10,11] and longer-term change [12,13,14,15,16].

Our broad-scale understanding of variation in the timing of spring and autumn heavily relies on satellite remote sensing-based measures of land surface phenology (LSP) [17]. LSP time series efficiently capture coarse resolution vegetational change, but the specific vegetational attributes that are being captured by this technology are not always clear [18]. While field observations strive to capture subtle differences in the response of individual plants, gridded measures of LSP provide a community-level indicator that captures the collective response of individual trees and species that are visible from above. This makes gridded observations inherently difficult to accurately and representatively validate with field data or other remote sensing approaches [19]. What gridded approaches lack in biological precision, they make up for with temporal consistency and spatial coverage; LSP is particularly well-suited for tracking broad-scale questions, such as how vegetation responds to year-to-year temperature and precipitation variation along environmental gradients.

One key challenge for understanding the causes of seasonal variation in LSP across complex topographic gradients is the importance of microclimate. Topographic factors affect how much solar radiation is received at different times of year, and cold air drainage, aspect and elevation affect the risk of frost and summer heat stress [20,21]. While topography is constant across years, the importance of topographic attributes that affect exposure to drought or frost often is not, since the importance of topography may emerge only during seasonal extremes. For example, low elevation forests may be less capable of extending growth during drought, because they are more prone to late growing season moisture deficits than are higher elevation sites [22,23]. While years with warmer springs can enhance green-up and carbon assimilation, exposure to late growing season drought could limit those gains [24]. Understanding these local exposures to seasonal weather variation is critical, as microclimate variability can mediate the phenological response of vegetation.

A second key challenge for understanding seasonal variation in LSP in complex environments is the variable response of species and vegetation to day length and seasonal weather. Variability in spring temperature explains most year-to-year variation in green-up, but the response of individual species and even individuals within a species is often not uniform [16,25]. Among hardwood trees, ring-diffuse genera such as Betula, Prunus, Acer and Liriodendron are generally earlier to achieve budburst and maximum leaf expansion, while semi-ring-porous to ring-porous genera such as Quercus, Carya, Fraxinus and Juglans develop later [26,27,28]. There are also strong differences in responsiveness among evergreen and deciduous components, and the ease of satellite green-up detection is related to the fractional cover of deciduous vegetation [29,30]. In deciduous forests, the earliest indication of green-up may be from the understory, weeks before budburst of canopy trees occurs [31]. Spring growth that is captured by remote sensing may transition from an understory signal to overstory at different times along its progression, depending on the canopy cover and understory attributes [32].

A third problem of monitoring LSP responsiveness to seasonal weather is that LSP is sensitive to disturbance and cover change, not just to temperature and precipitation. This sensitivity is so strong that changes in LSP are systematically monitored to recognize and monitor disturbances and subsequent recovery [3,33,34]. While all disturbances may not disrupt the sensitivity of LSP to temperature or precipitation, many disturbances or land cover changes could confound analyses, as locally disturbed vegetation becomes a less consistent measure of LSP sensitivity [2].

Despite the nuancing effects of microclimate and vegetation, our understanding of how climate drives LSP for temperate latitudes is evolving rapidly from experimental work and broad-scale, longer-term analyses [35,36,37,38]. It is widely recognized that late winter and early spring warmth accelerates budburst, with mixed evidence for cross-seasonal legacy effects [25,39,40,41]. The causes of variation in autumn LSP are less well understood, but it tracks the predictable decline in autumn day length and the less predictable summer and autumn weather [42,43,44]. Based on chamber experiments [37], field observations [6,45], and remote sensing [46], it is thought that drought accelerates leaf senescence, while warmer temperatures delay it. This generalization becomes problematic when droughts are warm, or when wet autumns experience cooler than average temperatures. Such warm-dry, cool-wet weather combinations occur frequently due to meteorological feedbacks related to cloudiness. Warmer autumn temperatures increase evapotranspiration rates and drought stress, while wet autumns often come with increased cloud cover which can reduce warmth from solar radiation. These complications confound our understanding of how autumn temperature and precipitation influence senescence.

In this paper we evaluate the contributions of topography, forest cover type, disturbances, and within and across-seasonal temperature and precipitation variation on spring and autumn LSP to clarify the role of climate dynamics. This approach addresses ways to overcome the monitoring challenges from microclimate, vegetational sensitivity and disturbance history that are common in many landscapes.

2. Materials and Methods

2.1. Study Area

We focused on the Great Smoky Mountains National Park (GRSM) in the Southern Appalachians, southeastern USA, as it supports particularly rich vegetational community types that occur across a highly complex topographic gradient (Figure 1). The 2114 km² (816 mi²) land unit is situated in the Blue Ridge Mountains of the Appalachian Mountain chain of North Carolina and Tennessee and ranges in elevation from 260 m (875 ft.) in stream valleys in the West to 2025 m (6643 ft.) on Clingman’s Dome. During the growing season, higher elevation sites are 10–15 °C cooler than are forests at lower elevation [47]. In part due to orographic precipitation, GRSM experiences some of the highest rainfall totals in eastern North America, though microclimatic differences are substantial [48]. Seasonal variation in temperature and precipitation is considerable from year to year and at longer time scales. Since 1980, lower elevation weather stations have warmed far more than have those at higher elevations [48].

GRSM was established in 1934 after well more than half of its area had been cleared by settlers, cutover by logging companies or burned by slash fires that radically altered the vegetation [49]. Substantial old growth forest remains, although most old growth has experienced mortality from non-native insects and diseases including chestnut blight and the balsam and hemlock woolly adelgids.

Topographic gradients drive most of the change in vegetation communities across GRSM and have inspired some of the earliest research on gradient analysis of vegetation [50]. Red spruce (Picea rubens)–Fraser fir (Abies fraseri) forests occur above 1400 m (4500 ft.), but the elevation of the ecotonal transition from the lower hardwood deciduous forest reflects slope aspect and site history [51,52]. During the 1960s, Fraser fir experienced high mortality from the balsam woolly adelgid, and while this species continues to recover, this non-native disturbance greatly impacted the spruce-fir forest. Eastern hemlock (Tsuga canadensis) is common on mesic sites, particularly along riparian zones of low and mid elevation, but eastern hemlock is often the dominant tree over wide areas on northerly aspects where it intergrades with the spruce fir ecotone [50]. Since the early 2000s, vast numbers of eastern hemlock have succumbed to the hemlock woolly adelgid where they were not biologically or chemically treated [53,54].

Much of the Park’s vegetation below the spruce-fir is deciduous hardwood forest, with dominant species varying by elevation, aspect, slope position and site history. Many of the Park’s deciduous forests have a complex deciduous understory comprised of saplings, shrubs, and ephemeral herbaceous plants. Evergreen understory shrub cover is also common, particularly in riparian areas and on mesic slopes where Rhododendron spp. is common and on xeric fire-prone slopes where mountain laurel (Kalmia latifolia) is common [55].

Pines provide a common evergreen cover type on xeric low and mid elevation slope. Fire-adapted table mountain pine (Pinus pungens), pitch pine (P. rigida) and shortleaf (P. echinata) prefer ridgelines and southern xeric slopes, often in association with oak (Quercus spp.) and hickory (Carya spp.) [56,57]. Virginia pine (P. virginiana) is locally common, especially in xeric, low elevation sites that are regenerating from past agricultural clearing. Numerous wildfires and low-intensity prescribed fires have occurred within the Park since 2001, with the two most extensive fires having occurred during extreme drought. The 2001 Sharp Fire burned over 24 km² (6000 acres), then in 2016, over 4000 ha (10,000 acres) of drought-stressed forests within GRSM burned by the Chimney Tops 2 Fire. Since 2000, prescribed fires of low severity have occurred on low elevation sites in the western portion of the Park on the Tennessee side. This same western section of GRSM experienced significant damage by a tornado in 2011 [58].

2.2. Land Surface Phenology

We used the 232 m ForWarn dataset of the US Forest Service and NASA that was derived from Moderate Resolution Imaging Spectroradiometer (MODIS) imagery with 8-day time steps or 46 values of Normalized Difference Vegetation Index (NDVI) per year [3,59]. For the 2000–2015 period of analysis, this includes 736 NDVI values for each of the 38,347 MODIS cells in the Park.

The processing steps used to develop the ForWarn data used in this analysis consisted of the following [33]: the MOD13 and MYD13 composites from the Terra and Aqua satellite data streams were decomposed to extract the specific day of observation used by those 16-day maximum value products. These two 16-day data streams are offset by 8 days. Extracted NDVI daily values were then filtered to remove atmospheric contamination, excess obliqueness, and aberrantly high or low outlying values. The two data streams were merged and combined into 8-day maximum value composites, followed by interpolation and curve fitting using a Savitzky-Golay interpolating technique [60]. Because of the 8-day maximum value compositing technique used to minimize clouds and atmospheric effects, there is an inherent week-long delay in autumn or whenever NDVI is trending downward, but not normally in spring when the highest NDVI during green-up is at the end of the 8-day period.

2.3. Comparison with Independent Data

Ground-based web-cameras that target natural vegetation (i.e., “PhenoCams”) provide an alternative daily record of seasonal changes that can be compared with satellite-based LSP time series [61,62,63]. In this study, the average midday green chromatic coordinate (gcc) was calculated by 8-day periods for two stations adjacent to the study area. The Smokylook PhenoCam (DB_0002) lies at 801 m elevation near the western edge of GRSM with data from 2003–2008, 2010–2015 (Latitude: 35.6325°, Longitude: −83.9431°). Second, lying 85 km (53 mi.) east of Smokylook and immediately east of GRSM, the Smokypurchase PhenoCam (DB_0003) tracks a forest stand at 1550 m (5000 ft.) elevation for the years 2008–2015 (Latitude: 35.5900°, Longitude: −83.0775°). See https://phenocam.sr.unh.edu/ for more information.

The trees tracked by these two PhenoCams are few in number and likely are unrepresentative of the species that provide the NDVI in a single 5.4 ha MODIS pixel. To see how well these PhenoCams represent the MODIS NDVI variation of the broader landscape, we extracted the mean 8-day MODIS NDVI for mid-spring (ending 16 May) and mid-autumn (ending 31 October) of those MODIS grid cells having purely deciduous forest, based on the 30 m-resolution 2006 National Land Cover Database (NLCD) (https://www.mrlc.gov/) [64]. Only MODIS cells that fell within a 100 m elevational band centered on the PhenoCam’s elevation across GRSM were used. We compared these mean landscape NDVI values with the 8-day mean gcc value for the respective prior period (i.e., for the periods ending 8 May and 23 October) to compensate for the maximum-value compositing used to process the MODIS time series. Comparisons were made using linear regression.

2.4. Pheno-Topographic Relationships

To describe the influence of topography on the timing of spring and autumn phenological transitions, elevation and slope parameters were calculated for every MODIS cell. Solar radiation, a surrogate for slope aspect and its associated dryness, was calculated for the summer equinox (21 June) with a 30 m digital elevation model using the Solar Radiation tool in ESRI’s ArcMap and assigned relative values from 1 to 100 [65]. Indexed solar radiation values were up-scaled to MODIS resolution using the mean neighborhood value in a 240 m moving window, then binned into three evenly represented classes based on the distribution of values within the Park (i.e., relatively xeric, moderate, or mesic). Elevation was binned into 152 m (500 ft.) classes from 457 to 1676 m (1500 to 5500 ft.). Because NDVI varies by vegetation type and topography, values were transformed into a Percent Completion metric, expressed as the percentage of the range between the average 2000–2015 summer maximum and winter minimum NDVI represented by the current NDVI for each MODIS cell. These Percent Completion NDVI values were summarized for each combination of elevation band and xeric/mesic class through spring and autumn to show the relative impacts of slope and elevation on phenology as these seasons progressed.

2.5. Pheno-Climatological Relationships

The strong temperature and precipitation gradients across GRSM mean that station data could misrepresent the local microclimate of areas captured by MODIS-resolution LSP. We used the 1-km² daily gridded Daymet v.3 dataset [66] to better reflect these gradients and to more precisely tailor the precipitation and temperature data to its corresponding MODIS cell. For each period, growing degree days (GDD) were accumulated using daily TMAX values that exceeded 10 °C. Chilling degree days (CDD) were accumulated using TMIN values below 5 °C. To capture the cumulative influence of GDD, CDD and precipitation (PCP), values were summed to 24-day periods. GDD and CDD may not capture the influence of short-term frost events, so we calculated the last spring frost and first autumn frost from gridded daily Daymet data, 2000–2015, using a threshold of 0 °C on TMIN.

Climate-NDVI relationships were explored quantitatively using correlation analysis that compared GDD, CDD and PCP with the mid-spring and mid-autumn NDVI across the 16 years of record. While most LSP researchers attempt to use a slope- or derivative-based algorithm to describe the dates for cell-specific yearly variations in spring and autumn seasonal transitions [67], we instead selected two fixed periods that were determined empirically from results of our pheno-topographic analysis. These a priori fixed dates correspond to mid-spring (ending 16 May) and mid-autumn (ending 31 October) for all cells in GRSM and are based on outcomes of the pheno-topographical analysis described above. This uniformity of assignment permitted a more straightforward cross-seasonal comparison of phenological differences across different elevations, aspects and vegetation composition mixtures.

For summary analysis and interpretation, MODIS grid cells were filtered according to cover type using the 30 m 2006 NLCD [64]. This allowed us to isolate the phenological responses of deciduous forests from that of conifers (hemlock, pine, spruce, fir) and mixed forest types independently. The response of conifers can differ from hardwoods, because the NDVI of conifer forests is less sensitive to change [2]. A strict pure hardwood forest filter was used to separate these two forest groups; if even a single 30 m cell of the 64 NLCD cells per MODIS cell was evergreen or mixed forest, the MODIS cell would still not be considered to be the deciduous forest type.

Forest disturbances and recovery that occurred between 2000 and 2015 may have altered the NDVI-microclimate relationship. To address this, records of recent disturbances within the Park were obtained from GRSM (https://irma.nps.gov). These data included areas of known hemlock occurrence to capture probable mortality from the hemlock woolly adelgid, a 2011 tornado and wildfires.

3. Results

3.1. Land Surface Phenology

Mean NDVI across the Park showed a strong, but lagged relationship with the pattern of DAYL, GDD and CDD during both spring and autumn (Figure 2). At the latitude of the Park (35.5°N), daylength peaks on the summer solstice, reaches its minimum on the winter solstice in late December, and is at 82% of the maximum on the vernal and autumnal equinoxes. The actual solar radiation received on sites depends on the way that daylength and sun angle relate to the particular slope, aspect, and topographic exposure at each location through the year.

DAYL steadily increases from January, preceding the increase in GDD by about a month in early spring, then DAYL begins to fall 3–4 weeks prior to GDD in late summer. In turn, NDVI lags GDD by an additional 3–5 weeks in early spring, and about four weeks in early autumn. In late summer, CDD begin to accumulate about three weeks after GDD begins to decline, with NDVI falling as CDD rises through the autumn. CDD continue to accumulate through the winter, peaking in early February, roughly midway between the winter solstice and spring equinox, with the majority of its decline occurring well before the visible start of the growing season across GRSM, which is mid-April. Minimum NDVI typically occurs in the weeks immediately after the peak in CDD from late February or early March. Peak NDVI occurs in July, soon after maximum DAYL, but well before GDD peaks in early August. Average PCP is highly irregular at eight-day resolution, but PCP dips in October as NDVI is falling.

Considering just grid cells of pure deciduous forest cover, both elevation and solar radiation (site moisture) influence NDVI in spring and autumn, but elevation is most important (Figure 3). Forest green-up advances earliest on mesic slopes of the Park’s lower elevations, then progresses to sunny slopes, then this sequence is repeated progressively upslope. In the autumn, this pattern reverses, with shaded sites dropping in NDVI before sunny slopes, with the greatest slope difference expressed at high elevations. In contrast to the strong elevational association that was observed for spring green-up, autumn decline progresses quickly, with relatively weak elevational control.

In spring, hardwood forests at elevations below 610 m (2000 ft.) achieve 15–20 Percent Completion of their maximum summer NDVI by 22 April, and over half by 8 May (Figure 3). In contrast, hardwood forests above 1372 m (4500 ft.) achieve 22–28% of their total green-up by 8 May, with their half-way mark about two weeks after that. This elevational lag in green-up changes over the course of spring, with the highest and lowest elevation bands differing in state of green-up by just one period (eight days) in early spring, but by twice that by late spring. Autumn decline commences at a mean NDVI value that lies about 95 Percent Completion of the maximum growing season NDVI (which had occurred in early July (see Figure 2), drops slowly at first, and then reaches a maximum rate of decline in late October.

Across years, variance in the timing and progression of both spring and autumn was substantial, with each season varying by about 20 days (Figure 4). Most springs that started early continued to progress as early, but 2005 and 2007 started early, then ended late. The cross-year distribution of the timing of springs is more symmetrical than for autumns, with 2002, 2005 and 2007 being particularly late autumn outliers, while in spring, only the early spring of 2012 appears as an outlier.

3.2. Comparison with Independent Data

Comparison of MODIS NDVI with PhenoCam data across years shows results that are directionally consistent in magnitude, with spring more highly correlated than autumn (Figure 5). In spring, the high-elevation Smokpurchase PhenoCam data were more strongly correlated with NDVI (p = 0.003; R² = 0.914, N = 6) than at the lower elevation Smokylook site (p < 0.001; R² = 0.777, N = 12). In autumn, MODIS NDVI comparisons with Smokylook were also weaker (p = 0.012; R² = 0.418, N = 14) than with Smokypurchase (p = 0.051; R² = 0.496, N = 8).

3.3. Pheno-Climatology

Dates for the mid green-up and mid brown-down were empirically determined from the multi-year LSP curves shown in Figure 3 and Figure 4; the NDVI periods ending 16 May and 31 October were chosen for comparison with three temperature and precipitation factors, as outlined below. The distributions of correlations that relate spring and autumn NDVI with GDD, CDD and PCP differ for deciduous and conifer-mixed vegetation types (Figure 6). Deciduous forests usually exhibit stronger positive or negative correlations than do their paired conifer-mixed forest counterparts.

Correlations between mid-spring NDVI (for the period ending 16 May) and lagged GDD are positive for the lagged weather periods ending 29 March, 14 April and 30 April, meaning that greener mid-spring conditions are associated with warmer temperatures in the weeks prior to mid-spring (Figure 6A). This correlation becomes inverse for lagged GDD that accumulate for the periods ending 8 and 24 January, suggesting a possible association between years with earlier spring green-up and January coldness. Correlations with CDD are generally the reverse of those with GDD, with a colder early spring being associated with a delay in green-up (Figure 6B). Results also suggest a possible longer-reaching cross-seasonal relationship between early springs and coldness during the prior autumn (i.e., for the 24-day period ending 7 October) (Figure 6A,B). Correlations with spring NDVI and lagged PCP are weak, but positive for 24 November and 9 February (Figure 6C).

Correlations between mid-autumn NDVI (for the period ending 31 October) and lagged GDD are strongly positive for 29 September and 15 October, with a predictable inverse pattern with CDD (Figure 6D,E). A cross-seasonal association is suggested between late autumns (i.e., higher autumn NDVI) and cool springs (i.e., more CDD) for the period ending 8 May. As seen with lagged PCP in the spring, fall PCP shows a weak relationship with autumn NDVI (Figure 6F).

When selected lagged correlations are mapped, the statistical distribution of Figure 6 exhibits coherent spatial structure in terms of correlation strength and sign. Correlations of mid-spring NDVI (16 May) and GDD accumulated from 7 April to 30 April show a stronger positive relationship for lower elevations, particularly to the south and west (Figure 7A). Similar spatial patterns occur with spring NDVI and lagged spring CDD, but correlations are reversed (Figure 7B), consistent with Figure 6A,B. Correlations of spring PCP with mid-spring NDVI show strong geographic differences. The northern and southern foothills are strongly positive, while low to mid-elevations of the west and east, and higher elevations more generally, exhibit weak to negative correlations (Figure 7C).

Mapped autumn correlations with mid-autumn NDVI (for the period ending 31 October) and GDD (accumulated from 6 to 29 September) are particularly strong and positive over nearly all of the Park (Figure 7D). Exceptions are the grass-dominated Cades Cove in the west, higher elevations to the east, and similar lower elevations in the east that were anomalous in spring. Autumn NDVI and autumn CDD largely mirror the pattern observed with GDD, though reversed (Figure 7E). Autumn NDVI and PCP correlations are weakly, but broadly negative across most lower and mid elevations, implying that drier conditions relate to delayed senescence. PCP correlations are consistently uncorrelated or only weakly correlated in many of the same locations as observed by other correlations (Figure 7F).

The combined influence of temperature and precipitation on spring and autumn NDVI are conveyed by Figure 8. Cool-dry springs are more likely to have a lower NDVI in May (represented by smaller diameter circles in Figure 8A) than are warm-wet springs, although moderately warm, but dry springs are occasionally abnormally green (2004 and 2010, at center-left of Figure 8A). Extreme warm-wet or cool-dry conditions are not represented within this 2000–2015 measurement period in autumn, but late senescence (indicated by larger diameter circles in Figure 8B) is more likely when conditions are warm and dry rather than cool and wet.

In spring, the last frost date varied by 21 days with a Park-wide average date of 11 April (

Table 1. Average gridded first and last frost dates (<0 °C) for the Park, 2000–2015. Rank indicates position out of 16 years arranged in decreasing sorted order. Season is length in days.

Year	Spring	Rank	Autumn	Rank	Season	Rank
2000	16 April	13	13 October	1	198	6
2001	19 April	16	13 October	1	182	14
2002	16 April	12	3 November	13	191	8
2003	1 April	2	3 November	13	202	3
2004	15 April	11	7 November	16	190	10
2005	14 April	9	30 October	10	206	1
2006	5 April	3	17 October	3	189	11
2007	12 April	8	1 November	11	185	13
2008	18 April	15	26 October	8	199	5
2009	11 April	6	23 October	5	191	8
2010	29 March	1	1 November	11	192	7
2011	14 April	9	24 October	6	200	4
2012	10 April	5	4 November	15	203	2
2013	11 April	6	24 October	6	188	12
2014	16 April	12	22 October	4	174	15
2015	5 April	3	27 October	9	171	16
Mean	11 April	-	26 October	-	191	-

). The earliest frost date in autumn varied by 25 days, with a Park-wide average of 26 October. Growing season length, measured as the number of days between the last and first frost, ranged from 171 days (in 2015) to 206 days (in 2005).

Places with anomalous correlations in Figure 7 usually occur where disturbance has occurred over the 2000–2015 era (Figure 9). A 2011 tornado scar shows up in the west in all three spring correlation maps (Figure 7A–C), but not in the autumn maps (Figure 7D–F). Wildfires show no discernable impact on correlations, but most of the higher elevation and eastern forests with low or anomalous correlations are where hemlock is dominant or co-dominant and where the non-native hemlock woolly adelgid has led to widespread hemlock mortality. Mortality is widespread because treatment, where it occurs, targets only a fraction of the hemlock trees that are present.

4. Discussion

The temporal and spatial resolution of MODIS LSP successfully conveys the influence of environmental gradients on spring and autumn land surface phenology across GRSM. Temperature and precipitation’s fundamental influence on LSP is suggested by Figure 2, whereby the late winter and spring increase in day length and solar radiation drive changes in temperature, which, in turn, leads to an abrupt increase in NDVI in early April. NDVI persists through summer after GDD and DAYL decline, then the onset of autumn is associated with increased CDD in late September.

Topography provides a secondary influence on LSP across GRSM, with elevation, in particular, controlling the onset and progression of spring and autumn (Figure 3). The impacts of aspect are secondary, but it is difficult to know if aspect differences relate to microclimate or to differences in species composition that also conform to aspect. The importance of elevation in this research is consistent with results from prior research in mountainous terrain [22,68,69,70,71,72,73], but there is a clear need for more research regarding how topographic factors beyond elevation affect spring and autumn timing.

Our documented spring-to-autumn reversal in the order of seasonal change across slopes is notable, particularly because it represents the long-term average response over a large landscape (Figure 3). In spring, mesic slopes are first to advance, while in autumn they are the first to senesce. This pattern could be explained by microclimate or by the attributes of the dominant species that occur across this moisture gradient. While dry south-facing slopes warm relatively early in the spring, and may experience relatively early water stress in autumn, our observed LSP patterns do not show these expected LSP responses—sunny slopes advance after shaded slopes in spring, and water stress would be expected to be observed first on drier slopes. Instead, these patterns may relate to species’ phenological attributes. Xerophytic species often leaf out weeks after mesic species in spring [27]. In autumn, a combination of differential species responses and microclimates having greater exposure to cooler autumn temperatures may explain these long-term average patterns.

Prior studies have measured spring and autumn using NDVI slopes or inflection points to estimate local dates of transitions, while in this research a less complex approach is used that relies on empirically selecting fixed dates for spring and autumn, in order to directly compare seasonality across climatic and topographic gradients. In Figure 3, comparability across sites within elevation bands is gained by calculating values as the Percent of Completion. This conversion reveals the extent of the summer reduction in NDVI from leaf aging, moisture stress or disturbance that occurs prior to the autumn season—a phenomenon that has been described in the northeastern US, but is not widely recognized [74]. In GRSM, neither low summer PCP nor warm summer temperatures (GDD) are strongly associated with delayed senescence (Figure 6 and Figure 7), although they may increase drought stress in concert (Figure 8).

Lower elevation forests may be more vulnerable to drought-induced late-growing season stress than are forests at mid- and higher elevations, with implications for future forest productivity [22]. In general, autumn warmth delays senescence, while drought advances the onset of autumn [6,37]. While late-season hydrologic stress and productivity may not be consistently measured by autumn NDVI, here we find no strong support for summer or autumn drought-induced senescence (Figure 6F and Figure 8B). While neither extreme warm-wet nor cool-dry autumns occurred during our 2000–2015 analysis period, warm-dry autumns were generally associated with delayed senescence and cool-wet autumns showed earlier brown-down, as expected (Figure 8B). The earliest autumns—2011 and 2012—appear as outliers to this trend, and they may have been caused by factors other than GDD or PCP as measured, but they are not explained by unusually early autumn frost dates (Table 1).

An early spring brings an increased hazard of frost damage, and growing season duration affects community productivity and competitive relationships. The extremely early spring of 2012 also experienced an early senescence. Compare the 2012 season with 2013, where both spring and autumn were late, and also with 2002, which experienced a particularly long growing season (Figure 4). LSP-based phenology results are often inconsistent with growing season expectations from the last and first frost dates alone (Table 1), and this highlights a need for better understanding the mechanisms that regulate the growing season in agricultural and wildland systems. The spring of 2000, 2001, 2007 and 2014 started early, but lost momentum (Figure 4A), and this helps explain their relative position in Figure 8A, since early green-up was caused by the early accumulation of GDD that then slowed. Each of these years had either a late-ranking frost (Table 1) or, in the case of 2007, a severe April freeze [20].

Our results are directionally consistent with a growing body of evidence regarding within-season and cross-seasonal temperature effects on spring and autumn LSP [24,36,45]. Recognition of within-season lags has long been acknowledged and captured through use of cumulative GDD and cumulative CDD to predict spring and autumn phenology [75]. Our strongest relationships are for within-season GDD in autumn, but we also show that on some sites, spring may be accelerated by cold winters and autumn may be delayed by cold springs. These cross-seasonal correlations do not, however, necessarily imply some causative physiological mechanism, since chance combinations of observed seasonal weather over the short 16-year period or some inherent cross-seasonal climate mechanisms could generate these results, with LSP only being responsive to shorter-term climatic factors. If, however, there are causal biological mechanisms at work, such cross-seasonal complexities may make it more difficult to accurately predict how LSP will respond to climate change.

Improved precision in monitoring provides an important path forward. By filtering our analysis according to cover type, we demonstrate that deciduous and conifer/mixed forest types respond differently to antecedent seasonal weather conditions (Figure 6). This filtering was particularly important in GRSM because dying evergreen hemlock trees and a tornado during the period of analysis likely reduced the statistical strength of correlations between LSP and climatic and topographic variables, although fire had no clear impact, presumably because severity was relatively low, with most mortality confined to the understory (Figure 7 and Figure 9). In disturbed areas, the timing of green-up may increasingly reflect the understory signal, which can differ from that of the overstory [31,32,76]. With improved understanding of where disturbances have occurred or how vegetational responses have changed because of disturbance, phenological observations can be refined to improve landscape monitoring.

5. Conclusions

In this study, we report on the behavior of MODIS LSP along environmental gradients among years having different combinations of seasonal weather. We document stronger effects from elevation than from slope or vegetation type. We found that deciduous hardwood forests exhibit a much stronger relationship with climatic variation than conifer and mixed-forest types. We also found a particularly weak 16-year temperature and precipitation relationship with LSP in areas of recent acute disturbance, including areas affected by the hemlock wooly adelgid and a recent tornado.

Across years, we found that spring NDVI was strongly correlated with accumulated warming, with up to a six-week lag at lower elevations, and that autumn NDVI was extended when autumn is warm more broadly. We also describe possible cross-seasonal relationships. Spring green-up was sometimes correlated with colder winters. Autumn greenness was correlated with colder prior-springs. The influence of precipitation on seasonal NDVI was weak. Wetter springs, but drier autumns, may extend the growing season, but this was dependent on the dominant influence of temperature.

Remote sensing provides an efficient all-landscape means for tracking growing season dynamics of vegetation. However, just as vegetation typologies do not perfectly convey the species assemblages within a given vegetation type, so the phenological responses of individual species may not be discernable from LSP gridded data, particularly in species-rich forests such as those of the Southern Appalachians. With numerous species and individuals contributing to the LSP behavior of a grid cell, observed LSP responses may not equate to ground-based phenological observations of particular species any more than characteristics of individuals relate to communities.

This research describes remarkable variation in spring and fall vegetational phenology among years and across environmental gradients. For the plant and animal species of these mountains, these spatial and temporal attributes constitute a rich diversity of phenological habitat that enriches the evolutionary environment of this place. Researchers have long-argued that the topographic and microclimatic complexity of the Great Smoky Mountains, and of the Southern Appalachian Mountains more broadly, have contributed to the region’s rich biodiversity [50,77,78]. This research shows that, while vegetation phenology responds to those fixed environmental attributes, it also captures a dynamic year-to-year and cross-seasonal flux that is superimposed on those stable gradients. This flux includes the variable vegetative responses to spring and autumn weather across years, and the impacts of localized biotic and abiotic disturbances that alter phenological conditions at longer time scales. While such phenological complexity can be difficult to monitor, filtering observations with disturbance history and vegetation type can clarify how landscape vegetation responds to seasonal forcing.