Fine-Scale NDVI Reconstruction Back to 1906 from Tree-Rings in the Greater Yellowstone Ecosystem

: Global warming and related disturbances, such as drought, water, and heat stress, are causing forest decline resulting in regime shifts. Conventional studies have combined tree-ring width (TRW) and the normalized difference vegetation index (NDVI) to reconstruct NDVI values and ignored the influences of mixed land covers. We built an integrated TRW-NDVI model and reconstructed the annual NDVI maps by using 622 Landsat satellite images and tree cores from 15 plots using point-by-point regression. Our model performed well in the study area, as demonstrated by significant reconstructions for 71.14% ( p < 0.05) of the area with the exclusion of water and barren areas. The error rate between the reconstructed NDVI using the conventional approach and our approach could reach 10.36%. The 30 m resolution reconstructed NDVI images in the recent 100 years clearly displayed a decrease in vegetation density and detected decades-long regime shifts from 1906 to 2015. Our study site experienced five regime shifts, markedly the 1930s and 1950s, which were megadroughts across North America. With fine resolution maps, regime shifts could be observed annually at the centennial scale. They can also be used to understand how the Yellowstone ecosystem has gradually changed with its ecological legacies in the last century.


Introduction
Forests cover 30% of the world's land surface and support ecosystem functioning and human societies [1]. Vegetation plays a vital role in stabilizing the ecosystem, yet its interaction with climate, human activities, and related disturbances trigger ecosystems to drift from the standard thresholds resulting in regime shifts [2]. The ecological regime shifts definition, according to Biggs et al. [3], emphasized three characteristics: large, abrupt, and long-lasting changes in the ecosystems, which often have considerable impacts on human economies and societies. "Abrupt" refers to the instant variation possibly caused by short-term disturbances, such as fire, insect outbreak, or drought. These ecological disturbances break down the balances in the previous ecological system. "Large" and "long-lasting" refer to persistent observation.
Currently, the temporal observation of forest change on a broad scale comes from satellite imagery, such as Landsat series data [4,5], MODIS data [6,7], and SPOT data [8,9]. The earliest Landsat series satellite was launched in the 1970s [10]. Before the 1970s, there were no consistent high-resolution vegetation density data which influenced our longterm observations to detect regime shifts in the forest. Tree-ring width (TRW) is a sensitive index representing the growth of the xylem as a response to canopy structure [11]. Annual tree-ring chronologies are a widely used proxy that allows us to reconstruct past climate, ranging from annual to centennial time scales [12,13]. The NDVI is a measure of the density of the vegetation based on the reflection of red light and infrared light [14]. The rela-tionship between NDVI and TRW is based on the assumption that the total leaf area represented by NDVI can determine the light interception and CO2 assimilation that goes into tree-ring widths. With enough nutrients, the xylem in tree rings can grow large during high NDVI years [15].
Currently, the NDVI reconstruction with tree-ring data is still in the phase of pointby-point regression; however, the spatial resolution is coarse [16][17][18]. The reconstruction can reflect the regime shift on a broad scale and the global drivers causing the shift. However, there are still some subtle forest regime shifts at a finer scale and some regional drivers are ignored by coarse resolution images. One coarse grid can contain different land cover types, such as forest, grassland, barren land, and water, while tree-ring data will just reflect the growth of the forests. There could be a mismatch between the NDVI and tree-ring data which could lead to some inevitable errors in the reconstruction. Brehaut and Danby [19] found an inconsistency between NDVI and TRW in Southwest Yukon and Canada but the reason was unclear. Correa-Díaz et al. [20] reconstructed correlation maps between tree-ring width and winter NDVI at the 250 m resolution, which displayed more landscape detail than the previous studies. We anticipated that more landscape structure and vegetation variation could be shown if we could reconstruct NDVI using finer resolution images (30 m). Liu et al. [21] reconstructed cumulative NDVI with a single linear regression model, and they found that the climate factors influencing NDVI, such as precipitation and temperature, were also driving tree-ring width. To have a more robust model, it is also important to consider slope, elevation, and lag year effects when modeling. Thus, monitoring the forest regime shift on a finer scale and considering both topographic and temporal factors allows us to build more accurate models at a finer resolution. In this study, we not only reconstruct the NDVI maps but also explore the factors influencing the correlation between NDVI and TRW.
Based on our current understanding, we aimed to assess forest regime shifts on a longer time scale and at a finer spatial resolution than previous studies. The specific objectives of the study were to i) make a 100-year series of NDVI maps with tree-ring width to observe regime shifts; ii) explore the factors influencing the correlation between NDVI and TRW. We envision our results as a series of NDVI maps where we could vividly "see" the development of the forest structure and ecological succession of the Greater Yellowstone Ecosystem in the recent one hundred years.

Site Description
Our study area ( Figure 1) is located on the northeast of the Greater Yellowstone Ecosystem. The average temperature is around 8 °C in summer (June through August) and −7 °C in winter (December through February). The annual rainfall is 510 mm. The plot sites were located on the slopes of mountainous areas in the Greater Yellowstone Ecosystem. The study site is laid out along an elevation gradient in a 12 km × 12 km area. The dominant trees are lodgepole pine (Pinus contorta Douglas) and whitebark pine (Pinus albicaulis Engelm.), whose mortality increased a lot when the temperature reached around −10 °C in the growing season [22,23], and due to the mountain pine beetle (MPB, Dendroctonus ponderosae Hopkins), and white pine blister rust (Cronartium ribicola J.C. Fisch.) outbreaks since 2000. Other common species include Engelmann spruce (Picea engelmannii Parry ex Engelm.), subalpine fir (Abies lasiocarpa (Hook.) Nutt.), and Douglas fir (Pseudotsuga menziesii (Mirbel) Franco). The study area was selected to represent the vegetation types and elevational changes in the Greater Yellowstone Ecosystem, with elevations ranging from 1170 m to 3087 m. Conventional studies have just used one pixel to represent the whole study area, and conventional research will regard the pixel as the forest where trees take up more than 50% of the area. However, instead of a pure forest pixel, the area also includes lakes, rock, and grass. We chose the 7 July 2015 Landsat 8 image from the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/, accessed on 2 August 2020) as our base image to map the land use and land cover distribution. We classified the whole area into four land covers (forest, grassland, barren land, and water) using the random forest method, which is an ensemble learning method to generate a series of decision trees for one pixel and assign it as one specific land cover with highest probabilities [24]. The whole classification had four land covers: forest, grassland, barren land, and water, with 200 samples for training and 100 samples for validation. Figure 1. The study area is located in the mountainous area in the Greater Yellowstone Ecosystem. The green circles represent the final 15 plots that we chose for this study.

Tree-Ring Measurement, NDVI mapping, and Climate Data Collection
Along the southwest aspect of South Bird Mountain, we placed a plot every 250 m along an elevational gradient and collected tree cores in 16 plots (ECO100, ECO200, ECO400, ECO500, ECO600, ECO800, ECO900, ECO1000, ECO1100, ECO1300, SBM200, SBM300, SBM400, SBM500, SBM600, and SBM700). Within each stand, we selected 10 trees closest to the center point whose diameters at breast height were more than 10 cm. We extracted two cores from each tree with an increment borer at 30 cm height. Samples were glued to prefabricated core mounts with water-soluble white glue, sanded with a progressively finer sandpaper, (120, 220, 320, and 400 grit and hand sanded with 30-micron sanding film), and then skeleton plotted for age dating [13].
In addition to the 16 plots, we also used tree-ring data from other studies from the North American Dendroecological Fieldweek on four sites (CBS, BTP, TOW, and RCCF1) [25,26]. We used the skeleton plot method to cross-date the samples with an internally derived master chronology for each species [13]. Once the cores were accurately dated, we used a TA4030H1-S6Velmex measuring system, (Velmex Inc., New York, NY, USA) to measure each core with a 0.001 mm precision, and measurement and dating quality were checked with the computer program COFECHA version 6.06P, (Laboratory of Tree-Ring Research, University of Arizona, Tucson, Arizona, AZ, USA) [27], and segments with low correlations were examined and corrected at the microscope. We detrended and standardized the series with the age-dependent spline in ARSTAN version 48d2_ (Win, Tree-Ring Laboratory, Lamont Doherty Earth Observatory of Columbia University, New York, NY, USA) for each plot to remove the age-related tendency and the growth difference between each tree [28]. ARSTAN also calculated the autocorrelation and the Expressed Population Signal (EPS), and we cut off the species-level chronologies at an EPS of 0.85. It also calculated the rbar, which is a running correlation showing the synchrony of the species-level chronologies in our study. We set up a threshold to choose the plots which considered two aspects: One is to cover as many plots as possible despite the narrow temporal range. The other is to cover as long a temporal range as possible, but it would exclude more tree cores with young plot age. Considering these two thresholds, we preferred to have old-aged trees over more numbers of plots. The tree-ring series were subset to span the time period from 1905 to 2015 thus, we excluded ECO100, ECO200, SBM200, SBM400, and SBM500 plots with young-aged trees. Finally, 15 plots were selected out of a total of 20 plots (Figure 1). The autocorrelation showed that most plots had a first-order model, which indicated the previous year's climate may influence the tree growth in the current year.
We downloaded the Landsat series (Landsat 3-Landsat 8) images (622 images) from March to October from 1980 to 2015 from the USGS (https://earthexplorer.usgs.gov/, accessed on 2 August 2020). The images have a fine spatial resolution (30 m), enabling the identification of land cover on the landscape. We calculated NDVI using Equation (1). NDVI represents the density and health of vegetation and is a good indicator of the local climate, such as drought or precipitation [29,30]. Each grid (30 m × 30 m) in each year could have more than one NDVI value on different dates. We took the median of them to get annual NDVI for the growing season. We also collected the annual temperature, precipitation, and Palmer Drought Severity Index (PDSI) at https://wrcc.dri.edu/wwdt/time/spanning from 1905 to 2015 (accessed on 2 August 2020).
where the InR and Red are the infrared and red bands respectively from the Landsat sensors.

Model Building
We developed a NDVI model derived from the forests and applied the model to each grid cell in the study area. We split the period of our study into two sub-periods of the calibration sub-period  and reconstruction sub-period . TRW and NDVI both reflect vegetation growth and the local climate. The NDVI values in each grid could be correlated with plots with similar slopes and elevations. Conventional studies used latitude and longitude to define the distances between the grid and selected plots and ignored elevation effects. Our study area is small but has a large elevation gradient. We defined the distance using elevation and slope with the Digital Elevation Model (DEM) data from the USGS. To unify the various scales for elevation and slope, we standardized the units of elevation and slope with Min-Max scaling, which rescaled the two features (slope and elevation) into the same range [31,32], calculated the distances between each grid using Equation (2), and chose the five "closest" chronologies for each grid as potential candidates. The five potential candidate plots had similar elevation and slope with the grid, and their vegetation distributions were also similar.
where d is the comprehensive distance considering the relevant elevation and relevant slope, E and S are the rescaled elevation and slope for each grid, E_plot and S_plot are the rescaled elevation and slope for each plot.
To make full use of those five chronologies, we reduced the redundancy and noise using principal component analysis (PCA) to extract the first component (FC), the second component (SC), and the third component (TC). Previous studies and the autocorrelation in our research suggest that reconstructing NDVI should also take into account a time lag (at least one year) [33,34]. In our study area, there are multiple types of land cover, so we built three models with increasing complexity to match the NDVI to different land covers (Equations (3)- (5)). The three sub-models were applied to each pixel. We used the AIC criteria to judge which sub-model was the best among those three and calculated the pvalue using the best sub-model.
NDVI = a1 × FC + a2 × FC_pr + a3 × SC + a4 × SC_pr + b (4) NDVI = a1 × FC + a2 × FC_pr + a3 × SC + a4 × SC_pr+ a5 × TC + a6 × TC_pr + b (5) where NDVIs are the simulated NDVI value; FC, SC, and TC are the first component, the second component, and the third component from the PCA in the current year, FC_pr, SC_pr, and TC_pr are those components in the previous year. a1, a2, a3, a4, a5, a6, and b are the parameters for each variable and constant values. After classifying the land cover and conducting linear regression for each pixel, we determined the NDVI values in each grid cell from 1906 to 1979. Among cells in the study area, we averaged annual NDVIs for each type of land cover (forest, grass, rock, and water) and compared the overall NDVIs with them. Time series analysis for the NDVI from 1906 through to 2015 was also performed using the strucchange package in R [35], the most widely used method to identify regime shifts. The program identifies breakpoints (step change in the mean) in the NDVI time series via an estimation based on the minimization of the residual sum of squares (RSS) within the series and the Bayesian Information Criterion (BIC) score for scoring the optimum number of breakpoints [35] along with the calculation of confidence intervals based on equations from [36]. We then explored the relationships between NDVI and local climate using precipitation, temperature, and PDSI in our model.

Evaluating the Model Performance
After computing regression for each grid cell (16,000,000 cells in total), we could observe the annual summer NDVI variation of the whole area and each land cover. Then we conducted three evaluations on the model and the simulated values: (1) We evaluated the model performance on the whole area and each land cover. We calculated p-values for each cell to judge whether the p-value meets a significant level (p < 0.05, 0.05 < p < 0.1, and p > 0.1). This statistical analysis was conducted for each land cover. (2) We evaluated the model performance across the whole study time using leave-oneout cross-validation. We chose one specific year as an independent dataset and involved all years' datum, excluding the specific year in the analysis. We calculated the p-value for the grid without that year's data and judged whether the p-value met the significant level. Our study time covers 36 years (1980-2015), so we built 36 separated datasets to analyze the influence of the excluded year on the whole regression. If we exclude one year and the areas with low p-values (p < 0.05) increase, this indicates the data quality in that year is poor, and exclusion may improve the whole regression. If the areas with low p-values (p < 0.05) decrease, this indicates the data quality in that year is good, and exclusion may weaken the whole regression. The ideal situation is that for each exclusion, the size of low P areas is stable, which means the model performance will be spoiled with low p-values or not too dependent on some years. (3) We evaluated the gap between NDVI in the whole area and NDVI in the forest areas.
The TRW reflected the growth of the trees or forests but not the water or barren land, so it should match NDVI values in the forest. The previous studies [21,37] used coarse resolution images where all land cover was mixed and directly built the regression between the NDVI in the whole area and the TRW in the tree plots. We introduced the error rate to calculate the gap between those two NDVI values (NDVI for forest and NDVI for whole areas) (Equation (6)) and evaluated whether the gap will influence further analysis, such as the r-value, with other climate factors.
Error rate = |Approximate value − Exact value| |Exact value| × 100% where approximate values include NDVI values for the whole area, r-values for the whole area, and the exact value include NDVI values for forest.

Reconstructed NDVI Maps and Land Cover Map
With the tree-ring series from 1905 to 2015 and the NDVI images from 1980 to 2015, the regression could simulate the NDVI images from 1906 to 1979 (without imagery from 1904, the time lag regression could not simulate the NDVI image for 1905). The annual NDVI maps in the growing season are available in Supplementary Material. From the 1906 simulated image (Figure 2A), we found that most areas still had comparatively high vegetation density whose overall NDVI reached 0.32 ( Figure 2E).  1906-1921, 1922-1937, 1938-1953, 1954-1983, and 1984-2015 (Figures 2E and 3) with the four break dates (1921, 1937, 1953, and 1983) as indicated by the lowest BIC for the NDVI. The average NDVI decreased from 1921 to the mid-1930s, then increased from 1940 through 1950, then declined and reached a mean level equal to 0.2628 until early 1980, and it again increased and reached a mean level equal to 0.3113 through 2015.  As we previously assumed, most areas in our study field were forests. The conventional classification would interpret those two pixels as the pure forest in the coarse-scale map and ignored other minor land covers. The whole classification had a satisfying overall accuracy (95%). The accuracies of the water and rock areas are good, but that of the forest is a little lower (80%), which contributed to the mixture of the tree, brush, and grass (Table 1).

Spatial Effect to the Model Selection
NDVI maps from 1906 to 1979 were reconstructed with the integrated model, and the performance of the simulation needs to be evaluated. In the p-value maps including all land cover types (Figure 5a), we classified all p-values into three levels (p < 0.05; 0.05< p < 0.1; p > 0.1). We excluded barren land and water areas and generated the exclusion map ( Figure 5b). The p-value in the first level (71.14% for exclusion map) was less than 0.05, which met the conventional statistical requirement. The p-value in the second level (11.31% for the exclusion map) was between 0.05 and 0.1. If we excluded barren land and water, the ratio whose p-value (< 0.05) increased by 17% to 71.14%.
The three-component sub-model fitted 92.81% of the forest and 53.20% of the grassland even though AIC had severe penalties to the redundant parameters ( Figure 6 and Table 2). About 75.37% of the water and 66.92% of the barren land were described by the first component sub-model reflecting their simple characteristics.  Figure 6. The sub-model preference for each land cover.
From the validation of the models, if we did not exclude the barren land and water, 53.97% of the areas met the 0.05 standards. If we excluded the barren land and water, nearly 71% of the area met the 0.05 alpha level. We conducted a leave-one-out validation to verify that the good fit was not due to one good year increasing the general goodness of the regression fit (Table 3). If we excluded 2012 data, the ratio of good match areas resolved to be 76.97% because, in 2012, the only available Landsat series satellite was Landsat 7, but there were some data gaps in the images due to the failure of the Scan Line Corrector (SLC). However, on the whole, the fitness of each year was very even, which verifies that the model works well for the entire time series.

The Similarities and the Difference between the Whole Area and Forest
In most cases, overall, NDVI shares a similar pattern with the forest NDVI. In the precipitation analysis, overall NDVI and forest NDVI both reached significant positive relationships in June (Overall_Jun = 0.25, Forest_Jun = 0.25; Figure 7A). In temperature analysis, overall NDVI and forest NDVI both showed a significant negative correlation with temperature in June (Overall_Jun = −0.20, Forest_Jun = −0.21; Figure 7B). In the PDSI analysis, overall NDVI and forest NDVI both displayed significant positive correlations with PDSI in July, August, and September (Overall_Jul = 0.20, Forest_Jul = 0.20; Over-all_Aug = 0.22, Forest_Aug = 0.21; Overall_Sep = 0.23, Forest_Sep = 0.22; Figure 7C). The absolute difference between overall NDVI and forest NDVI was not large, and the absolute gap between overall correlation and forest correlation was also small. However, the range of NDVI (−1, 1) and r (0, 1) was small, so we introduced the relative error rate. The NDVI error rate between overall NDVI and forest NDVI was 10.36%. If we bring those NDVI values with high error rates into a correlation analysis with other climate factors (precipitation, temperature, and PDSI), it will cause high annual error rates for r values which average the correlations from January to December. The annual error rates for precipitation, temperature, and PDSI were 20.25%, 66.92%, and 11.00%. In some extreme months, such as precipitation in April and temperature in November, the error rates for correlations were 71.89% and 462.18%. Even though r values in these months did not meet significant levels, the error rates still indicated the big gap between forest NDVI and overall NDVI and spoiled the quality of our whole dataset. Water may play an important role in generating these big gaps because its correlation was usually opposite to those in the forest, such as precipitation and temperature in June (Water_preci = −0.10, Forest _Aug = 0.25; Water_temp = 0.23, Forest_temp = −0.21; Figure 7A,B).

Regime Shift in NDVI from 1906 to 2015
By analyzing century-long NDVI images from 1906 to 2015, our study site experienced five periods of regime shifts in vegetation density during some specific dates : 1921, 1937, 1953, and 1983 (the blue dashed lines in Figure 2E). The assumption that the vegetation density is decreasing (the red dashed line in Figure 2E) in the short term cannot tell the whole story because the available satellite images started from the 1980s. However, with our reconstructed NDVI data, we could identify a subtle increasing tendency (the gray solid line in Figure 2E, p < 0.05) over the long term.
The periodic changes in the landscape mosaic, particularly the vegetation composition and structure in the Greater Yellowstone Ecosystem, have been attributed to multiple disturbances occurring through time and across spaces [38]. Disturbances range from finescale, such as geomorphological changes, to broad-scale drought events driven by climate and intensive land-use management interventions in this region and across the western United States through centuries [39,40]. The study site did not experience any severe canopy replacing fire events after 1900 during the period of this study, as shown by a fire history study by Brown et al. [41].
With the motive of maintaining ideal wildlife populations in the park, the park has a history of intervening in conservation and hunting approaches to increase herds or to reduce them depending on the resources available [42]. Due to cold and lack of food availability during the winter of the late 1910s, the park experienced an unusually high die-off of the large grazing and browsing herbivores which recovered during the 1920s following care and winter-feeding programs by the park [42].
In the northern region of the Greater Yellowstone Ecosystem, increased browsing of the stem and root sprouts by a large number of elk and other ungulates and herbivores from 1886 to 1930 was found to affect forest regeneration of aspen [43,44]. The artificial reduction in elk populations was extensive during dry periods after the 1930s Dust Bowl drought. Intentional hunting of an overabundant number of elk and ungulate populations reduced their number to as low as 25% of their original populations by the late 1960s, resulting in the reduction in browsing, which favored regeneration of the vegetation in the area [43][44][45]. In addition, widespread mountain pine beetle outbreaks occurred across the USA and the Greater Yellowstone Ecosystem, heavily affecting whitebark pine stands in the 1930s [46]. The combined effects from a large number of elk, mountain pine beetle outbreaks, and the Dust Bowl drought contributed to the low vegetation density in 1930 (NDVI = 0.09) in our study area. Then the drought passed, and the number of elk declined. This instigated forest recovery during the 1940s and maintained a high vegetation density in the 1950s ( Figure 2E). The average NDVI values for our study areas in 1906, 1930, and 1945 (Figure 2A-C) experienced large NDVI variations. From the 1970s to the 2010s, the vegetation density was stable with some fluctuations. Then in the 2000s, there was a slight decline. Occurrences of the outbreaks of the mountain pine beetle and the white pine blister rust in the study site since 2000 have affected over 100,000 acres of forest dominated by whitebark pine, lodgepole pine, and limber pine in the study site [47].

Evaluating the Combined Model
The integrated model took three measures to fit the data: reducing the redundant signals and signal-to-noise ratio, applying three sub-models on various land cover types, and taking into account the growth lag. We chose the "closest" five plots using their standardized elevation and slope and conducted a Principal Component Analysis on those plots to reduce the data dimension, extract useful information, and avoid overfitting of the regression. Unlike the conventional studies, which just regarded the whole area as a couple of pixels at coarse resolution, we took into account the mixed land cover, analyzed the images at a finer resolution (30 m), and each pixel was fit with three sub-models: one principle component sub-model, two principal components sub-model, and three principal components sub-model. The model with the lowest AIC value was selected as our final model fit. A high percentage of the low p-value (<0.05) areas indicated the combined (including three sub-models) model was robust in the study area and performed better in the vegetated areas. The gap (17%) between the ratios whose p-value was less than 0.05 of including and excluding water and barren land suggested that we should not neglect this spatial effect on the reconstruction of NDVI. This process demonstrated that forest and grasslands were more complex, needing a three-component model, while barren land and water were more basic, only needing a one-component model. We recommend this approach for future modeling attempts as it allows for a better fit depending on the land cover type, similar to what Brehaut and Danby [19] found. Cook et al. [48] put forward the point-by-point regression with the AIC standard to judge the performance for the dendrochronology-climate model. In our study, the model with fewer variables fit the barren land and water pixels with low AIC values, which could indicate that barren land and water are not sensitive to tree-ring width variations. Barren land and water do not reflect the vegetation changes very well. However, a few water pixels tended to fit the threecomponent model because of the inclusion of swamps in this land classification. Similarly, some barren land pixels also tended to choose the three-component model because of the presence of some sparse vegetation. Over 90% of the forest areas fit the three-component models because the NDVI and TRW are both incorporated very well in the forest areas. The coupling relationship between those two variables was very stable, which is supported by previous studies [37,49,50]. The percentage of grassland described by one and three-component models was found to be intermediate between the water and forest pixels. Half of the grassland pixels (53.2%) fit the three components model, and slightly less (45.25%) fit the one-component model because the grassland could present the NDVI variation but could not completely reflect the TRW changes.
Our model also considered the lag effect. According to the autocorrelation in our analysis and the results from Nothdurft et al. [51], the growth lag for pine usually lasts 19 months. Nothdurft et al. [51] also claimed that pine in pure stands and mixed stands could both benefit from precipitation during the spring and summer months of the previous year. The possible explanation for the time lag is the stomatal closure during drought in the previous year. To avoid water loss in the summer, the stomata close, and meanwhile, the CO2 cannot be absorbed by the vegetation. Then the rate of photosynthesis would drop dramatically, and carbohydrate storage would reduce, which influenced the foliage growth the next year. The reconstructed NDVI and precipitation in the previous June showed a significant positive relationship (r = 0.27, p < 0.05; Figure 7A).
The correlation for the forest and the correlation for the whole area in June were different (r_forest = −0.22; r_whole = −0.20; Figure 7B). The NDVI maps with fine resolution generated by our models were affected by the local land cover and local climate. The main reason is the inclusion of water areas whose correlation reached 0.23 even though the water areas took up only 2.38% of the whole area (Table 1). However, water presented a negative relationship with NDVI decreasing the skill of the model. If our study applied a coarser resolution and did not separate the water and forests in some areas where there were more lakes or rivers, the error between the forest and the whole area could be worse. Because of the finer resolution of the imagery that we used, the NDVI derived from the forest was not confounded by the inclusion of water or barren land in the pixel composition.

Limiting Factors for the Coupling with NDVI and TRW
For the NDVI simulation in the Greater Yellowstone Ecosystem, most areas (71.14%) achieved a significant level (p < 0.05). However, the robustness of the model is still under some constraints if we intend to apply the model in other seasons and other areas. Some studies failed to find consistent correlations between NDVI and TRW [19,52,53]. In our study, we found out that the stable relationship between NDVI and TRW only occurred in some scenarios, some of which are discussed below. The coupling between the NDVI and TRW should meet two conditions. The first one is that the forest areas should be under some regular climate stress. The stable correlation for NDVI and TRW was found in our research, which is consistent with Wang et al. [50], whose study areas are in the Western US and also experienced a lack of precipitation. A lot of studies successfully simulated NDVI from TRW in Russia and Alaska [16,[54][55][56], which experienced enduring cold temperatures. With correlations between NDVI and tree rings, Babst et al. [57] detected the outbreak of autumnal moth (Epirrita autumnata) in Northern Sweden, which caused defoliation but did not directly kill the trees. Other patterns of stressors were detected from NDVI using TR, respectively, heavy metal pollution in a heavily industrialized area [58]. With Leibig's law of the minimum, the most limiting factor controls the main stress on forest growth and synchronizes the stress to the leaves (recorded by NDVI) and the xylem (recorded by TRW) [13]. In this case, the growth of the canopy and the cambium have a similar response resulting in a relationship between NDVI and ring width. The most common controlling factors are low temperature or lack of precipitation. When we analyzed the reconstructed NDVI and the PDSI in Figure 7C, the correlations were positive from June to September, which indicated that the drought was the main limiting factor. During drought events, the stomata close, and CO2 cannot enter into the leaf, which reduces the photosynthetic potential for the leaves and, therefore, the carbohydrates for ring growth. The NDVI and TRW both dropped during that time.
Another pre-requisite is that the magnitude and the active period of the responses from the foliation and cambium should be similar. Between September and October, cambium cells still grew, and positive correlations between NDVI derived in the forest and temperature were also found in Figure 7B even though the significant level did not meet 0.05, which indicated that the low temperature was one of the stresses but not the major one. With various sensitivities to drought and low temperatures, the response magnitude to the multiple stresses to the cambium and leaves may be different. The main stress could be the decreasing temperature after October and the frost tolerant xylem cells suffered extracellular freezing, which resulted in the rupturing of the cell [59]. Although the cambium was inactive, there was no defoliation of a great number of needles or yellowing of the needles, but the trees still photosynthesized at a low rate during that time (the magnitude of the response to decreasing was different between leaves and cambium affecting NDVI and TRW differently). The cambium cells and leaf cells showed an asynchronous active period, even though the cold period may cause the same stress for those two types of cells. The lags for the active time can occur in the early spring when the cambium cells are reactivated in springtime, but the new leaves are not matured enough to conduct photosynthesis effectively, and that is when there will not be a good fit between NDVI and TRW. Other factors could also lead to a poor fit, such as the presence of mixed stands or mixed land cover types in the study area. The ratios of chosen tree cores for each species did not necessarily match the areas covered by the leaf for each species. In addition, there was the presence of mixed land covers in the study area. Brehaut et al. [19] verified that forest, tundra, and shrubland could have differing spectral signal characteristics, which could also influence the reconstructions, especially for the water areas, which will dramatically affect the final model fit.

Conclusions
From the reconstructed NDVI spatial distributions in 1906 through 2015, we found the tendency of mean NDVI was increasing from 1906 to 2015. Multiple disturbances and driving factors have influenced this forest over the past century, with drought dominating the vegetation dynamics. Similarly, browsing intensity determined by ungulate population density could be attributed to the regime shift in vegetation density found in our study.
Our models removed noise and redundancy, chose three sub-models for each land cover, and added the time lag into the equation. However, there are still some constraints that need to be considered before the model can be applied to other areas. The target areas should respond to a single and persistent stress, and during the stress, foliage and treering width should grow synchronously and react similarly to the stress. If the satellite images used in the research are coarse resolution images, we need to consider the effects of the mixed forest species and mixed land cover, especially for the water areas.
Our model simulated the variations and made NDVI annual maps with fine resolution (30 m) one hundred years ago, which was the first time to consistently "see" regime shifts during that time. The reconstructed NDVI even rectified our assumption that the NDVI in our study areas was decreasing. The annual NDVI series of maps  could vividly display how some ecological disturbances modify the local forest construction and lead to vegetation succession. As the model could be applied in more areas, the effects of global climate change on forest ecosystems may show up, and the anthropogenic effects on the whole forest ecosystem could be quantitatively evaluated.