Spatio-Temporal Variability in Remotely Sensed Vegetation Greenness Across Yellowstone National Park

The study’s objective was to quantify the responses of vegetation greenness and productivity to climate variability and change across complex topographic, climatic, and ecological gradients in Yellowstone National Park through the use of remotely sensed data. The climate change signal in Yellowstone was pronounced, including substantial warming, an abrupt decline in snowpack, and more frequent droughts. While phenological studies are increasing in Yellowstone, the near absence of long-term and continuous ground-based phenological measurements motivated the study’s application of remotely sensed data to aid in identifying ecological vulnerabilities and guide resource management in light of on ongoing environmental change. Correlation, time-series, and empirical orthogonal function analyses for 1982–2015 focused on Daymet data and vegetation indices (VIs) from the Advanced Very High-Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS). The study’s key questions address unique time scales. First, what are the dominant meteorological drivers of variability in vegetation greenness on seasonal to interannual time scales? Key results include: (1) Green-up is the most elevationand climate-sensitive phenological stage, with La Niña-induced cool, wet conditions or an anomalously deep snowpack delaying the green-up wave. (2) Drought measures were the dominant contributors towards phenological variability, as winter–spring drought corresponded to enhanced April–June greening and spring–summer drought corresponded to reduced August–September greening. Second, how have patterns of productivity changed in response to climate change and disturbances? Key results include: (1) The park predominantly exhibited positive productivity trends, associated with lodgepole pine re-establishment and growth following the 1988 fires. (2) Landscapes which were undisturbed by the 1988 fires showed no apparent sign of warming-induced greening. This study motivates a systematic investigation of remote-sensing data across western parks to identify ecological vulnerabilities and support the development of climate change vulnerability assessments and adaptation strategies.


Introduction
The Western United States, a region renowned for water scarcity, is a hotspot for recent and future projected climate change.Recent trends of warming, amplified droughts, diminished snowpack, and expanded fire season [1,2] have generated concern among land managers and private citizens.The projected warming and drought-related water shortages for the Western United States are expected to alter terrestrial ecosystems [3][4][5] and induce pronounced shifts in species ranges and biodiversity patterns [6], although with such studies placing a greater focus on the water-limited Southern Rockies than the temperature-limited Middle Rockies [7,8].Climate change and associated global-change-type droughts [3] pose an unprecedented risk to all lands of the Western United States, including the many resource-based national parks therein [9][10][11].The National Park Service (NPS) was created in 1916 and is charged with the mandate to protect these parks "in such a manner and by such means as will leave them unimpaired for future generations" [10].However, recent changes in climate and disturbance regimes are affecting ecosystem functions, such as net primary productivity, phenology, and post-disturbance succession.Phenological variability and trends affect park resource management through impacts on predator-prey, plant-herbivore, and plant-pollinator interactions [12][13][14]; species' phenotypic mismatch and resulting food web disruption and biodiversity loss [15][16][17]; large-scale movements of herbivores tracking green-up-associated resource waves [18]; spread of invasive species with flexible phenologies and need for mechanical removal or herbicide application [19,20]; and altered timing of cultural festivals linked to flower displays and migrations [14,21].According to the green wave hypothesis, plant phenology critically shapes the resource landscape and migratory ungulates are found to track the green wave, with implications of climate change affecting the persistence of migratory taxa [18].In order to effectively steward national parks into this century of climate change, land managers need the best possible information about how the ecosystems they manage respond to climate variability, particularly changes in temperature, snowpack dynamics, and drought.Here, we quantify ecosystem-level responses to climate variability in an effort to support land manager decision-making within an important NPS asset, establishing a framework for assessing climate-ecosystem interactions that is both location-specific and broadly extendible.
The present study focuses on Yellowstone National Park, which is the oldest official national park in the United States and the second largest park in area within the lower 48 states.Its largely pristine environment; pronounced topographic, climatic, and ecological gradients; and well-documented fire and bark beetle (Dendroctonus spp.) disturbances make Yellowstone National Park an ideal focal region to study ecological responses to climate variability and change.The park's heterogeneous elevational gradients lead to complex patterns of temperature and precipitation, which in turn control vegetation growth and phenological responses.In order to quantify ecosystem responses to climate and disturbance, we employed a variety of remotely sensed and gridded data products.Remotely sensed data products provide many advantages, including wall-to-wall coverage of Yellowstone National Park, a temporal range that encompasses the 1988 fires, and a high sensitivity to ecologically relevant reflectance properties.These advantages are particularly important for NPS units in the Western United States, where on-the-ground data can be sparse and temporally discontinuous.Of primary interest, vegetation indices (VIs) serve as radiometric measures of the amount of photosynthetically active radiation absorbed by chlorophyll within canopies' green leaves, and thus, valuable surrogate measures of the level of physiologically functioning surface greenness across a region [22][23][24].Vegetation indices reveal changes in biomass, biodiversity, carbon sequestration, and plant water stress [25][26][27], all key metrics for land resource managers.
There have been a limited number of remote sensing-based vegetation monitoring studies specifically focused on Western United States' national parks, including Bandelier [28], Big Bend [29], Crater Lake [30], Glacier [31], Saguaro [32], Yellowstone [33][34][35][36][37][38], and Yosemite [39,40].The Yellowstone studies of Franks et al. [33], Jakubauskas and Price [34], Zhao et al. [35], Emmett et al. [36], Potter [37], and Garroutte et al. [38] focused on relating remotely sensed VIs to observed leaf area index (LAI) and net primary productivity, forest structure, post-disturbance fire recovery, vegetation greening trends, analysis of vegetation cover changes, and relating remotely sensed VIs to grassland biomass and quality, respectively.Based on regression analysis, Franks et al. [33] determined that Landsat-derived Normalized Difference Vegetation Index (NDVI) and band-5 reflectance metrics could describe about 60-70% of the variability found in ground-based measurements of LAI and annual net primary productivity over a six-year fire recovery period; this led to the conclusion that spatial variability in regrowth rates can be inferred from temporal trends in satellite vegetation metrics.Jakubauskas and Price [34] demonstrated the value of remote sensing for examining the biotic factors of gross physical structure for Yellowstone's lodgepole pine forest canopies by regressing forest overstory and understory field data against Landsat radiance values and transformed data to quantify the reliability of the remote sensing.By assessing long-term forest spectral recovery following disturbances in Landsat data via the application of the Landsat Time Series Stack-Vegetation Change Tracker algorithm, Zhao et al. [35] found that forest recovery is sensitive to forest type, soil type, and elevation.Applying linear least squares regression and the Mann-Kendall significance test, Emmett et al. [36] focused on the greening trend across the Greater Yellowstone Ecosystem and the relative importance of disturbance and climate in explaining trends in United States Geological Survey (USGS) Earth Resources Observation and Science (EROS) NDVI during 1989-2014.Consistent with Potter [37], the study concluded that disturbance history was a major driver of the park's NDVI trends and that it is essential to analyze spatial patterns of trends in plant productivity using very high-resolution data.The aforementioned studies proved that VIs can be effective measures of ecosystem attributes, such as LAI [33] and gross physical structure [34].While these individual studies contribute to our understanding of discrete phenomena, NPS leadership has identified the need for synthesizing the complex interplay of multiple climate and vegetation signals in support of management decision-making [41].
Of the many previously cited papers which use space-borne remote sensing to investigate ecological questions within individual western national parks, it is surprising that only two focus on vegetation phenology specifically: O'Leary et al. [30], who stressed the impact of snowmelt timing on green-up, and Wallace et al. [32], who mapped invasive buffelgrass and its phenology.This highlights a notable gap in the literature, considering that the importance of long-term, continuous remote sensing data was emphasized by White and Swint [29] and Soulard et al. [40], while most of the remote sensing studies of the western national parks have analyzed only a limited number of snapshot images [28,31,39] or 3-5 years of satellite data [29,32,37].As pointed out by O'Leary et al. [30], conifer regions have especially received insufficient attention in remote sensing phenological studies.Furthermore, such studies have often considered minimal or no climate datasets in their analysis of VI patterns, thereby neglecting the importance of spatio-temporal variations in meteorological conditions in driving phenological responses; these studies used topographic variations in elevation, slope, and aspect as proxies for temperature, precipitation, and solar radiation [28,30,31,40] and often concluded that elevation is the main driver of NDVI's spatial heterogeneity [31].Therefore, there is a need to quantify the complex interactions between the many important climate variables and wildfire disturbance as drivers of vegetation phenology and productivity, particularly within an important natural resource where changes in climate and disturbance patterns are expected into the future.This study aims to address this need by quantifying the spatio-temporal variability in the vegetation greenness of Yellowstone National Park in response to changing climatic and disturbance conditions across multiple spatio-temporal scales and ecosystems.
Here, we identify the impacts of climate variability/change and wildfire on vegetation greenness and productivity at seasonal, annual, and decadal timescales.Analyses presented here address the following questions: (1) What are the dominant meteorological drivers of variability in vegetation greenness across Yellowstone National Park on the seasonal to interannual time scales?(2) How have patterns of productivity changed across the park in recent decades in response to climate change and disturbances?The study's findings offer new insight into Yellowstone's vulnerability to climate variability and change relevant to park leadership, wildfire management, ecology, and a highly invested public.

Study Area
Yellowstone National Park is characterized by pronounced topographic, climatic, and ecological gradients across its vast 8983 km 2 area, located primarily in Wyoming and extending into Idaho and Montana [42][43][44].The park is comprised of five geologic provinces according to bedrock type, namely the Absaroka Range, Gallatin Range, Southwest Plateaus, Central Plateaus, and Yellowstone-Lamar Valleys [42].The park's mean elevation is 2453 m, ranging dramatically from 1610 m at Reese Creek near Gardiner, Montana to 3466 m at Eagle Peak within the Absaroka Range, which is a sub-range of the Rocky Mountains along the park's eastern boundary (Figure 1).The lower-elevation Lamar Valley spans the park's north and contains the Lamar River, which is a tributary of the Yellowstone River, and is known for abundant wildlife.Yellowstone Lake is the largest North American high-elevation (above 2.1 km) lake by area.Yellowstone National Park is largely dominated by sub-alpine fir landscapes, with the most abundant habitats consisting of lodgepole pine (Pinus contorta, 64.0% of area) due to fire history and soil conditions [45], followed by non-forested landscapes (13.4%), whitebark pine (Pinus albicaulis, 12.8%), Douglas fir (Pseudotsuga menziesii, 6.0%), Engelmann spruce and sub-alpine fir (Picea engelmannii and Abies lasiocarpa, 3.2%), aspen (Populus tremuloides, 0.2%), and krummholz (various species of stunted trees near tree line, 0.1%) (Figure 2) [37,42,[45][46][47].The distribution of vegetation types in the park generally follows an elevational gradient, with aspen at the lowest elevations and krummholz at the highest elevations (Figure 2).Krummholz and aspen are typically confined to a narrow range of suitable elevations, while non-forested landscapes are distributed across a diverse range of elevations.The climate is characterized by short, cool summers and long, cold winters [42]. 5

Applied Topographic, Cover Type, Soil, Disturbance, and Climatic Datasets
This study applied commonly used datasets that were individually calibrated and validated, with details about quality control found in their respective cited literature.Elevation was retrieved, and slope and aspect were computed based on each grid cell and its eight neighbors (using the arctangent function for aspect in the National Center for Atmospheric Research (NCAR) Command Language (NCL, www.ncl.ucar.edu)),from the 800-m digital elevation model dataset from Oregon Yellowstone National Park has experienced distinct environmental changes in the last century, including rising air temperatures [48][49][50], earlier spring onset [51,52], declining snowpack and earlier snowmelt [37,48,53], expanding growing season and fire season [48,50], shrinking wetlands [54,55], declining amphibian species richness [56], an eruption of mountain pine beetle outbreaks and whitebark pine mortality (Pinus albicaulis) [49,[57][58][59][60][61], upslope shift of lodgepole pine (Pinus contorta) [10,52,62], and upslope expansion of seedling aspen (Populus tremuloides) after the 1988 fires [63,64].The 1988 Yellowstone fires, which were believed to be facilitated by La Niña dry conditions [65], summer drought, lightning, and strong winds, burned approximately 4850 km 2 of the Greater Yellowstone Ecosystem, primarily conifer-dominated forests [32,[66][67][68][69]. Approximately 250 different fires burned across Yellowstone and nearby national forests during June-August 1988, scorching about 36% of the park with the greatest impact across its western North Fork region [66][67][68][69][70][71].These fires provide a natural experiment, allowing for the paired comparison of disturbed and undisturbed sites across multiple vegetation classes, as well as the quantification of ecosystem succession, illustrating a variety of ecosystem responses to climate and disturbance which are all documented in the available satellite record.This study applied commonly used datasets that were individually calibrated and validated, with details about quality control found in their respective cited literature.Elevation was retrieved, and slope and aspect were computed based on each grid cell and its eight neighbors (using the arctangent function for aspect in the National Center for Atmospheric Research (NCAR) Command Language (NCL, www.ncl.ucar.edu)),from the 800-m digital elevation model dataset from Oregon State University, as used by the Parameter-elevation Regressions on Independent Slopes Model (PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu)[72].The primary plant type distribution was determined from the NPS' YELL 1999 Cover Type dataset, as provided by the NPS, which consists of a geodatabase of land cover and habitat types layers based on 1:15,840 scale color aerial photography and field surveys [42].The cover type dataset was processed into eight categories: aspen, Douglas fir (combining three classifications, namely climax, post-disturbance, and successional Douglas fir), Engelmann spruce and sub-alpine fir, krummholz, lodgepole pine (combining four classifications, namely climax, post-disturbance, successional, and pygmy lodgepole pine), non-forested, water, and whitebark pine (combining three classifications, namely climax, post-disturbance, and successional).A fire extent perimeter polygon dataset, containing shapefile perimeters of fires larger than 0.4 km 2 , for 1881-2016 was acquired from the Yellowstone National Park Spatial Analysis Center (provided by Alex Zaideman).
Monthly Palmer Drought Severity Indices (PDSI) [73,74] and Standardized Precipitation Evapotranspiration Indices (SPEI) [75] were obtained, averaged across Wyoming's Park County, through the West Wide Drought Tracker website (https://wrcc.dri.edu/wwdt) of the Western Regional Climate Center [76].The SPEI combines the advantages of the Standardized Precipitation Index and PDSI, namely, the consideration of different timescales and the simultaneous consideration of both temperature and precipitation, respectively, leading to a potentially more meaningful measure of drought impacts on vegetation [77,78].Specifically, SPEI was considered here at time scales of one (SPEI1mn) and six months (SPEI6mn) to consider both short-and long-term drought and pluvial impacts.Analyses applied air temperature, precipitation, incoming surface solar radiation, and snow-water equivalent (snow depth × density) on a 1 km × 1 km grid from the Daymet Daily Surface Weather and Climatological Summaries [79,80], which are based on the Global Historical Climatology Network (GHCN)-Daily dataset distributed by the National Centers for Environmental Information (NCEI) [81].The Daymet V3.0 data were acquired through the Oak Ridge National Laboratory website (https://daymet.ornl.gov).

Remotely Sensed Vegetation Indices
In an effort to assess the robustness of the results, this study considers two remotely sensed VIs, namely the NDVI and the Enhanced Vegetation Index (EVI), the former of which is the most commonly applied VI (with an extensive time duration, not requiring blue-band data as in the EVI which was unavailable on earlier platforms) and the latter of which aims to reduce known deficiencies in the NDVI related to atmospheric and background effects.The NDVI and EVI are commonly used as proxies [24] for phenology [82,83], aboveground biomass [84,85], primary productivity [86][87][88], forage quality [89,90], and chlorophyll content [91,92].The NDVI is computed based on the reflectances in the near-infrared and visible channels, ranges from −1 to 0 for bare ground to 1 for dense vegetation and is sensitive to chlorophyll concentration [93].The EVI is calculated based on the near-infrared, red, and blue bands; is sensitive to canopy structural variations; and is intended to represent an improvement over NDVI by minimizing canopy background variations, enhancing sensitivity over densely vegetated canopies, and removing residual atmospheric contamination associated with thin cirrus clouds and smoke [94][95][96][97].However, EVI is more vulnerable than NDVI to contamination associated with heterogeneous terrain [24,95,97,98].The phenological development of the plant canopy is closely linked to NDVI seasonality [99].As the plant canopy shifts from early spring growth to late-season maturity and senescence, leaf reflectances change, as captured through remotely sensed NDVI.
The Advanced Very High-Resolution Radiometer (AVHRR) 3rd generation Global Inventory Monitoring and Modeling System (GIMMS) NDVI (NDVI3g) dataset from Boston University, with 1/12th degree spatial resolution, covers the period of 1982-2015 [100].The NDVI3g dataset, retrieved from the National Aeronautics and Space Administration (NASA) Monitoring, Modeling, and Forecasting Ecosystem Change website (https://ecocast.arc.nasa.gov)was generated by Zhu et al. [101] based on multiple AVHRR sensors using a feed-forward neural network algorithm, while accounting for calibration loss, orbital drift, and volcanic eruptions.The Moderate Resolution Imaging Spectroradiometer (MODIS) MOD13Q1 Vegetation Indices level-3 dataset from the NASA Earth Observation System, with a 250-m spatial resolution, covers the period of 2000-2017 [82,94,102]; this dataset was obtained through the Reverb website of the NASA Earth Observing System Data and Information System (https://earthdata.nasa.gov)for horizontal tile 10 and vertical tile 4. Values of MODIS NDVI range from −2000 to 10,000, with a scale conversion factor of 10,000.For each 16-day period, the MOD13Q1 algorithm selects the most reliable pixel value from all data acquisitions based on low cloud abundance, aerosol loading, view angle, and peak NDVI value.Values over water or outside of the shapefile boundaries of the park were masked.By averaging values that fall within a given month, these 16-day VIs were converted to monthly values for the current study (as was the daily climate data) to focus on sub-seasonal to seasonal time-scales.Based on comparison of MODIS VIs with field observations of canopy reflectance, Miura et al. [103] concluded that mean uncertainties are ±0.01VI units for NDVI and ±0.02 VI units for EVI under typical atmospheric conditions with 2% reflectance calibration uncertainty.The MODIS data reliability and quality are described in the Supplementary Materials.

Summary of Methodology
All meteorological and topographic datasets were bilinearly re-gridded to the 250-m MODIS VI grid using the ESMF_regrid function in NCL, which applies the Earth System Modeling Framework software.According to bilinear interpolation, each destination point was mapped to a location within the source grid, while the position of the destination point relative to the source points was used to compute the interpolation weights.Further testing of the study's results, instead using nearest neighbor re-gridding (not shown), reveals nearly identical findings.The spatial range of the study region was 44.13 • N-45.11 • N, 111.15 • W-109.83• W. The VI imagery was examined for 1982-2015 for AVHRR data and 2000-2017 for MODIS data.Applied statistical methods, as expanded upon below, include the Theil-Sen estimate of linear trend, Mann-Kendall test of trend significance, Pearson's correlation, Spearman's rho test, and unrotated empirical orthogonal function (EOF) analysis, using NCL's trend_manken, rgrid2rcm, rtest, escorc, spcorr, and eofunc functions.The applied significance level was 0.05.
Here, more specific details are provided on the applied methodology used to address the two primary scientific questions (Figure 3).The first question focuses on the dominant meteorological drivers of variability in vegetation greenness on the seasonal to interannual timescales.The pronounced spatial heterogeneity of Yellowstone's climate is demonstrated by compositing climatological maps of Daymet summer daily maximum temperature, winter precipitation, annual solar radiation, and annual maximum snow water equivalent.Likewise, the distinct spatial heterogeneity of vegetation greenness is shown through composites of annual mean MODIS NDVI and EVI, along with mean peak-month NDVI and EVI.When computing spatial correlations among topographic, ecological, and climatic datasets, one must address the potential for large spatial autocorrelation that can artificially inflate the correlation coefficient.As a consequence, semi-variogram graphs were generated using the "gstat" R package [104] to quantify semi-variance, which measures the spatial dependence between observations as a function of the distance between them.Specifically, plots of distance versus semi-variance were generated for eight select comparisons: annual mean air temperature versus elevation, annual mean precipitation versus elevation, annual mean EVI versus annual mean NDI, annual mean NDVI versus annual mean air temperature, annual mean NDVI versus annual mean precipitation, annual mean EVI versus annual mean air temperature, annual mean EVI versus annual mean precipitation, and the annual seasonal maximum of EVI versus the annual seasonal maximum of NDVI.The semi-variograms from the residuals of these fitted linear models provide the effective radius of remaining spatial autocorrelation; for the eight aforementioned comparisons, the maximum radius was 750 m.Data were randomly sampled at a greater distance than that radius to remove the effects of spatial autocorrelation, leading to an adjusted sample size (N) of 17,481; this newly determined N was applied when assessing statistical significant for all spatial correlations.
The seasonality and dominant modes of variability in vegetation phenological growth, and their sensitivity to climate and topography, were examined through EOF analysis of MODIS detrended NDVI data during the non-winter months of April-November.Also known as principal component analysis, EOF analysis permits a simplified interpretation of complex data across the space-time domain (ideal for the vast, high-resolution Yellowstone datasets) by decomposing a spatio-temporal field into dominant orthogonal spatial patterns and their associated uncorrelated time-series, or principal components, while reducing noise in the data [105][106][107][108][109]. Typically, analysis focuses on the leading EOF modes, as they contain most of the data's original variance [109].Through further EOF analysis, the leading three principal components were examined for monthly MODIS NDVI anomalies (after removing the mean seasonal cycle), both un-detrended and linearly detrended, in order to understand the dominant modes of phenological variability on the seasonal to interannual time scale.By averaging the loading pattern (magnitude of the spatial pattern of variability) of the dominant EOF modes across 200-m elevational bins, the preferred elevation for specific patterns of variability was determined per calendar month.The amount of variance explained by the principal components was a metric for evaluating the efficacy of EOF analysis.For the EOF analysis, the study did not attempt to mask snow-covered regions and periods, which were characterized by low NDVI values.The primary climatic regulators of variability in vegetation greenness were identified, by month, through Pearson correlations between monthly detrended AVHRR NDVI anomalies and antecedent atmospheric conditions, leading by 0 to 7 months, from Daymet and Western Regional Climate Center (WRCC) data.These atmospheric variables include PSDI, 6-month SPEI, vapor pressure deficit, vapor pressure, mean daily air temperature, snow-water equivalent, maximum and minimum daily temperature, precipitation, and incoming surface shortwave radiation.
The second question focuses on changing patterns of park productivity in recent decades due to climate change and disturbances.In order to identify potential climate change signals, linear trends were computed according to the non-parametric Theil-Sen slope estimator, both spatially by pixel and for the park average.While least squares regression determines the slope based on the weighted mean, the Theil-Sen estimator determines the slope based on the median, thereby making the latter approach more robust against outliers.Linear trends were assessed for significance via the two-sided Mann-Kendall test [36,[110][111][112] for NDVI and EVI from AVHRR and MODIS (as proxies of productivity), WRCC PDSI and SPEI, and precipitation, snow-water equivalent, air temperature, vapor pressure deficit, incoming surface solar radiation, and vapor pressure from Daymet.The Mann-Kendall Trend Test is a non-parametric test that assesses if there is a monotonic upward (consistently increasing, although not necessarily in a linear fashion) or downward trend over time.The Theil-Sen slope estimator and Mann-Kendall test are flexible in that neither make assumptions about the distribution of the analyzed variable.By overlaying shapefiles of the extent of the 1988 fires and more localized fires from 2003,2007,2008,2009, and 2013 onto the map of NDVI trends, changes in productivity were related to post-disturbance forest regrowth, while also examining the undisturbed landscapes for potential responses to greenhouse warming.
In order to ensure the robustness of the study's findings, all correlations and trend analyses were repeated with the rank-based, non-parametric Spearman's rho (ρ) test.This test is advantageous over Pearson's correlations as it operates on the data's ranks, rather than the raw data, is unaffected by the population's distribution, is mostly insensitive to outliers, and is not limited to linear dependencies.However, it has certain disadvantages to the Pearson's correlation, namely, the loss of quantitative information when the data is converted to ranks and less power for normally distributed data.In terms of trend analysis, a trend is considered significant with the Spearman's rho test if there is a significant correlation between the ranks and time-steps.Elevation was negatively correlated with annual mean air temperature and positively correlated with annual mean precipitation, with significant (p < 0.05) Pearson spatial correlations (adjusted N = 17,481 MODIS grid cells) of −0.90 and 0.47, respectively (Spearman's rho  = −0.89and 0.53, respectively).The park contains a sharp north-south climatic gradient (recognized by Whitlock and Bartlein [113] and Bartlein et al. [44]) in annual mean precipitation.As discussed by Kokaly et al. [45], the park's high-elevation temperate forests receive abundant wintertime precipitation, while the valleys are drier and largely support grassland and sagebrush communities.The impressive spatial gradient in annual mean precipitation is characterized by a 6:1 ratio between its wettest and driest points, far exceeding the 2:1 ratio in elevation between the highest and lowest points.The parkaverage precipitation seasonal cycle contains two maxima in November-January and May-June and a minimum in August-October.Mean precipitation was greatest during December across 78% of the park and May across 22% of the park (e.g., Lamar Valley).This is consistent with Despain [42], who Elevation was negatively correlated with annual mean air temperature and positively correlated with annual mean precipitation, with significant (p < 0.05) Pearson spatial correlations (adjusted N = 17,481 MODIS grid cells) of −0.90 and 0.47, respectively (Spearman's rho ρ = −0.89and 0.53, respectively).The park contains a sharp north-south climatic gradient (recognized by Whitlock and Bartlein [113] and Bartlein et al. [44]) in annual mean precipitation.As discussed by Kokaly et al. [45], the park's high-elevation temperate forests receive abundant wintertime precipitation, while the valleys are drier and largely support grassland and sagebrush communities.The impressive spatial gradient in annual mean precipitation is characterized by a 6:1 ratio between its wettest and driest points, far exceeding the 2:1 ratio in elevation between the highest and lowest points.The park-average precipitation seasonal cycle contains two maxima in November-January and May-June and a minimum in August-October.Mean precipitation was greatest during December across 78% of the park and May across 22% of the park (e.g., Lamar Valley).This is consistent with Despain [42], who identified two major climatic types in the park, consisting of a valley type with a spring precipitation peak and a mountain type with a winter precipitation peak.Summer to early-autumn is the dry season, with lowest mean precipitation in July across 74% of the park and August across 21% of the park.

Park Summary Statistics on Vegetation Greenness Patterns
High-resolution (250-m) maps of basic statistics on MODIS-based vegetation phenology indices (VIs), namely, NDVI and EVI, were developed for Yellowstone National Park, including the annual mean and peak-month average during 2000-2017 (Figure 5).Annual mean NDVI (EVI) ranged from −0.04 to 0.59 (−0.03 to 0.33), while peak-month average NDVI (EVI) ranged even more widely from −0.01 to 0.82 (0.00 to 0.65).As seen in Figure 5, peak-month greenness was generally greatest across low-elevational meadows, such as Cascade Corner, Hayden Valley, fringes of Lamar Valley, Pelican Valley (near the north shore of Yellowstone Lake), and the south entrance near Snake River, and least across Madison Plateau, Pitchstone Plateau, Hebgen Basin, and peaks of the Absaroka Range.While low elevation typically supported higher peak greenness in the park, Hebgen Basin was a clear exception to this relationship.That low-elevational basin, known for its diverse array of flowering plants, also exhibited a low amplitude of the seasonal cycle of greenness.Spatial correlations between the NDVI-and EVI-based maps shown in Figure 5 indicate extremely high consistency between the two MODIS VIs for annual mean (Pearson's correlation r = 0.91, p < 0.05; ρ = 0.92) and moderate consistency for peak-month average (r = 0.69, p < 0.05; ρ = 0.65).Furthermore, the NDVI and EVI anomaly time-series were highly correlated with each other across non-forested landscapes and moderately correlated across Douglas fir-dominated landscapes, particularly around Lamar Valley.Spatial correlations between VI and climatic maps indicated significant positive correlations between annual NDVI/EVI and air temperature (NDVI: r = 0.26, ρ = 0.29; EVI: r = 0.31, ρ = 0.34) and significant negative correlations (NDVI: r = −0.28,ρ = −0.27;EVI: r = −0.28,ρ = −0.26) between annual NDVI/EVI and precipitation, all at p < 0.05.This may be interpreted as climatologically warm, dry (cold, wet) regions corresponding to greater (reduced) annual greening.More specifically, springtime (April-June) greening occurs later in cold, wet regions that are typically characterized by an extensive snowpack that must melt prior to vegetation emergence.Overall consistency of results between NDVI and EVI, which are computed using different reflectance channels, enhances the confidence in the study's findings and the capacity of remote-sensing methods to extract spatio-temporal variations in greenness.
regions corresponding to greater (reduced) annual greening.More specifically, springtime (April-June) greening occurs later in cold, wet regions that are typically characterized by an extensive snowpack that must melt prior to vegetation emergence.Overall consistency of results between NDVI and EVI, which are computed using different reflectance channels, enhances the confidence in the study's findings and the capacity of remote-sensing methods to extract spatio-temporal variations in greenness.

Spatio-Temporal Variability in Vegetation Greenness
The spatio-temporal variability in vegetation growth across Yellowstone National Park is largely determined by complex patterns of topography and associated climatological air temperatures and snowpack, as seen in Figure 6.Here, the seasonal cycle of the primary (first) mode of spatio-temporal variability (EOF1) in vegetation greenness was assessed to understand the favored timing and location of high variability in vegetation greenness, and thus phenological activity.This was addressed by applying unrotated EOF analysis to detrended (to remove long-term trends and focus instead on interannual variability) NDVI anomalies during individual calendar months from April through November of 2000-2017 (Figure 6a-h).This dominant mode, EOF1, explains the highest percentage of NDVI variability during spring to early summer, ranging from only 23% in July during the timing of maximum greenness to 66% in June during the timing of peak green-up.For each month, EOF1 was best characterized as a monopole, with the entire park generally exhibiting the same sign anomaly in NDVI.These monopole patterns were largely homogeneous in sign, but not magnitude, with distinct areas of enhanced NDVI variability that largely follow the park's vegetation classes and elevational gradients.Reflecting the influence of elevation on air temperature, and thus the seasonal timing of greening, the peak variance in NDVI anomalies occurred early in April at the milder low elevations of around 2000-2200 m in the Lamar Valley and Hebgen Basin (Figure 6a), in May at moderate elevations of around 2400-2600 m (Figure 6b), and late, in June, at cold high elevations of around 2800-3000 m on Madison and Pitchstone Plateaus and Absaroka Range (Figure 6c), representing interannual variability in the timing of growing season onset (Figure 6i-p); this elevational transition reflects the shift from valley meadows to lodgepole forests to alpine meadows.The timing of green-up is highly dependent on, and positively correlated with, elevation.Green-up is the most elevation-sensitive phenological stage, compared to the phases of maturity and senescence (Figure 6i-p).As evidence, the average EOF1 loading pattern exhibits a distinct elevation-based peak during April (Figure 6i), ranging from <0.0013 for 2800-3200 m to 0.0045 for 2000-2200 m, and June (Figure 6k), ranging from <0.0004 for 1600-2200 m to 0.0043 for 2800-3000 m.In contrast, there is minimal signature of the effect of elevation during August (Figure 6m), when the loading pattern modestly varies between 0.0012 and 0.0027 across elevational bins.(d-f, i-j) time-series of (a,d,g,i) empirical orthogonal function mode 1 (EOF1), (b,e,h,j) EOF2, and (c,f) EOF3 of year-round monthly Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) anomalies across Yellowstone National Park during 2000-2017, both (a-f) without and (g-j) with linear detrending applied to the NDVI data.EOF3 of the detrended NDVI anomalies is excluded here, as it is a noisy, incoherent spatial pattern with trivial explained variance.The explained variance associated with each eigenvalue in EOF1-3 was 31.0%,13.1%, and 7.7%, respectively, for the un-detrended data and 31.8%,9.3%, and 3.4%, respectively, for the detrended data.The explained variance for each EOF mode was computed by dividing that mode's eigenvalue by the sum of all modes' eigenvalues and multiplying by 100%, where an eigenvalue expresses how much of the variance can be explained by its associated eigenvector.
Next, attention shifts to the topic of interannual to decadal variability in vegetation productivity.An EOF analysis was performed on year-round monthly NDVI anomalies during 2000-2017, rather  (d-f,i-j) time-series of (a,d,g,i) empirical orthogonal function mode 1 (EOF1), (b,e,h,j) EOF2, and (c,f) EOF3 of year-round monthly Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) anomalies across Yellowstone National Park during 2000-2017, both (a-f) without and (g-j) with linear detrending applied to the NDVI data.EOF3 of the detrended NDVI anomalies is excluded here, as it is a noisy, incoherent spatial pattern with trivial explained variance.The explained variance associated with each eigenvalue in EOF1-3 was 31.0%,13.1%, and 7.7%, respectively, for the un-detrended data and 31.8%,9.3%, and 3.4%, respectively, for the detrended data.The explained variance for each EOF mode was computed by dividing that mode's eigenvalue by the sum of all modes' eigenvalues and multiplying by 100%, where an eigenvalue expresses how much of the variance can be explained by its associated eigenvector.
Next, attention shifts to the topic of interannual to decadal variability in vegetation productivity.An EOF analysis was performed on year-round monthly NDVI anomalies during 2000-2017, rather than the seasonal cycle time scale of Figure 6.The focus was on the leading three principal components, which explain the greatest percentage of NDVI variance (Figure 7a-f).According to Figure 7, the dominant signal of variability in vegetation greenness during the early 21st century was that of forest regrowth following the 1988 fires.The pattern of EOF1 appears to consist of a largely park-wide re-establishment of productivity since the 1988 fires, although with notable, climate-associated disruptions in October 2004, May-June 2011, May 2008, October 2000, and September 2013 (Figure 7a,d).The decline in productivity associated with fires during 2003-2009 in proximity to Yellowstone Lake, Absaroka Range, and Central Plateau was suggested through EOF2 (Figure 7b,e), with a close correspondence between the shapefiles of the 2003-2009 fire distributions and the dark-brown colors in Figure 7b.The pattern of EOF3 appears to represent contrasting spring-time phenological variability between low-and high-elevational regions, with no trend in the loading pattern (Figure 7c,f).As noted in Figure 6, when the valleys are seasonally greening up, the high mountains remain mostly the same, and vice versa.Repeating the EOF analysis, but for linearly detrended monthly NDVI anomalies (Figure 7g-j), allows for the exclusion of the positive NDVI trend associated with recovery from the 1988 fires.The park-wide phenological responses (mostly uniform in sign) to climate variability are captured by EOF1 of the detrended NDVI anomalies.The second dominant mode of variability of un-detrended NDVI anomalies closely matched the elevation-driven pattern seen in Figure 7c,f.

Climate-Phenology Correlations
The sensitivity of park-averaged vegetation greenness to climatic anomalies was examined through temporal correlations between detrended AVHRR NDVI anomalies and antecedent atmospheric detrended anomalies during 1982-2015 (Figure 8).
As seen in Figure 8, NDVI was most climate sensitive during the green-up period of May-June, particularly to drought indices and vapor pressure deficit, and least sensitive during the peak growing season month of July.As a general metric of the importance of individual atmospheric forcings, the number of significant Pearson's correlation values in the correlation matrices shown in Figure 8, from greatest to least, was 32 (out of 96) for vapor pressure deficit, 29 for PDSI, 25 for SPEI6mn, 25 for precipitation, 22 for maximum temperature, 17 for snow-water equivalent, 17 for mean temperature, 17 for minimum temperature, 9 for vapor pressure, and 7 for solar radiation.As an alternative metric, the peak lag-0 month correlation value in these correlation matrices was greatest for PDSI (r = −0.81 in May), SPEI6mn (r = −0.80 in June), vapor pressure deficit (r = 0.79 in May), mean temperature (r = 0.70 in May), snow-water equivalent (r = −0.69 in May), and maximum temperature (r = 0.68 in June), with generally very low correlation values for solar radiation, vapor pressure, and precipitation.These two metrics suggest that PDSI, vapor pressure deficit, and SPEI6mn, namely, drought measures, were the dominant contributors towards phenological variability in Yellowstone National Park, while solar radiation and vapor pressure played a minimal role.Winter-spring drought, as evident by negative anomalies in PDSI, SPEI (both drought indices take snowfall into account), and vapor pressure deficit, typically corresponds to enhanced greening in April-June, while spring-summer drought was often associated with negative NDVI anomalies, reduced greening, and potential senescence during the dry season months of August-September.The moderate-duration drought index, SPEI6mn, exhibited substantially greater correlations with NDVI anomalies than did the short-duration SPEI1mn, with comparable correlation strengths between SPEI6mn-NDVI and PDSI-NDVI.Temperature was most positively correlated with NDVI and EVI in May-June and most negatively correlated in August.Higher maximum and minimum temperatures were associated with enhanced greenness in March-June and April-June, respectively.An anomalously deep snowpack in February-June delays the onset of greening in April-June.The negative correlation between late spring NDVI and both preceding and concurrent precipitation reflects delayed spring-time growth due to (1) anomalously deep and persistent snowpack from cumulative cold-season precipitation and (2) lower air temperatures generally associated with wet conditions.September is the only month in which precipitation was significantly positively correlated with NDVI.July-September is Yellowstone's dry season, producing only 13% of the annual mean precipitation, so soils were typically driest by September [42] and vegetation becomes moisture-limited.Park-average anomaly time-series of climatic (precipitation, snow-water equivalent, PDSI, air temperature, SPEI, vapor-pressure deficit, incoming surface solar radiation, and vapor pressure) and phenological (NDVI) variables were assessed for 1982-2015 (Figure 9).The tendency towards a warmer climate with diminished snowpack was dramatic in this three-decade period (Figure 9).All of these environmental variables, except for precipitation, experienced a statistically significant trend (p < 0.05) during the 34-year analysis period, based on the Theil-Sen estimate of linear trend and Mann-Kendall trend significance test; the trends were likewise significant (p < 0.05) for all of the variables, except precipitation, according to the Spearman's rho test.Our comprehensive analysis of the park's trends indicates that linear trends were positive for temperature (net change: +2.5 • C), vapor pressure (+0.95 hPa), and vapor pressure deficit (+0.24 hPa), and negative for NDVI (−0.02),SPEI6mn (−0.65), snow-water equivalent (−221.9kg m −2 ), PDSI (−1.42), and solar radiation (−37.0W m −2 ).The warming trend was most pronounced at night, with a net change of +3.4 • C in minimum temperature compared to +1.6 • C in maximum temperature.Mean snow-water equivalent declined abruptly [37,53,114]  By relating the climatic and NDVI time-series in Figure 9, it is evident for largely temperature-limited Yellowstone National Park that cool, wet conditions generally delay vegetation growth (e.g., 2010-2011 La Niña), while warm, dry conditions enhance greenness (e.g., 2014-2015 El Niño).During the MODIS era, there is evidence that the annual growth of the young post-1988 lodgepole pines was delayed by anomalous browning episodes during October 2004, May 2008, and May 2011, each associated with relatively cold conditions.Unlike the delayed spring-time green-up episodes in May of 2008 and 2011, the anomalous late-season browning in October 2004 represents earlier senescence, ending the growing season, due to persistently low maximum temperatures during May-October 2004, including a rapid decline in air temperature and occurrence of snowfall by mid-October.
According to the Theil-Sen slope estimator, positive NDVI (EVI) trends were detected across 75% (83%) of the park, with negative trends found over only 25% (17%) of the region-largely coinciding with recent fire activity.If only statistically significant trends (p < 0.05) are considered based on the Theil-Sen slope estimator and Mann-Kendall test, then positive NDVI (EVI) trends were identified over 13.8% (28.2%) of the park, as opposed to a very limited extent of negative NDVI (EVI) trends over 0.7% (0.9%) of the park.Specifically, within the 1988 fire burn polygons, 92.2% (96.4%) of the landscape experienced positive NDVI (EVI) trends.When the calculations were repeated using the Spearman's rho test, it was determined that positive NDVI (EVI) trends were found across 74% (82%) of the park, with negative trends located over 26% (18%) of the park; furthermore, significant positive NDVI (EVI) trends were detected over 27.4% (41.6%) of the region and significant negative NDVI (EVI) trends were found over 4.0% (3.3%) of the region.
the MODIS era, there is evidence that the annual growth of the young post-1988 lodgepole pines was delayed by anomalous browning episodes during October 2004, May 2008, and May 2011, each associated with relatively cold conditions.Unlike the delayed spring-time green-up episodes in May of 2008 and 2011, the anomalous late-season browning in October 2004 represents earlier senescence, ending the growing season, due to persistently low maximum temperatures during May-October 2004, including a rapid decline in air temperature and occurrence of snowfall by mid-October.The study's spatio-temporal analysis of remotely sensed phenology yielded key insights into Yellowstone National Park's terrestrial ecosystems.Peak seasonal greenness is achieved across low-elevational meadows in Cascade Corner, Hayden Valley, Lamar Valley, Pelican Valley, and Snake River Valley.Green-up is the most elevation-sensitive phenological stage in Yellowstone National Park, with the variance in NDVI anomalies peaking in April at relatively low elevations, May at moderate elevations, and June at high elevations.Snow often covers the highest elevations of Yellowstone through June, delaying the arrival of green-up [53].The finding that the timing of green-up more closely tracks elevational gradients than does the timing of senescence is consistent with studies by Vitasse et al. [115], Hwang et al. [116,117], and Elmore et al. [118].The park achieves its peak MODIS NDVI and EVI in July, which is consistent with the study by Wilmers et al. [48] but slightly earlier than the identified July-August peak by Potter [37] using Landsat NDVI data.
Drought metrics are the dominant influence on phenological variability in Yellowstone National Park.This established importance of drought on Yellowstone's vegetation phenology makes it critical to better understand drought-phenology interactions under future climate change with a changing disturbance regime, as projected amplification of western drought events with anthropogenic climate change [119,120] will likely have critical implications for Yellowstone's vegetation productivity.Drought and high temperatures enhance spring-time greenness and reduce late summer greenness.During anomalously mild springs, snow melt occurs earlier, triggering the green-up wave that propagates from lower to higher elevations over the period of 1-2 months.Given that the soils are initially saturated from snow melt in spring [121], water is likely not limited and warm temperatures can trigger plant growth.Summer drought timing determines when herbaceous vegetation begins to cure and when conifers restrict stomatal conductance and transpiration to conserve water [122].Cool, wet conditions associated with La Niña generally delay green-up.Unseasonably low temperatures resulted in anomalous browning episodes in October 2004, May 2008, and May 2011, thereby slowing the growth of the post-1988 lodgepole pine.Higher temperatures were associated with enhanced greenness in March-June, which has potentially positive implications for climate change with greater forage for ungulates and an enhanced seasonal carbon sink.Moving-window temporal correlations (not shown) hint at an increase in the climatic sensitivity of NDVI in Yellowstone National Park since the 1988 fire, perhaps as younger, shorter trees that have become established are more vulnerable to climatic anomalies [123,124]  The climate change signal has been especially pronounced in recent decades across Yellowstone National Park, including substantial warming, an abrupt decline in snowpack, and more frequent drought episodes.The identified warming trend was consistent with findings from Chang and Hansen [49] and Sepulveda et al. [125].Specifically, Chang and Hansen [49] assessed historical climate change across the broader Yellowstone and Grand Teton National Parks for 1948-2010 and concluded that annual temperatures increased at +0.16 • C decade −1 .The absence of an identified precipitation trend in the current study was also noted by Gray et al. [126], Newman and Watson [127], and Wilmers et al. [48].The warming over recent decades has been associated with earlier spring snowmelt; expanded growing and fire seasons; and forest pest outbreaks and resulting tree mortality [2,[48][49][50]128,129].Consistent with the current study, Tercek et al. [53], Potter [37], and Tercek and Rodman [114] likewise reported observed declines in Yellowstone snowpack.
The AVHRR NDVI time-series was characterized by mostly positive anomalies during 1983-1987, record positive anomalies in May 1987, a vegetation collapse between 1987 and 1988 corresponding to the 1988 fires, largely negative anomalies during 1988-2011 corresponding to bark beetle outbreaks in 2008-2009, and peak whitebark pine mortality in 2008-2011, record negative anomalies in June 2011 (likely in response to persistent anomalously cold conditions during February-June 2011), and a recovery towards positive anomalies during 2012-2015 (when 73% of months averaged above-normal in temperature) (Figure 9).Throughout 1982-2016, Yellowstone experienced the greatest area burned in 1988, followed by 2016 (specific to Hebgen Basin), both during low-precipitation years.Franks et al. [33] noted that the summer of 1988 was the driest on record since 1886, although the region experienced a near-average preceding snow year [30].Consistent with the findings of Changnon [130], the present study found that the extreme 1988 fire season included the lowest PDSI value since 1948, −5.39 in October.Since the 1988 fires, lodgepole pines established and have been growing rapidly [68,[131][132][133], with an associated increase in productivity and biomass [132,134].The re-establishment of vegetation since the 1988 fires has been rapid yet spatially variable, with stand structure and function not yet converged [132,135].The regrowth is clearly evident in the predominantly positive productivity trends during 2001-2016.Note the close correspondence in Figure 10 between the spatial extent of the 1988 fire perimeters and areas of positive productivity trends, as evidence of continued re-establishment of the burned forest lands [136].Spatio-temporal variability in fire disturbances in the park incurs a distinct impact on trends in forest productivity (Figure 10).nine-point temporal smoother was applied to the time-series in (a-i).Time-series of monthly anomalies across the park during 2000-2017 in Moderate Resolution Imaging Spectroradiometer (MODIS)-based (j) NDVI and (k) Enhanced Vegetation Index (EVI).
According to the Theil-Sen slope estimator, positive NDVI (EVI) trends were detected across 75% (83%) of the park, with negative trends found over only 25% (17%) of the region-largely coinciding with recent fire activity.If only statistically significant trends (p < 0.05) are considered based on the Theil-Sen slope estimator and Mann-Kendall test, then positive NDVI (EVI) trends were identified over 13.8% (28.2%) of the park, as opposed to a very limited extent of negative NDVI (EVI) trends over 0.7% (0.9%) of the park.Specifically, within the 1988 fire burn polygons, 92.2% (96.4%) of the landscape experienced positive NDVI (EVI) trends.When the calculations were repeated using the Spearman's rho test, it was determined that positive NDVI (EVI) trends were found across 74% (82%) of the park, with negative trends located over 26% (18%) of the park; furthermore, significant positive NDVI (EVI) trends were detected over 27.4% (41.6%) of the region and significant negative NDVI (EVI) trends were found over 4.0% (3.3%) of the region.

Focal Question 1: What are the Dominant Meteorological Drivers of Variability in Vegetation Greenness across Yellowstone National Park on the Seasonal to Interannual Time Scales?
The study's spatio-temporal analysis of remotely sensed phenology yielded key insights into Yellowstone National Park's terrestrial ecosystems.Peak seasonal greenness is achieved across low- While Yellowstone National Park has exhibited predominantly positive productivity trends across lodgepole pine forests and non-forested landscapes, associated with re-establishment following the 1988 widespread fires, negative trends were identified in more recently disturbed landscapes in proximity to the Absaroka Range, Central Plateau, and Yellowstone Lake associated with fires in 2003, 2007, 2008, 2009, and 2013.Within the fire burn polygons for these five active fire years, 97.6% (64.6%) of the landscape experienced negative NDVI (EVI) trends.A negative NDVI trend for higher elevation areas dominated by whitebark pine was likely triggered by rising temperatures and resulting spread of blister rust and mountain pine beetle vulnerability [37,60,61,[137][138][139][140].Whitebark pine is a key food source for grizzly bears, thus partly determining suitable grizzly bear habitat [141], and its presence increases biodiversity in animal and plant communities [142].The survival of whitebark pine is threatened by infestations of mountain pine beetles and spread of whitepine blister rust [143], warranting continued monitoring through remote sensing, both in Yellowstone National Park and other additional, less-managed landscapes.
While the present study does not rigorously attribute identified NDVI trends to disturbances, it is complementary to the statistical analysis of Emmett et al. [36] which focused on vegetation greenness trends across the Greater Yellowstone Ecosystem and the relative importance of disturbance and climate in explaining NDVI trends.The current study's findings were generally consistent with those of Emmett et al. [36], who focused on trends in annual maximum NDVI during 1989-2014 across the Greater Yellowstone Ecosystem using USGS EROS NDVI data.Their study identified statistically significant (p < 0.05) positive and negative trends across 26.5% and 6.2% of the region, respectively, and linked these NDVI trends to disturbances and changes in annual temperature and summer precipitation.
Based on a separate analysis (not shown) of park regions which were largely undisturbed by the 1988 fires, there was no evidence during the MODIS era of the direct temperature effect of enhanced greenhouse effect on greening across these undisturbed regions, consistent with Potter [37].This can be partly understood by the fact that the park is dominated by lodgepole pine, which is found to be a generally temperature-insensitive species based on correlations between Daymet air temperature and MODIS NDVI/EVI across lodgepole pine-dominated landscapes in the park (not shown here).Conifer forests grow rapidly during the window up to about 90 years following disturbance.Disturbance resets the stands and initiates new periods of biomass aggradation.During a short window of time, it is expected that mature conifer forests would not notably respond to a warming trend.

Comparison of Results with Past Studies
A further comparison of the present study's findings and those of past studies is presented in Table 1, while noting that few studies have extensively analyzed remotely sensed VIs across Yellowstone National Park.The results were largely consistent, with disparities often associated with differences in analyzed time periods between studies.

Next Research Steps
Several Yellowstone-specific questions emerge in this study, warranting further investigation.(1) How does El Niño-Southern Oscillation affect spatio-temporal variability in precipitation, temperature, and greenness across the park? (2) Why is NDVI more sensitive to climate variability than EVI across Yellowstone National Park?(3) Why is the timing of green-up much more elevationally-sensitive than the timing of senescence?(4) How might future climate change affect the EOF-established spatio-temporal patterns of variability in greenness across Yellowstone National Park as phenological timing shifts earlier?(5) How do the individual tree species in the park differ in terms of sensitivity to variability in temperature, precipitation, snowpack, drought, and snowpack?

Conclusions
Spatio-temporal variability in climatic drivers and terrestrial phenology are examined across the complex topographic, climatic, and ecological gradients of Yellowstone National Park using AVHRR and MODIS VIs and Daymet climate data, aimed at demonstrating the value of remote sensing in exploring ecological sensitivity to climate variability and change.The dominant signals in Yellowstone's remotely sensed VI record are long-term, post-disturbance regrowth trends, particularly seen in lodgepole pine landscapes since the 1988 fires, and short-term episodes of delayed spring-time green-up/anomalous browning associated with exceptionally cool, snowy winters-springs.In the absence of extended phenological records in the park, remote sensing is valuable in exploring vegetation's responses to climatic forcings, and beyond what is demonstrated here, specific species' vulnerabilities and potential refugia, yet remote sensing has been underutilized to deepen ecological knowledge within Yellowstone National Park and most other parks.While the study focuses on a single United States national park, it demonstrates the value of terrestrial ecosystem remote sensing and motivates a systematic application across all major western national parks to aid in assessing ecological vulnerabilities to climate variability and change, and guiding future adaptation planning.
The current study supports the following recommendations for future exploration.
(1) The general dearth of continuous, long-term ecological measurements across Yellowstone National Park serves as a call to action, particularly given its long history of 147 years as an official national park, vast spatial extent, and priceless environmental and societal value.
There is a clear need for new and expanded, long-term initiatives across Yellowstone National Park, and likely many other national parks, for citizen science phenology data collection, installation of pheno-cams and soil moisture probes, and ground measurements of LAI.Such observations are valuable for identifying species' sensitivities to climate variability and change, elucidating remotely sensed phenological changes across diverse landscapes, and assessing the relative reliability of individual satellite-derived VIs.Such observational data needs are further motivated by concerns regarding the reliability of MODIS and Landsat NDVI and EVI across Yellowstone National Park from localized studies by Franks et al. [33], Garroutte [24], and Garroutte et al. [38].Further studies of Yellowstone National Park should consider multiple available gridded meteorological datasets, noting that Daymet is characterized by small biases in temperature, but large biases in precipitation across the conterminous United States according to Behnke et al. [144].It is also noted that all gridded solar radiation datasets, including Daymet, exhibit substantial seasonal biases that may affect regional-scale findings [145].(2) Remote sensing has been widely underutilized across the national park system, including Yellowstone, especially given the aforementioned deficiency in continuous ground-based phenological observations.The importance of monitoring protected landscapes, including national parks, through remote sensing was stressed by Pettorelli et al. [146] and Gillespie et al. [27] and is further emphasized by the present study.Specifically, there is a clear need for a systematic, remote-sensing assessment of species' and landscape vulnerabilities to climate variability and change, including ecological droughts operating on a spectrum of time scales, across the Western United States' national parks and the identification of potential drought refugia.Such insights would be highly valuable to park resource managers, aiming to adapt to a highly variable and rapidly changing climate.Further remote-sensing analyses should also expand to consider higher resolution Landsat and Visible Infrared Imaging Radiometer Suite (VIRRS) data and the potential value of solar-induced chlorophyll fluorescence [147-151] and MODIS-based chlorophyll/carotenoid indices of foliar pigment levels for tracking evergreen photosynthesis, which is a challenge with the standard NDVI [152].

Figure 1 .
Figure 1.Elevation (m) across Yellowstone National Park based on the Parameter-elevation Regressions on Independent Slopes Model (PRISM) 800-m digital elevation model data.Key park regions are identified in text.Lakes are displayed in blue.

Figure 1 .
Figure 1.Elevation (m) across Yellowstone National Park based on the Parameter-elevation Regressions on Independent Slopes Model (PRISM) 800-m digital elevation model data.Key park regions are identified in text.Lakes are displayed in blue.

6 Figure 2 .
Figure 2. (a) Primary vegetation type distribution based on the National Park Service's (NPS) Yellowstone cover types data for 1999, following the species' color scheme in (b).(b) Percent abundance (bar chart) of each vegetation type across the park and the mean elevation (m) for each type provided as text, along with the 5th and 95th percentiles of elevation in parentheses.

Figure 2 .
Figure 2. (a) Primary vegetation type distribution based on the National Park Service's (NPS) Yellowstone cover types data for 1999, following the species' color scheme in (b).(b) Percent abundance (bar chart) of each vegetation type across the park and the mean elevation (m) for each type provided as text, along with the 5th and 95th percentiles of elevation in parentheses.

Figure 3 .
Figure 3. Research schema, summarizing the method applied in the current study.

1 .QUESTION 2 :
Focal Question 1: What are the Dominant Meteorological Drivers of Variability in Vegetation Greenness across Yellowstone National Park on the Seasonal to Interannual Time Scales? 3.1.1.Climatic Heterogeneity across the Park QUESTION 1: What are the dominant meteorological drivers of variability in vegetation greenness across Yellowstone National Park on the seasonal to interannual time scales?How have patterns of productivity changed across the park in recent decades in response to climate change and disturbances?Method: spatial climate composites Data source: Daymet Question: Which climatic variables exhibit the most pronounced spatial heterogeneity across the park?Method: spatial VI composites, spatial correlations and semivariograms Data source: MODIS, Daymet, WRCC Question: How do the mean and variability in vegetation greenness vary across the park, and which topographic and climatic factors govern this spatial heterogeneity in growth?Method: unrotated EOF analysis of VIs Data source: MODIS, Daymet, WRCC Question: What are the dominant modes of variability in vegetation phenological growth, their patterns of seasonality, and their sensitivity to elevation and climate on the seasonal to interannual time scales?Method: lagged Pearson temporal correlation analysis Data source: AVHRR, Daymet, WRCC Question: What are the primary climatic regulators of variability in vegetation greenness, by calendar month?Method: time series analysis; Theil-Sen trend analysis and Mann-Kendall test Data source: AVHRR, MODIS, Daymet, WRCC Question: What signs of climate change impacts are detectable in recent decades across the park?Do trends in productivity reflect climate change or disturbances?

Figure 3 .
Figure 3. Research schema, summarizing the method applied in the current study.

1 .
Focal Question 1: What are the Dominant Meteorological Drivers of Variability in Vegetation Greenness across Yellowstone National Park on the Seasonal to Interannual Time Scales? 3.1.1.Climatic Heterogeneity across the Park Sharp elevational gradients yield dramatic spatial heterogeneity in the climate of Yellowstone National Park (Figure 4).Four climatological variables with noteworthy spatial heterogeneity across Yellowstone National Park include summer (June-August) maximum temperature (ranging from 13.7 • C to 28.2 • C), winter (December-February) precipitation (ranging from 40 mm to 646 mm), annual incoming surface solar radiation (ranging from 307 W m −2 to 443 W m −2 ), and annual maximum snow water equivalent (ranging from 22 kg m −2 to 879 kg m −2 ), based on Daymet data for 1982-2015 (Figure 4).

Figure 6 .
Figure 6.(a-h) Maps of the leading mode (EOF1) of Moderate Resolution Imaging Spectroradiometer (MODIS) detrended Normalized Difference Vegetation Index (NDVI) variability for April-November from 2000-2017 based on empirical orthogonal function (EOF) analysis for Yellowstone National Park.The percent explained variance in NDVI anomalies is provided by each map.(i-p) Average EOF1 loading pattern by elevation (in 200-m bins), with green shading indicating the preferred elevation for specific patterns of variability shown in (a-h).

Figure 6 .
Figure 6.(a-h) Maps of the leading mode (EOF1) of Moderate Resolution Imaging Spectroradiometer (MODIS) detrended Normalized Difference Vegetation Index (NDVI) variability for April-November from 2000-2017 based on empirical orthogonal function (EOF) analysis for Yellowstone National Park.The percent explained variance in NDVI anomalies is provided by each map.(i-p) Average EOF1 loading pattern by elevation (in 200-m bins), with green shading indicating the preferred elevation for specific patterns of variability shown in (a-h).

Figure 7 .
Figure 7. (a-c,g-h) Spatial pattern and(d-f,i-j) time-series of (a,d,g,i) empirical orthogonal function mode 1 (EOF1), (b,e,h,j) EOF2, and (c,f) EOF3 of year-round monthly Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) anomalies across Yellowstone National Park during 2000-2017, both (a-f) without and (g-j) with linear detrending applied to the NDVI data.EOF3 of the detrended NDVI anomalies is excluded here, as it is a noisy, incoherent spatial pattern with trivial explained variance.The explained variance associated with each eigenvalue in EOF1-3 was 31.0%,13.1%, and 7.7%, respectively, for the un-detrended data and 31.8%,9.3%, and 3.4%, respectively, for the detrended data.The explained variance for each EOF mode was computed by dividing that mode's eigenvalue by the sum of all modes' eigenvalues and multiplying by 100%, where an eigenvalue expresses how much of the variance can be explained by its associated eigenvector.

Figure 8 .
Figure 8. Pearson's temporal correlations between detrended monthly Advanced Very High-Resolution Radiometer (AVHRR) Normalized Difference Vegetation Index (NDVI) anomalies, averaged across Yellowstone National Park, and park-mean antecedent atmospheric detrended monthly anomalies (leading NDVI by 0 to 7 months) from Daymet and Western Regional Climate Center (WRCC) during 1982-2015.Only statistically significant correlations (p < 0.05) are displayed.Correlations are shown for (a) vapor pressure deficit, (b) Palmer Drought Severity Index (PDSI), (c) 6-month Standardized Precipitation Evapotranspiration Index (SPEI), (d) precipitation, (e) maximum daily temperature, (f) snow-water equivalent, (g) mean temperature, (h) minimum daily temperature, (i) vapor pressure, and (j) incoming surface shortwave radiation.For each panel, the number of significant correlation boxes (out of 96) is provided.For example, in column 5 of panel (f), the park's average NDVI anomalies in May exhibit a statistically significant negative correlation with

Figure 8 .
Figure 8. Pearson's temporal correlations between detrended monthly Advanced Very High-Resolution Radiometer (AVHRR) Normalized Difference Vegetation Index (NDVI) anomalies, averaged across Yellowstone National Park, and park-mean antecedent atmospheric detrended monthly anomalies (leading NDVI by 0 to 7 months) from Daymet and Western Regional Climate Center (WRCC) during 1982-2015.Only statistically significant correlations (p < 0.05) are displayed.Correlations are shown for (a) vapor pressure deficit, (b) Palmer Drought Severity Index (PDSI), (c) 6-month StandardizedPrecipitation Evapotranspiration Index (SPEI), (d) precipitation, (e) maximum daily temperature, (f) snow-water equivalent, (g) mean temperature, (h) minimum daily temperature, (i) vapor pressure, and (j) incoming surface shortwave radiation.For each panel, the number of significant correlation boxes (out of 96) is provided.For example, in column 5 of panel (f), the park's average NDVI anomalies in May exhibit a statistically significant negative correlation with snow-water equivalent from the preceding to concurrent period of January-May (as shown by the blue/purple boxes under lag 0 to lag −4).Dots indicate significant Spearman's rho (p < 0.05) as a further test of robustness, beyond the Pearson's correlations.The number of significant boxes for the Spearman's rho is25,27,24,22,20,18,17,16,9, and 8 for (a) through (j).

3. 2 .
Focal Question 2: How Have Patterns of Productivity Changed across the Park in Recent Decades in Response to Climate Change and Disturbances?
from 3453 kg m −2 during 1982-1999 to 1883 kg m −2 during 2000-2015, representing a reduction of 45%.Some of the trends during 1982-2015, such as declining solar radiation (likely linked to greater cloud cover or optical depth), may reflect the long-term shift in the Pacific Decadal Oscillation.The MODIS era was dominated by drought across Yellowstone National Park, with 55 months during 2000-2016 characterized by moderate to extreme drought and only 11 months with unusual to extreme moist conditions.The extreme 1988 fire season included the lowest PDSI value on record (since 1948), −5.39 in October and a record-low (since 1982) vapor pressure anomaly of −4.14 hPa in August.

4 . Discussion 4 . 1 .
Focal Question 1: What are the Dominant Meteorological Drivers of Variability in Vegetation Greenness across Yellowstone National Park on the Seasonal to Interannual Time Scales? .

4. 2 .
Focal Question 2: How have Patterns of Productivity Changed across the Park in Recent Decades in Response to Climate Change and Disturbances?

Figure 10 .
Figure 10.Change (shading) in Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI)-based productivity for May-September of 2001-2016, based on linear trend from the Theil-Sen slope estimator, with overlying red polygons indicating the spatial extent of the 1988 widespread fires and overlying blue polygons indicating the extent of the 2003, 2007, 2008, 2009, and 2013 localized fires according to fire perimeter data from the Yellowstone National Park Spatial Analysis Center.The negative NDVI trends in proximity to Yellowstone Lake are associated with fires in 2003, 2007, 2008, 2009, and 2013, as labeled.

Figure 10 .
Figure 10.Change (shading) in Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI)-based productivity for May-September of 2001-2016, based on linear trend from the Theil-Sen slope estimator, with overlying red polygons indicating the spatial extent of the 1988 widespread fires and overlying blue polygons indicating the extent of the 2003, 2007, 2008, 2009, and 2013 localized fires according to fire perimeter data from the Yellowstone National Park Spatial Analysis Center.The negative NDVI trends in proximity to Yellowstone Lake are associated with fires in 2003, 2007, 2008, 2009, and 2013, as labeled.

Table 1 .
Summary of nine papers on variability and change in Yellowstone's vegetation and climate, including the time period, applied datasets, key results, and comparison of results with the current study.
The following are available online at http://www.mdpi.com/2072-4292/11/7/798/s1, Figure S1: Seasonal cycle of quality assurance (QA) indices for MODIS vegetation indices (VI) across Yellowstone during 2000-2017, shown as percentage of data fitting into each color-coded category.Author Contributions: M.N. designed the project objectives, performed all analyses, and wrote the manuscript.K.E. and D.O. provided substantial insights and edits.The research performed by Michael Notaro was funded through the Fall Research Competition, as operated by the University of Wisconsin-Madison Office of the Vice Chancellor for Research and Graduate Education and supported by the Wisconsin Alumni Research Foundation.Kristen Emmett was supported by the National Science Foundation (NSF) Graduate Research Fellowship Program (GRFP) under grant number DGE-1632134.Donal O'Leary was supported by the NSF GRFP under grant number DGE-1322106.Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of NSF.