Disentangling Climate and Disturbance Effects on Regional Vegetation Greening Trends

Productivity of northern latitude forests is an important driver of the terrestrial carbon cycle and is already responding to climate change. Studies of the satellite-derived Normalized Difference Vegetation Index (NDVI) for northern latitudes indicate recent changes in plant productivity. These detected greening and browning trends are often attributed to a lengthening of the growing season from warming temperatures. Yet, disturbance-recovery dynamics are strong drivers of productivity and can mask direct effects of climate change. Here, we analyze 1-km resolution NDVI data from 1989 to 2014 for the northern latitude forests of the Greater Yellowstone Ecosystem for changes in plant productivity to address the following questions: (1) To what degree has greening taken place in the GYE over the past three decades? and (2) What is the relative importance of disturbance and climate in explaining NDVI trends? We found that the spatial extents of statistically significant productivity trends were limited to local greening and browning areas. Disturbance history, predominately fire disturbance, was a major driver of these detected NDVI trends. After accounting for fire-, insect-, and human-caused disturbances, increasing productivity trends remained. Productivity of northern latitude forests is generally considered temperature-limited; yet, we found that precipitation was a key driver of greening in the GYE.


INTRODUCTION
Increases in photosynthetic activity detected in northern latitudes are often attributed to a lengthening of the growing season from warming temperatures (Myneni and others 1997;Tucker and others 2001;Nemani and others 2003;Zhu and others 2013). Evidence suggests anthropogenic warming is also lengthening the fire season, especially in northern latitudes and resulting in increased forest fire area (Flannigan and others 2013;Abatzoglou and Williams 2016). Since wildland fire plays a major role in temperate and boreal forest dynamics, research on the recent changes in photosynthetic activity must consider disturbance recovery and potential changes in disturbance dynamics (Bond-Lamberty and others 2007; Sulla-Menashe and others 2018). Understanding vegetation dynamics in northern latitudes can aid in anticipating vulnerability of forests to climate change and impacts on terrestrial carbon cycling (Pan and others 2011;Le Qué ré and others 2017). The increased availability of satellite-derived observations of Normalized Difference Vegetation Index (NDVI), Leaf Area Index (LAI), and other vegetation indices has led to many recent global and regional studies of plant productivity others 2013, 2016;Potter 2015;Ju and Masek 2016;Rafique and others 2016;Sulla-Menashe and others 2018). Yet, the spatial variability in changes in plant productivity and its drivers are not well understood (Wang and others 2011;Zhu and others 2016).
Plant growth in northern latitudes is considered primarily limited by temperature and the length of the growing season (Whittaker 1975;Churkina and Running 1998). Evidence of spring time warming in the northern hemisphere is well documented and associated with detected increases in photosynthetic activity, or vegetation ''greening'' trends (Myneni and others 1997;Tucker and others 2001;Zhou 2001;Nemani and others 2003;Wang and others 2011;Xu and others 2013). For the northern hemisphere, the growing season onset is estimated to be 1.1-3.3 days per decade earlier (Wolkovich and others 2012). However, changes in summer precipitation and increases in water vapor pressure deficit (VPD) were associated with decreases in photosynthetic activity, or vegetation ''browning'' trends (Nemani and others 2003;Angert and others 2005;Piao and others 2011;Zhu and others 2013). Vegetation responses are also dependent on the bioclimatic context. Warming can cause greening trends in ecosystems that are not water limited, while also causing lowered productivity in drier areas where plants experience increased water stress (Sulla-Menashe and others 2018). These findings and recent tree-ring research support the idea that with warming temperatures, presumed temperature-limited boreal vegetation dynamics are becoming limited by other climatic factors (Angert and others 2005;others 2008, 2009;Williams and others 2011).
In contrast to these direct climatic drivers, growing evidence suggests that detected changes in plant productivity can also be attributed to disturbance recovery, especially post-fire regeneration (Wang and others 2011;Ju and Masek 2016;Sulla-Menashe and others 2018). Ju and Masek (2016) found that most of the greening and browning trends in areas of the eastern boreal forests of Alaska and Canada were related to the fire history. In a follow-up study, Sulla-Menashe and others (2018) determined that fire effects explained the magnitude and direction of NDVI changes, with results dependent on the stage of the disturbancesuccession cycle: greening trends were detected in areas that were regenerating after older fires (prior to 1990), whereas recently burned areas displayed browning trends (Sulla-Menashe and others 2018). Increases in the extent and severity of forest insect outbreaks, including mountain pine beetle (Dendroctonus ponderosae), are attributed to warming temperatures that expand suitable habitat and increase population growth for insects (Logan and others 2010;Coops and others 2012;Weed and others 2013). A severe mortality event following an outbreak can affect productivity for decades (Weed and others 2013). Extensive human-caused disturbances, including logging, also affect productivity trends (Zhu and others 2016). As with fire, for insect-and human-caused disturbances, detections of browning or greening trends depend on the time since the event and the productivity of secondary vegetation (Hicke and others 2012;Zhu and others 2016). Herbaceous vegetation initially dominates after high mortality disturbance events until the newly established tree seedlings grow to dominate the canopy. Because NDVI measures the overall amount of light absorbed and reflected by vegetation, it will indicate recovery of any vegetation cover after a disturbance, capturing the initial recovery of herbaceous vegetation with less sensitivity to later successional changes in forest structure (Franks and others 2013). Invasion by opportunistic exotic herbaceous species post-disturbance would also be detected as ''early'' recovery in NDVI.
Many studies evidence the rapid regeneration of forests post-disturbance in the U.S. Rocky Mountain forests of the Greater Yellowstone Ecosystem (GYE) others 2004, 2016;others 2005, 2013;Zhao and others 2016). Within the GYE, productivity and biomass increased rapidly following the large fires of 1988 others 2004, 2016). In a 300-year chronosequence Kashian and others (2013) found that carbon accumulation recovered to 80% of pre-fire carbon storage in only 50 years. Harvested areas in the region recovered even more quickly (50-90% forest spectral recovery) than burned areas (< 40% forest spectral recovery) since the 1980s (Zhao and others 2016). However, forest regeneration rates can vary widely, with post-fire stem density varying from zero to 344,067 stems ha -1 24-years postfire (Turner and others 2004). After accounting for differences in density, Turner and others (2004) found above-ground net primary productivity decreased with increasing elevation, primarily due to differences in climate. In a study of mountain pine beetles effects on forest productivity in the GYE, Romme and others (1986) found that productivity recovered to pre-outbreak levels in only 11-15 years, mostly due to increased growth in the understory. In Douglas fir-dominated forests in the GYE, understory cover was 24% greater in areas affected by high mortality bark beetle events 4-5 years prior than compared to undisturbed sites due to greater forb and grass cover (Griffin and Turner 2012).
The detection of local greening and browning areas is dependent on the spatial resolution of satellite data used for the study. Regional studies of recent changes in plant productivity show that trends are heterogeneous across landscapes, with some local greening, browning, and areas without a detectable change (Stow and others 2007;Forkel and others 2013;Sulla-Menashe and others 2018). In a comparison of a Local Area Coverage (1-km resolution) to a Global Area Coverage ($ 8-km resolution) NDVI dataset, many greening locations were lost in the coarser resolution dataset (Stow and others 2007). A study using 30-m resolution Landsat NDVI data also shows discrepancies in trend detection compared to 8-km AVHRR (Advanced Very High Resolution Radiometer) NDVI dataset (Ju and Masek 2016). In addition, gridded climate datasets, used as covariates to explain greening trends, are interpolated from weather station data that are only available at moderate spatial resolutions (800 m, 1 or 4 km).
Here we investigate potential greening trends and their drivers in an analysis of recent NDVI trends in the Greater Yellowstone Ecosystem (GYE) of the U.S. northern Rocky Mountains, a large intact, fire-prone temperate forested ecosystem. Our analyses addressed the following questions: (1) To what degree has greening taken place in the GYE over the past three decades? and (2) What is the relative importance of disturbance and climate in explaining NDVI trends? To address these questions, we analyzed 25 years  of vegetation greening (increasing NDVI) and browning (decreasing NDVI) within the GYE using 1-km AVHRR NDVI Composites and compared detected NDVI trends to meteorological data at the same resolution, while using Landsat 30-m resolution data to more precisely attribute disturbance histo-ries (Eidenshink 1992(Eidenshink , 2006Goward and others 2015;Thornton and others 2017).

STUDY AREA
The Greater Yellowstone Ecosystem (GYE) is considered one of the most intact ecosystems in the contiguous USA and encompasses parts of four national forests and two national parks (Yellowstone and Grand Teton National Park) (Parmenter and others 2003) (Figure 1). The climate of the GYE is highly influenced by the U.S. northern Rocky Mountain topography, with over 70% of the region at or above 1800 m (calculated from elevation data: Lehner and others 2008). Although ecologically considered a predominately temperate biome (according to Olson and others (2001) world map), three Kö ppen-Geiger climate classifications describe the majority of the GYE: boreal warm summer on the western side (36% of the GYE), boreal cool summer in the central more mountainous terrain (32%), and on the eastern side due to the rain shadow effect, arid cold steppe (28%) (calculated from Kö ppen-Geiger data: (Kottek and others 2006).
In western North America, cold winters and dry summers colimit vegetation productivity (Churkina and Running 1998;Nemani and others 2003). Within the GYE winters are long and harsh, with the growing season as short as two months in higher elevations and further limited by summer drying (Despain 1990). Fire also impacts vegetation patterns, especially in the GYE (Whitlock and others 2003). The historical range of variation of fire disturbance in the GYE is dependent on the vegetation type, climate, and elevation (Schoennagel and others 2004). Subalpine forests occupying boreal cool summer areas are characterized by infrequent (fire-return intervals of 150-300 years) high-severity fires (Romme and Despain 1989; Higuera and others 2010). Lodgepole pine (Pinus contorta) dominates much of the subalpine, with subalpine fir (Abies lasiocarpa) and Engelmann spruce (Picea engelmannii) as the late seral species, and whitebark pine (Pinus albicaulis) dominating at higher elevations. Mid-elevation mixed conifer forests consisting of predominately Douglas-fir (Pseudotsuga menziesii), with lodgepole pine, aspen (Populus tremuloides), and limber pine (Pinus flexilis), are characterized by a boreal warm summer climate and mixed-severity fire regimes, with infrequent high-and low-severity fires (Schoennagel and others 2004). Lower-elevation sagebrush steppe and shrublands are characterized by an arid cold climate with fire-return intervals ranging from Disentangling Climate and Disturbance Effects decades to centuries (Whisenant 1990;Bukowski and Baker 2013).
Analysis of gridded climate datasets for the GYE indicates annual temperatures have already increased by 0.02°C y -1 from 1928 to 2010 for the region (Chang and Hansen 2015). Despite annual precipitation increasing 0.48 mm y -1 (0.8%) for the same time period (Chang and Hansen 2015), stream discharge decreased, especially during summer months in the GYE (Leppi and others 2012). Leppi and others (2012) explain this paradox by suggesting air temperatures and evapotranspiration rates were stronger drivers of summer stream discharge than precipitation implying increased plant water stress via increases in VPD. Future climate change is expected to shift vegetation growth in northern latitude forests from being temperature to precipitation limited and drastically increase fire severity, fire occurrence, and area burned (Westerling and others 2011; Williams and others 2011).

DATA AND METHODS
For this study we identified the boundaries of GYE based on 10-digit hydrologic units as defined by the USGS Watershed Boundary Dataset (Figure 1). NDVI and covariate datasets were clipped to the GYE boundary and analyses carried out in the World Geodetic System (WGS 84) coordinate reference system. All analyses were carried out using the R statistical program version 3.4.3 (R Core Team 2017) or by using QGIS geoprocessing tools (QGIS Development Team).

Satellite-Derived Vegetation Index
To track changes in vegetation productivity, we used remotely sensed NDVI data. NDVI datasets derived from AVHRR sensors are the most widely used proxy for observing temporal and spatial changes in terrestrial plant productivity and its response to climate (Running 1990;Pettorelli and others 2005;Beck and others 2011). Satellite sensors measure solar radiation reflected by vegetation and NDVI is calculated from the red (RED, 0.58-0.68 lm) and near-infrared (NIR, 0.725-1.1 lm) reflected light channels (Tucker 1979;Running 1990;Myneni and others 1995) using the formula in equation 1: The chlorophyll in plants absorbs visible light including the red light spectrum and the mesophyll leaf tissue reflects near-infrared light in relation to plant water status (Myneni and others 1995). Thereby, generally the more dense the vegetation the greater the difference between near-infrared light reflected and visible light reflected, since the proportions of both the visible light absorbed and the near-infrared light reflected are greater. NDVI values range from -1 to + 1; values near zero indicate sparse vegetation; and values near + 1 indicate dense vegetation (Tucker 1979;Running 1990). Whereas single NDVI measures serve as ''snapshots'' of forest cover, analysis of changes in NDVI over time can indicate changes in productivity. The U.S. Geological Survey's (USGS) Earth Resources Observation and Science (EROS) Center used daily observations from Advanced Very High Resolution Radiometer (AVHRR) on National Oceanic and Atmospheric Administration (NOAA) satellites 11 (1989-1994), 14 (1995-2000), 16 (2001-2003), 17 (2004-2009), 18 (2010-2011), and 19 (2012-2014) to produce biweekly composites of maximum NDVI (Eidenshink 1992(Eidenshink , 2006. Time variant calibration coefficients for intra-and inter-sensor prelaunch and/or postlaunch calibration were derived using a piecewise linear fit of observations over desert, ocean or cloud observations areas (Vermote and Kaufman 1995;Kaufmann and others 2000). The USGS created the maximum NDVI biweekly composites for the conterminous USA from 1989 to 2015 at 1-km resolution (Eidenshink 1992(Eidenshink , 2006. We used the USGS EROS NDVI biweekly composites (data available from the U.S. Geological Survey), because we were interested in investigating moderate spatial resolution (1 km) landscape heterogeneity. The benefit of using a higher-resolution dataset is that localized areas of greening can be easily detected that may be missed in coarser datasets (Stow and others 2007). Data of each year were checked for missing, incomplete, or contaminated biweekly composites (from http://lpdaac.usg s.gov) ( Figure 2). Gaps were infilled by values from a fitted harmonic sinusoidal model (Verbesselt and others 2010). We created monthly maximums by weighting the biweekly NDVI composites by the number of days they covered within a particular calendar month. Monthly NDVI maximums were used to detect seasonal and annual trends to reduce the influence of incomplete or missing biweekly data. Maximum NDVI value compositing is a common preprocessing technique to ''unmask the vegetation signal'' by reducing cloud contamination, shadow effects, and water-vapor effects (Holben 1986;Myneni and others 1995). Annual maximum and integrated NDVI were calculated from the monthly maximum NDVI values. Breakpoint analysis on the NDVI time series was completed for each pixel in the GYE to test for artifacts of satellite sensor changes using the ''breakpoints'' function in the ''strucchange'' R package others 2002, 2003;Bai and Perron 2003); no systematic breaks were detected in the time series.
NDVI saturates in high biomass areas, but this was not an issue in the region since values over the saturation point of 0.9 NDVI (Huete and others 2002) accounted for 0.015% of the dataset.

Disturbance, Meteorological, and Land Cover Data
Disturbance, meteorological, and land cover datasets were considered as potential covariates in both filtering the NDVI dataset for spatial analysis and statistical modeling of NDVI trends (Appendix D: Table D1 lists all variables considered). Only three different disturbance agents were considered: fire, insect, and human, because of a lack of comprehensive datasets for other disturbance types. Fire occurrence and severity came from the USGS Monitoring Trends in Burn Severity (MTBS) dataset, a Landsat-based compilation of annual fire information from 1984 through 2015 for the USA (Eidenshink and others 2007). MTBS data has complete coverage of all public and private lands using consistent methodology (Dennison and Brewer 2014). We used the 30-m resolution burn severity mosaic gridded spatial dataset or rasters to identify burn areas in the GYE (MTBS Project (USDA Forest Service/U.S. Geological Survey)).
Data on insect disturbance came from the U.S. Forest Service's Insect and Disease Aerial Detection Surveys (ADS), which provides a semiquantitative assessment of annual bark beetle outbreak area boundaries from 1970 to 2014 for areas flown in the GYE. The USDA Forest Service Forest Health Protection Program and its partners created and provided the ADS data maps from visual assessment from aircrafts. About 80% of the GYE was flown as a part of the ADS. ADS polygons were combined from USFS Regions 1, 2, and 4 and transformed into a 30-m resolution raster, re-projected, and masked to forested areas.
History of human-caused forest disturbance was derived from the North American Forest Dynamics (NAFD) products. We used the NAFD-NEX timeintegrated map from Landsat images for the conterminous USA from 1986 to 2010 (Goward and others 2015). At 30-m resolution, the map classifies pixels as water, non-forest, ''undisturbed'' forest, outside of study area, or ''disturbance'' with associated year. The Vegetation Change Tracker forest change analysis algorithm analyzes each Landsat image to create spectral indices and track the spectral trajectory for each pixel to produce a forest disturbance map (Huang and others 2010). Human-caused forest disturbance was considered as Disentangling Climate and Disturbance Effects the residual of NAFD data filtered by the MTBS fire data and the ADS bark beetle disturbance data.
To reconcile the spatial resolution differences between our disturbance data (30 m) and other datasets (1 km), we performed a percent disturbance threshold sensitivity analysis (details in Appendix C), and 30% threshold was selected to resample the 30-m resolution disturbance data to 1 km using the ''raster'' package in R (Hijmans 2015). Fire, bark beetle, and human disturbance 1km resolution raster layers were then used to mask NDVI trends. For the empirical modeling, the three disturbance datasets were combined to create three explanatory variables: if disturbed, years since disturbance (since 2014), and percent of 1-km pixel disturbed. The created combined disturbance dataset is shown in Appendix C: Figures C1 and C2.
For meteorological covariates, we used 1-km gridded monthly Daymet V3.0 temperature, precipitation, shortwave radiation, vapor pressure, and snow water equivalent datasets produced by Oak Ridge National Laboratory (Thornton and others 2017). The more recent TopoWx v 1.2.0 (''Topography Weather'') 800-m gridded temperature dataset was also used to compare results with Daymet (Oyler and others 2014). Monthly data were used to calculate mean annual, mean seasonal, delta annual, and delta seasonal variables using CDO functions (CDO 2015). Daily minimum (T min ) and maximum temperature (T max ) data were used to calculate annual aggregated growing degree days (AGDD) as: where base temperature (T base ) equals 5°C and D is the number of days in the year. VPD was calculated using monthly mean/min temperature and vapor pressure data following the American Society of Civil Engineers standardized equations (Walter and others 2005). The re-analyzed Kö ppen-Geiger climate classifications which are based on vegetative types, precipitation, and air temperature to produce a world map with 31 climate classes at 5-arc-min resolution, were also considered as covariates (Kottek and others 2006). We did not account for increased carbon dioxide in this study and suggest this is an important driver of plant productivity to test using process-based modeling. Land cover-type data layers were used to filter NDVI trends to limit spatial analysis to only natural vegetation. The National Land Cover Database 2001 (NLCD 2001 categorizes land cover for the USA into 16 categories based on a decision tree classification of 2001 Landsat satellite data (Homer and others 2007). The NLCD 2001 data were resampled from 30-m to 1-km resolution using the Resample process with the Majority technique in ArcGIS 10.3.1 (ESRI 2016). For filtering NDVI trends, the data were then reclassified into two classes including ''Natural Vegetation'' and an ''Other'' using the ArcGIS Spatial Analysis tool (Appendix B: Table B1).

NDVI Trend Analysis
NDVI trend slope estimates were calculated using a linear least squares regression for each pixel within the GYE using Climate Data Operators (CDO 2015). Using annual aggregated time series has been shown to perform better than removing seasonal cycles in a time series to overcome the adverse effects of inter-annual variability in NDVI time series data on trend estimation performance and produces Tau values that range -1 to 1 (Forkel and others 2013). Thereby, NDVI trend slope estimates were calculated from the annual maximum time series and also calculated using integrated NDVI (Appendix A). The two-sided Mann-Kendall test was used to evaluate the significance of monotonic trends in NDVI (Mann 1945). The Mann-Kendall trend test corrected for autocorrelation is commonly used for analyzing meteorological time series as an alternative to linear regression tests because Mann-Kendall does not require a normal distribution of the data and accounts for temporal autocorrelation (Hamed and Ramachandra Rao 1998;Henebry 2004, 2010). We conducted Mann-Kendall tests using the ''gimms'' R package Tucker 2014, 2016) and used a p < 0.05 significance level threshold.
To isolate the influence of human land use and disturbance, individual pixel NDVI trends were summarized for the GYE before and after filtering out pixels that had land cover layers that were not natural vegetation (using NLCD data), burned since 1984 (using MTBS data), areas affected by bark beetles (using ADS data), and logged stands (using NAFD data) using functions from the ''raster'' R package (Hijmans 2015). NDVI trends were compared across different vegetation types with pixels classified based on NLCD categories using an AN-OVA with a post-hoc Tukey's Honest Significant Difference method using the ''stats'' R package (R Core Team 2017).

Empirical Modeling
Linear regression models were fit using a suite of environmental variables to explain maximum annual NDVI trend slopes. We created a training dataset from a random sample of pixels within the GYE comprising 50% of the total, resulting in a sample size of 61,162 pixels. The remaining 50% of the data were reserved to evaluate how effectively the selected model explains NDVI. We used linear regression to model maximum annual NDVI trend (hereafter NDVI trend) as a function of the meteorological, topographic, and disturbance history explanatory variables listed in Appendix D: Table D1. Each explanatory variable was standardized by subtracting their mean from each value and then dividing by their standard deviation. Due to the short growing season within the GYE and the susceptibility of volcanic soils to summer drying, summer abiotic factors heavily influence plant growth and were thereby the focus of presented analysis, with summer considered as the months of June, July, and August (JJA). To better represent disturbance, it was included as two continuous variables: percent of pixel disturbed (Perc Dist) and years since disturbance since 2014 (YSD), and one binary variable: if observed disturbance (> 0%, If Dist) during the time period 1984-2014. Both continuous variables were included in the models paired with the if disturbed binary, so for example the years since disturbance component was only used to estimate the response for pixels that were disturbed during the time period analyzed.
To check for linear and nonlinear multicollinearity issues among potential covariate variables, a matrix of calculated correlation estimates and scatter plots were plotted using the ''pairs. panels'' function from the ''psych'' R package (Revelle 2016). Covariates were then reduced to those that were less strongly correlated (< 0.75 of the correlation coefficient) and after exploratory analyses of the different explanatory variables. All seasonal and annual meteorological variables were considered, and the variable with the highest explanatory power was selected for use in the analysis (for example, mean summer precipitation for precipitation). Individual linear regression models were fit, coefficient estimates plotted, and sums of squares F-tests calculated to test the importance of the different explanatory variables using the ''stats'' R package (R Core Team 2017). In addition, models of combinations of explanatory variables were fit and compared to determine the model that explained the highest percentage of the variance in NDVI slope based on adjusted Rsquared values. A semivariogram was created from the residuals of the fitted selected model to check for remaining spatial autocorrelation. AIC values were compared for all of the possible model combinations using the variables in the selected model and to the mean-only model (using dredge() from ''MuMIn'' R package; Bartoń 2016).

Greening and Browning in the GYE
A mean greening trend of + 0.0033 NDVI/yr (p < 0.05) was detected in annual maximum NDVI using the USGS EROS NDVI dataset from 1989 to 2014, for 26.5% of pixels within the GYE (Table 1). A mean browning trend of -0.0035 NDVI/yr (p < 0.05) was detected in annual maximum NDVI for 6.2% of GYE pixels (Table 1). The variance in annual maximum NDVI was 0.015. Analysis of growing season integrated NDVI showed similar trends (Appendix A: Table A1). Annual maximum NDVI was used for further analysis since all values were observations, not fitted values. Removing pixels with land cover other than natural vegetation reduced the detected greening trend spatial extent by 32%, to 17.9% of the GYE, and the detected browning spatial extent by 13%, to 5.4% of the GYE (Table 1). The highest concentration of pixels with statistically significant increasing trends in NDVI (greening trends) was in central GYE within Yellowstone National Park ( Figure 3A) and to the west in the Caribou-Targhee National Forest. Decreasing trends in NDVI (browning trends) were detected sporadically on the eastern half of the GYE predominately within Shoshone National Forest, with a few localized areas in the southeast within Bridger-Teton National Forest and Shoshone National Forest ( Figure 3A).
The response in NDVI within different vegetation types varied widely, with areas of greening and browning within each vegetation type (Figure 4). Comparison of the mean NDVI trends between vegetation types indicated unique responses for all categories except for mixed forest, which categorized only a few pixels in the GYE based on NLCD. The shrub/scrub, evergreen forest, and grassland/ herbaceous land cover classes describe over 90% of the GYE. Evergreen forest had a mean greening response (+ 0.001 SD 0.002), whereas shrub/scrub areas showed a weak mean greening response (+ 0.0002 SD 0.002) and grassland/herbaceous areas showed a weak mean browning response (-0.0004 SD 0.002).

Spatial Patterns of NDVI Trends and Disturbance
Filtering the NDVI trends in the natural vegetation raster for areas that were disturbed removed 38% of the greening pixels and 33% of the browning pixels (Table 1). Of the disturbance variables considered, removing pixels within burned areas caused the largest reduction in both the detected greening annual trend in maximum NDVI Impact of removing pixels affected by fire, insect, and human disturbances on detected mean trends in maximum annual NDVI. Values given for annual mean trend, decadal trend, standard deviation (SD), and portion of the GYE (%). Annual maximum NDVI variance was 0.015.
(+ 0.0029/yr) and the fraction of pixels considered (reduced by 3.9 to 14.0% of the GYE) (Table 1). Removing human and insect disturbed areas reduced the fraction of greening pixels by another 1.5 and 1.4%, respectively, without affecting the magnitude of the trend. Removing burned pixels reduced greening area extent by 65% within the National Park boundaries  and removed the strong browning areas to the east within the National Forest ( Figure 3B). Removing human disturbed areas reduced greening areas west of the parks and removing insect disturbed areas reduced browning areas to the east. After filtering out all disturbances, most of the remaining of the greening area (75%) was within National Forest designations, with 30% within Caribou-Targhee National Forest ( Figure 3B).

The Role of Disturbance and Climate in Explaining NDVI Trends
Coefficient estimates for the separate models fits of climate, growing season length, and disturbance variables with NDVI trends aided in identifying key drivers of change in NDVI ( Figure 5). The main effects of change in temperature and mean temperature for both mean annual (MA) and summer (JJA, for the months of summer: June, July, and August) variables showed negative relationships with NDVI trends, while their interaction effect was positive ( Figure 5). Aggregated growing degree days showed a negative relationship with NDVI trends. For the precipitation models, the main effects of change in precipitation and mean precipitation for both annual and summer variables showed positive relationships with NDVI trends, whereas their interaction effect was negative. VPD for both MA and JJA showed positive relationships with NDVI trends, while change in VPD and the interaction coefficients were negative. The disturbance model showed a positive relationship between disturbance terms and NDVI trends. The largest coefficient estimates were for change in summer precipitation (Delta JJA Precip) and disturbance (YSD : Perc Dist). Plotted density curves of covariate distributions for all GYE pixels and greening cells show associations with NDVI trends (Figure 6). The average greening pixel had a higher mean summer precipitation compared to all the GYE ( Figure 6A). Although most of the GYE had decreasing mean summer precipitation, greening pixels experience less drying ( Figure 6B). Greening pixels were limited in their temperature range, with fewer pixels in the hotter temperature range compared to all of the GYE ( Figure 6C). Greening pixels distribution also had a higher kurtosis in change in mean annual temperature ( Figure 6D).
The most parsimonious explanatory model explained 29% of the variation in NDVI slope (r 2 = 0.29, F = 2505, p < 0.0001). Exploratory statistical analysis was performed on all the considered variables, but issues with independence required selecting between collinear variables such as summer precipitation and VPD or mean annual temperature and AGDD. The explanatory model included mean summer precipitation (JJA Precip), change in mean summer precipitation (Delta JJA Precip), mean annual temperature (MA Temp), change in mean annual temperature (Delta MA Temp), years since disturbance (YSD), percent of  Table 2 provides units and data sources for these explanatory variables. The estimated coefficients and 95% CI for each explanatory variable in the model are shown in Figure 7. The second most parsimonious model dropped the interaction term between change in mean summer precipitation and delta mean summer precipitation and had a delta AIC of 652, providing support for including this interaction term (Table 3). Increases in summer precipitation had relatively the strongest positive relationship with NDVI slopes (Figure 7, + 0.0006 (SE 9 9 10 -6 )). Mean summer precipitation also had a positive relationship (+ 0.00035 (SE 1 9 10 -5 )), meaning that wetter sites had a larger greening effect. The interaction between summer precipitation and change in precipitation had a negative relationship (JJA Precip * Delta JJA Precip, -0.00022 (SE 8 9 10 -6 )). Areas in the GYE that were relatively dry had increased NDVI trends with increased summer pre- cipitation compared to wetter areas ( Figure 8B). Areas in the GYE that were relatively wet had increased NDVI trends with decreased summer precipitation while NDVI trends decreased for dry areas experiencing the same reduction in summer precipitation ( Figure 8B).
The second strongest positive relationship was between NDVI slopes and the disturbance threeway interaction (YSD * Perc Dist * If Dist, + 0.00047 (SE 9 9 10 -6 )), with disturbed pixels associated with higher NDVI slopes or greening trends. Disturbed cells that had longer to recover (higher YSD) and were more disturbed (higher Perc Dist) were associated with higher NDVI slopes ( Figure 8A).
Higher mean annual temperature and higher change in mean annual temperature had a negative association with NDVI slopes. The interaction be- Meteorological data are Daymet V3.0 (Thornton and others 2017). Disturbance history data are a combination of three disturbance datasets: North American Forest Dynamics (Goward and others 2015), Monitoring Trends in Burn Severity (Eidenshink and others 2007), and Arial Detection Surveys (Goward and others 2015). JJA June, July, and August for the months of summer. MA mean annual.  tween mean annual temperature and change in mean annual temperature had a positive relationship with NDVI slopes ( Figure 8C, + 0.00037 (SE 8 9 10 -6 )). Moderate warming (0-2°C) for all areas was associated with greening trends. Areas experiencing a greater increase in annual temperature (> 2°C) had a browning response in NDVI. Cooler areas were more sensitive to changes in annual temperature than warmer areas of the GYE, with the coolest areas experiencing the most extreme increases in temperature and were associated with lower NDVI slopes than areas with less warming ( Figure 8C).

Greening and Browning in the GYE
Significant greening and browning trends were limited in their spatial extent to only 17.9 and 5.4% of naturally vegetation areas in the GYE, respectively, with the majority of pixels within the region not showing strong trends during the studied time period of 1989-2014. The order of magnitude difference between the regional average NDVI trend (+ 0.0004 NDVI/yr) and the greening/browning localized trends (+ 0.0033 NDVI/yr, -0.0035 NDVI/yr, respectively) highlights the importance of analyzing plant productivity at finer spatial resolutions that otherwise are aggregated in coarse scale analyses (for example, Zhu and others 2016). Spatial analysis of NDVI slopes at 1-km resolution also revealed the heterogeneity of plant productivity trends across the GYE. Strong greening trends tended to be in the more moist central and western regions of the GYE within Caribou-Targhee National Forest and Yellowstone National Park, whereas strong browning trends tended to be in the more arid eastern regions of the GYE within Shoshone National Forest and Bridger-Teton National Forest. The fine-grained spatial pattern of productivity is consistent with previous research in Yellowstone National Park others 2004, 2016). The average NDVI trend across the entire GYE (+ 0.0004 NDVI/yr) is of the same sign, and general magnitude of NDVI trends detected in boreal climate regions from global studies (+ 0.00038 NDVI/ yr) (Rafique and others 2016). The magnitude of the detected greening trends, or greening areas in the GYE (+ 0.0033 NDVI/yr, p < 0.05) was comparable to trends detected in forests of Canada and Alaska (+ 0.003 NDVI/yr) (Ju and Masek 2016). For detected browning areas in the GYE, the mean trend of -0.0035 NDVI/yr, p < 0.05 was also comparable to detected browning trend values in forests of Canada and Alaska (-0.002 to -0.004 NDVI/yr) (Ju and Masek 2016).
Comparison of NDVI trends between different vegetation types based on NLCD categories indicated different responses by vegetation type, although the range in responses within each category varied widely. Evergreen forests showed a mean greening response and comprised the areas of strongest detected greening trends in the GYE. In contrast, more of the browning areas were categorized as shrub/scrub or grassland/herbaceous. The broad NLCD categories limit inference can be made from this analysis. For example, ''evergreen forests'' includes the low, mid, and high elevation forests of the GYE that have unique vegetation compositions and potentially unique responses in NDVI. Although comparison of NDVI responses at the species level is beyond the scope of this paper, it remains an important area of research for illuminating the underlying mechanism of changes in plant productivity.

The Role of Disturbance in Explaining NDVI Trends
Filtering the NDVI trends in the natural vegetation for areas that were disturbed removed 38% of the greening pixels and 33% of the browning pixels, providing evidence for the importance of disturbance in explaining detected changes in plant productivity. Filtering for fire disturbance alone removed 22% of the areas with greening trends, with human disturbance removing 8% and bark beetle infestations removing 8% of the greening trends. Statistical modeling of disturbance effects on NDVI trends showed that disturbed pixels that had longer to recover (higher YSD) and were more disturbed (higher Perc Dist) were associated with higher NDVI slopes or greening trends. Recently disturbed pixels were associated with negative NDVI slopes or browning trends. These results are consistent with studies of rapid post-disturbance recovery in the GYE others 2004, 2016;others 2005, 2013;Zhao and others 2016) and changes in productivity in northern latitude forests in Alaska and Canada (Ju and Masek 2016;Sulla-Menashe and others 2018) and reinforce the need for considering landscape history in analysis of plant productivity trends. As previously mentioned, recovery of NDVI values post-disturbance is sensitive to herbaceous plant recovery, including potential opportunistic exotic species (for example, Canada thistle (Cirsium arvense) and prickly lettuce (Lactuca serriola) (Turner and others 2008)). Thereby, positive NDVI trends do not indicate whether or not the forest canopy has regenerated specifically nor indicate the composition of native vs. exotic species without further analysis.

The Role of Climate in Explaining NDVI Trends
While statistical modeling supported the importance of disturbance in explaining greening and browning trends in the GYE, it also provided some evidence for the influence of meteorological variables on NDVI trends. Summer precipitation and annual temperature emerged as the meteorological factors with the strongest influence on NDVI trends. The importance of the interaction between mean summer (JJA) precipitation and change in JJA precipitation demonstrates the importance of the bioclimatic context of climate change. Increases in summer precipitation benefitted areas that were relatively dry. It follows that vegetation in the wetter areas were not water limited and thereby did not benefit from increased summer precipitation. The linear model fit instead suggests reductions in summer precipitation benefitted wetter areas, yet precipitation type (snow vs. rain) or nonlinear effects not represented might explain this result. Previous studies suggest extension of the growing season length and water stress, represented here by AGDD and VPD, respectively, drive plant productivity in the boreal climate region (Angert and others 2005;Boisvenue and Running 2006;D'Arrigo and others 2008;Williams and others 2011). Although both AGDD and VPD showed statistically significant relationships with NDVI trends, temperature and precipitation explained more of the observed trends than these variables.
Increased temperatures in boreal climate regions can facilitate plant growth through an elongated growing season without precipitation limitation (Boisvenue and Running 2006), and this is supported by the location of detected greening trends primarily in forests in the boreal climate of western and central GYE. Indeed, process-based model simulations also indicate that warming temperatures will favor positive growth responses for lodgepole pine dominated forests in Yellowstone National Park (Smithwick and others 2009). Yet, the coolest areas in the GYE, the Absaroka Mountain Range and the Wind River Range are also areas that experienced the greatest increase in mean annual temperatures and high-severity mortality due to mountain pine beetle (MacFarlane and others 2013) and were associated with browning trends. Climate warming benefits mountain pine beetles by increasing overwinter survival leading to increased outbreak populations (Bentz and others 1991). MacFarlane and others (2013) used an alternative aerial survey method and found whitebark pine mortality caused by mountain pine beetle to be more widespread than as reported in the Aerial Detection Survey data. Thereby, these detections of browning trends may be interpreted as signals of mountain pine beetle mortality that were not captured in the Aerial Detection Survey data.
Recent studies suggest the decoupling of northern latitude forest growth and warming temperatures, with other factors such as fire disturbance, insect and disease outbreaks, and water stress driving detected trends in forest productivity (Goetz and others 2005). Our results provide further evidence of the important effects of disturbance on productivity trends and the importance of precipitation during the growing season for tree growth. This may suggest that the greening trends previously detected and attributed to warming temperatures and lengthening of the growing season are not spatially consistent across boreal climate zones, nor are these relationships expected to continue under future climate projections. Growing evidence already suggests that tree growth in northern latitudes is being limited by drought stress (Angert and others 2005;others 2008, 2009;Williams and others 2011). For undisturbed forests during their study period, Sulla-Menashe and others (2018) found that the direct effects of climate change showed minimal influence on NDVI trends overall, but with varying effects based on local bioclimatic differences. Similarly, we show that the effects of precipitation and temperature changes on forest productivity are highly dependent on the local bioclimatic context, driving the sign and magnitude of detected NDVI trends.

Further Considerations
Global studies suggest that multiple factors drive productivity trends in northern latitudes (Nemani and others 2003;Zhu and others 2016). Beyond the climatic and disturbance factors explored in this paper, other possible mechanisms for detected increases in NDVI include a fertilization effect from increases in atmospheric carbon dioxide or nitrogen deposition. Increases in atmospheric carbon dioxide are expected to increase plant growth. Global ecosystem modeling simulations suggest carbon dioxide fertilization effects may explain more of the observed greening trends in satellite data than climate (Piao and others 2006;Zhu and others 2016). Although accounting for carbon dioxide fertilization effects was outside the scope of this paper, it remains a potentially important driver of detected NDVI trends. Process-based ecosystem modeling offers an opportunity for testing and quantifying the relative effects of carbon dioxide increases and climate change on plant productivity in the GYE. Increases in nitrogen deposition can also stimulate plant growth. Soils vary within the GYE from nutrient-poor rhyolitic soils of the Yellowstone Caldera to andesitic soils of higher-elevation forests that are higher in nutrients and water holding capacity (Despain 1990). The productivity of lodgepole pine forests on rhyolitic soils is generally considered nitrogen limited (Fahey and others 1985); thereby, increased nitrogen availability could act as a fertilizer. Furthermore, increased annual temperatures and precipitation could stimulate plant growth indirectly through increased soil temperature and moisture facilitating nutrient and water uptake (Reddell and others 1985).
Finally, it is worth noting that the inferences in this study are limited by the temporal extent and quality of the disturbance data. While the NDVI time series was from 1989 to 2014, the disturbance data did not represent even this discrete time frame. For example, extensive logging occurred in the Caribou-Targhee National Forest between 1950 and 1990 (Parmenter and others 2003), with only logging starting in 1986 and beyond captured in the NAFD disturbance data. The inclusion of U.S. Forest Service forest stand maps could aid in controlling for disturbance more fully in analyses of changes in productivity (Tinker and others 2003). In fact, century-long time scales are required to understand the regrowth legacies of disturbances Disentangling Climate and Disturbance Effects (Williams and others 2016), with carbon stocks and fluxes stabilizing after 100 years in subalpine Rocky Mountain forests (Bradford and others 2008).
Larger and more synchronous fires have occurred during hot and dry periods, suggesting climate is the dominant driver of fire behavior in the region (Balling and others 1992;Morgan and others 2008;Marlon and others 2012). Predicted climate induced changes in the fire return interval in the GYE could drastically alter forest regeneration, from forest extent, species composition, to age-class distribution (Romme and Turner 1991;Westerling and others 2011). Climate warming is also predicted to increase mountain pine beetle outbreaks causing forest die-off (Logan and Powell 2001;Logan and others 2010). With disturbance playing such a significant role in forest dynamics in the GYE and other northern latitude temperate and boreal forests, further research is needed on the interactive effects of temperature, precipitation, and disturbance on forest productivity and how these relationships may change under future climate conditions.

AC KNOWLEDGEMENT S
We thank Steve Cherry, Dave Roberts, and Mark Greenwood for discussions on statistical approaches. Thank you to Christopher Williams, Huan Gu, and Yu Zhou for sharing processed ADS and NAFD data and Arjun Adhikari for sharing NLCD. The article was improved by helpful suggestions from Robert Keane, Monica Turner, and Andrew Hansen. The research was supported, in part, by NSF Grant GSS-1461590; Grant No. G15AP00073 from the United States Geological Survey; NSF Graduate Research Fellowship Award No. DGE-1632134; and the Montana Institute on Ecosystems with support from NSF-IIA-1443108 and EPS-1101342.

OP E N A C C ES S
This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.