Mapping the Normalized Difference Vegetation Index for the Contiguous U.S. Since 1850 Using 391 Tree-Ring Plots

Abstract

The forests and grasslands in the U.S. are vulnerable to global warming and extreme weather events. Current satellites do not provide historical vegetation density images over the long term (more than 50 years), which has restricted the documentation of key ecological processes and their resultant responses over decades due to the absence of large-scale and long-term monitoring studies. We performed point-by-point regression and collected data from 391 tree-ring plots to reconstruct the annual normalized difference vegetation index (NDVI) time-series maps for the contiguous U.S. from 1850 to 2010. Among three machine learning approaches for regressions—Support Vector Machine (SVM), General Regression Neural Network (GRNN), and Random Forest (RF)—we chose GRNN regression to simulate the annual NDVI with lowest Root Mean Square Error (RMSE) and highest adjusted R2. From the Little Ice Age to the present, the NDVI increased by 6.73% across the contiguous U.S., except during some extreme events such as the Dust Bowl drought, during which the averaged NDVI decreased, particularly in New Mexico. The NDVI trend was positive in the Northern Forest, Tropical Humid Forest, Northern West Forest Mountains, Marin West Coast Forests, and Mediterranean California, while other ecoregions showed a negative trend. At the state level, Washington and Louisiana had significantly positive correlations with temperature (p < 0.05). Washington had a significantly negative correlation with precipitation (p < 0.05), whereas Oklahoma had a significantly positive correlation (p < 0.05) with precipitation. This study provides insights into the spatial distribution of paleo-vegetation and its climate drivers. This study is the first to attempt a national-scale reconstruction of the NDVI over such a long period (151 years) using tree rings and machine learning.

1. Introduction

The global cumulative net CO₂ emissions from 1850 to 1989 were 1400 Gt, while from 2010 to 2019, CO₂ emissions reached 410 Gt, meaning that emissions in the recent decade were around one-third of those in the previous 140 years [1]. Meanwhile, global warming and the increasing frequency of extreme weather events have severely affected forests and grasslands, as well as social and economic development [2,3]. Vulnerable biological ecosystems in the U.S. have suffered from excessive environmental degradation [4,5]. Researchers [6] investigated tree mortality in four national forest parks in California and Nevada and found out that around 48.9% of trees died during a three-year-long severe drought. Before and after Winter Storm Uri in 2021, researchers [7] compared the abundance of bufferflies in Willacy and Kenedy Counties in Texas and highlighted that the extreme cold temperatures caused a significant decrease (p < 0.001) in the abundance of local butterflies. The short-term influences on the local ecosystem can be evaluated accurately, but the long-term effect on large biological ecosystems requires an extended historical spatial database.

Even though optical sensors can monitor ecological changes across the country and present clear and reliable vegetation metrics, most of them fail to display the spatial distribution of vegetation density prior to the 1970s, when their carriers were not launched [8]. To overcome this gap, tree-ring data stand as ideal proxies to retrieve annual time series of national vegetation maps over the past 150 years with moderate spatial resolution. There are three reasons to use tree rings to reconstruct the spatial distribution of vegetation densities. Firstly, trees have long lifespans and are sensitive to climate or ecological changes. It is common to see trees over a century old. The ages of some Bristlecone pines (Pinus aristata) in the White Mountains, California can be over 4000 years [9]. Secondly, tree rings have well-resolved annual resolution [10]. The consistent ring index could support annual time-series reconstruction. Lastly, forests predominantly cover 2,965,619 km² in the contiguous U.S. in 2010 [11].

Numerous paleoclimate research studies used tree rings as indicators to explore past climates [12,13]. Researchers [14] reconstructed annual temperatures in North America using maximum latewood density from conifer trees and identified the summer temperature of 1810 as the lowest since 1760. Researchers [10] reconstructed the Palmer Drought Severity Index (PDSI), a drought index, for all of North America from A.D. 951, emphasizing that the greatest mega-drought occurred in the 12th century. Tree rings serve not only as a proxy for climate but also for ecological processes [15,16]. Researchers [17] collected tree cores from four sites in the Alatau Mountain, Kazakhstan, and reconstructed the normalized difference vegetation index (NDVI) during the growing season and ring-width index (RWI) using simple linear regressions. Similarly, researchers [18] discovered a positive linear relationship between the grassland NDVI and RWI in eastern Inner Mongolia, China, reporting that early growing season precipitation was crucial for capturing signals of grassland growth reflected in tree-ring width. The NDVI of one study area in those studies was simplified as one point, and Li et al. [19] went further, reconstructing the annual NDVI maps in the Greater Yellowstone Ecosystem with 15 tree-ring plots. They did not expand their study areas to a larger size due to lack of tree-ring collection.

However, few studies refer to the reconstruction of vegetation metrics with moderate resolution using tree-ring data due to two challenges. Firstly, the training data are limited. Our NDVI data source was derived from satellite data from 1981 to 2010, which provided us with only 30 years of data to train machine learning algorithms. Therefore, we chose three machine learning approaches which are highly adaptive to small training sets. Secondly, the number of tree-ring plots is limited. Following previous studies [20], we used point-by-point regression to reconstruct a time-series NDVI map at a moderate scale (5.5 km resolution), which requires a substantial number of tree-ring plots. We retrieved available ring-width data from tree stands across the contiguous U.S. from the International Tree-Ring Data Bank (ITRDB) and collected the additional cores ourselves. We conducted two rounds of reconstruction to simulate a longer time series of annual NDVI maps with acceptable accuracies. With these NDVI time-series maps, we can vividly observe variations in vegetation across the contiguous U.S. since the Little Ice Age (1850).

Our study on NDVI reconstruction was guided by three hypotheses. Firstly, we hypothesized that our model could simulate the NDVI values across the contiguous U.S. with tolerable errors. Meanwhile, developing a model to fit millions of pixels over a large spatial extent with heterogeneous landscape was a challenge. Secondly, we hypothesized that there are significant temporal and spatial differences in the NDVI between the present and 150 years ago across the contiguous U.S. While some studies reported regional vegetation shifts [20,21], it is questionable whether vegetation density has changed at a national scale over 150 years. Lastly, we hypothesized that temperature or precipitation are drivers of the NDVI in some states. Some fine-scale research has illustrated the influence of precipitation and temperature influences on the local NDVI [22]. However, we need more evidence to confirm the effects of precipitation and temperature on vegetation across the U.S..

2. Study Areas and Data Sources

Our study area covers contiguous U.S., covering 3,119,884.69 square miles (8,080,464.3 km², Figure 1). The study site ranges from 49.38N° to 25.84N° in latitude and from 66.95W° to 124.67W° in longitude. There are 10 ecoregions (level 1) in the contiguous U.S.: Eastern Temperate Forests, Great Plains, Marine West Coast Forests, Mediterranean California, North American Deserts, Northern Forests, Northwestern Forested Mountains, Southern Semi-Arid Highlands, Temperate Sierras, and Tropical Humid Forests. Our study period ranged from 1850 to 2010. We chose 1850 as the starting year to compare vegetation densities during the Little Ice Age, when periods of low temperatures were prevalent, with the present. For 161 years (1850–2010), the U.S. experienced warming temperatures, substantial urbanization, and industrialization. We intended to observe how vegetation shifted in the changeable era. We used two data sources: tree-ring data and vegetation remote sensing data.

The tree-ring data (391 plots) were retrieved from the international tree-ring data bank (ITRDB) and our field collection (Figure 1). Some of our collections were younger than 1850, so we divided our data into two sub-periods: Sub-period 1 (1850–2010) and Sub-period 2 (1955–2010). If tree-ring data were available for 1850 or earlier, the plot was assigned to both sub-periods (Sub-period 1 and Sub-period 2). If the tree-ring data were only available between 1850 and 1955, the plot was assigned to Sub-period 2. The spatial distribution of the involved tree stands was not even. In the Midwest, there were four states (South Dakota, Nebraska, Kansas, and Texas) without available plots. In those states, some pixels with long distances to the surrounding tree-ring plots might suffer a poor reconstruction because the climate (temperature and precipitation) of the target pixel and the tree-ring plots might be different. We used the NDVI as our vegetation density indicator to access vegetation health and growth. The remote sensing data source was the Advanced Very High-Resolution Radiometer (AVHRR) NDVI version 5 data derived from the NOAA AVHRR Surface Reflectance product (https://catalog.data.gov/dataset/noaa-climate-data-record-cdr-of-avhrr-normalized-difference-vegetation-index-ndvi-version-52, accessed on 5 May 2024). The NDVI dataset spanned from 1981 to 2010, and its temporal and spatial resolutions were one day and 5.5 km, respectively. We averaged the available NDVI images for the growing season (June, July, and August) as our annual NDVI data. Climate factors such as precipitation and temperature could affect vegetation growing (vegetation density), especially in the growing season; thus, the monthly temperature and precipitation data for the contiguous U.S. and each state from 1895 to 2010 were used from the National Centers for Environmental Information (https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/national/time-series, accessed on 5 May 2024). Before 1895, no instrumental records of temperature, precipitation, or PDSI records at the national level were available.

3. Methods

The whole study is built on the relationship between the vegetation densities and tree radial growths (Figure 2A). During the extreme weather periods, precipitation and temperature limit the photosynthesis process in chloroplasts (grasslands and forests) causing the vegetation to generate less energy (carbohydrate) which results in reduced leaves (or grass blade) growth and thus less vegetation density. Furthermore, the energy deficit also impairs tree radial growth resulting in the narrower tree-ring widths [23]. Hence, energy produced during photosynthesis is crucial to support a higher density of leaves, grass growth, and substantial ring growth. Vegetation density and tree radial growth could be two reliable proxies for the climate, where the two proxies could have a high correlation with each other, mediated by climate conditions. Based on the relationship, we used point-by-point regression [19] to reconstruct vegetation density (NDVI) with surrounding tree radial growth (RWI, Figure 2B). We computed the correlations between climate factors (temperature and precipitation) and the NDVI. We explored the spatial–temporal variations in the relationship between the NDVI and climate factors for the contiguous U.S. from 1850 to 2010 and examined the NDVI drivers for each state.

3.1. Tree-Ring Data Collection and Processing

We collected a total of 391 plots, of which 390 plots were from the ITRDB and 1 plot from our field collection (Figure 1). In our stand, we selected 10 trees that were closest to the center of the plot with diameters at breast height (DBH) larger than 100 mm. We extracted two cores from each tree using an increment bore at the height of 0.3 m [24]. In one plot, we collected a total of 20 tree cores from 10 trees which were then dried, glued to a wooden mount, and sanded progressively with different grades of sand papers (800, 400, and 150 grit). We measured and dated the tree-ring width with the Velmex measuring system (Velmex, Inc., Bloomfield, NY, USA) with a precision of 0.001 mm. The COFECHA software version 6.06P, (Laboratory of Tree-Ring Research, University of Arizona, Tucson, Arizona, AZ, USA) was used to cross date and check the measurements. The raw ring width from our field collection and from ITRDB (https://www.ncei.noaa.gov/products/paleoclimatology/tree-ring, accessed on 2 May 2024) were detrended using ARSTAN version 48 (Win, Tree-Ring Laboratory, Lamont Doherty Earth Observatory of Columbia University, New York, NY, USA) [25] with an age-dependent spline. We used the standard tree-ring index as the metric to describe tree width.

If the Expressed Population Signal (EPS) of plots were below 0.85, the ring-width indices (RWI) were cut off. After all tree-ring plots (391) were processed, we had 187 tree-ring stands in Sub-period 1 and 391 plots were in Sub-period 2. The detailed information of all collected plots could be found in the Supplementary Table S1. Among 391 plots, the 5 most frequent species in our study (the numbers of the involved plots are in the parentheses) included: Ponderosa pine (Pinus ponderosa, 39), Bur oak (Quercus macrocarpa, 32), Douglas-fir (Pseudotsuga menziesii, 26), White oak (Quercus alba, 24), and Eastern hemlock (Tsuga canadensis, 24).

3.2. Remote Sensing Data Processing

Tree-ring width has an annual resolution, so we averaged all the available NDVI images during a growing season from 1 June to 30 August (growing season). We collected NDVI images from 1981 to 2010 using the Google Earth Engine (https://earthengine.google.com/, accessed on 2 May 2024). Though the NDVI images after 2010 were still available, we excluded them because our tree-ring index collections were measured earlier. We also downloaded the 2010 land use and land cover classification map from Moderate Resolution Imaging Spectroradiometer (MODIS) product, MCD12C1(https://lpdaac.usgs.gov/products/mcd12c1v006/, accessed on 2 May 2024), whose spatial resolution was 5.5 km. We reclassified the contiguous U.S. into four land covers: water (including wetland), forest, grassland (including cropland and prairie), and impervious areas (including desert, built-up land, and barren land). Previous studies found out that the relationships between tree-ring indices and the NDVI were ambiguous in water and impervious areas [19], so we only built regressions for the forests and grasslands. During our study time (1850–2010), the U.S. experienced substantial land use changes with intense urbanization and industrialization, which rarely results in the impervious areas reverting back to forests or grasslands. Therefore, if one pixel was classified as forest or grassland in 2010, it was highly likely to be permanent forest or grassland during our study time.

3.3. Point-by-Point Regression

Researchers [19] reconstruct the Palmer Drought Severity Index (PDSI) using the point-by-point regression where the entire study areas are divided into grids. The PDSI in each grid was reconstructed using signals from its neighboring tree stands which showed the relationship between PDSI and ring-width index (RWI). For instance, when there is adequate precipitation and suitable temperature (high PDSI), vegetation, including forests and grasslands, could grow well, thus resulting in better vegetation density (NDVI) and ring-width indices (RWI) of surrounding plots. However, in the presence of some extreme climate events, during heat stress or lack of precipitation in the growing season (low PDSI), vegetation density and ring width in the plots are highly affected. We assume the relationship between the NDVI and RWI could be represented by Equation (1) and we intend to build the model for each pixel following this relationship.

NDVI = f (RWI1, RWI2, RWI3, … RWIn)

(1)

3.3.1. Nested Reconstruction (Defining a Temporal Range)

In general, the reconstruction only has one round of training where all selected RWI shared one common simulation period. Our study applied multiple regressions (two rounds) and divided our study time into two sub-periods: 1850–2010 and 1955–2010 to include as many plots and RWI of each plot as possible to ensure a long time series of annual reconstructed NDVI maps with the least number of errors. We had two training sets with 187 tree-ring stands in Sub-period 1 (1850–2010) and all 391 plots in Sub-period 2 (1951–2010). Both training sets had a common training period (1981–2010). During this period, every pixel (5.5 km × 5.5 km) had annual NDVI values, and each selected plot had annual RWI.

3.3.2. Choosing Surrounding Plots (Defining a Spatial Range)

From Equation (1), we assumed that the NDVI values in any vegetation pixel could be retrieved from the RWIs of the surrounding tree stands, so we should declare the definition of the “surrounding” tree stands. We made sensitivity maps with various searching radii: 825 km (150 pixels × 5.5 km), 1100 km (200 pixels × 5.5 km), and 1375 km (250 pixels × 5.5 km) for 1850–2010 (187 plots, Figure 3). Even for the most remote pixels, they could still find at least one surrounding plot within a 1100 km radius (Figure 3B). If the radius was too small, the number of the involved tree stands would also be small in some Midwest states. The plot density in the Midwest (especially in Nebraska and Texas, Figure 3A) was small and the reconstruction with less stands might impair model accuracy. If the range was too large, the NDVI in the target pixel and the tree stands might face different ecological environments and its correlation could be very weak (Figure 3C). Therefore, 1100 km was the best option for the searching radius (Figure 3B). To tease out the plots far from the target pixel, we ranked the plots from spatially closest to furthest away. If the number of surrounding plots was more than 20 within the search radius, we only chose the 20 closest plots to the target pixel.

3.4. Regression Selection and Evaluation

Our study chose three machine learning approaches as our regression algorithms: Support Vector Machine (SVM) regression, General Regression Neural Network (GRNN) regression, and Random Forest (RF) regression.

The SVM regression is a widespread machine learning approach put forward by previous studies [26]. The regression sets a hyperplane in high dimension and searches the best fits for the data in a continuous space. The approach has a promising performance when the size of the training set is small, which exactly fits our scenario where for each pixel regression, we only had 30 pairs (1981–2010) of NDVI values and RWIs in our training set. We assume that SVM might perform well with a limited training dataset. The GRNN regression is an improved nonparametric neural network created by researchers [27]. Unlike traditional neural networks with iterative training, the algorithm only has one single pass learning with the Gaussian functions, which efficiently trained our non-linear models. The model can also keep its high accuracy and resilience with some noise, and like SVM, the model performance cannot be affected by a small sampling size, which exactly matched our situation. RF regression is an ensemble approach combing multiple regression trees put forward by previous studies [28]. The algorithm generates multiple uncorrelated regression trees (in our study, we generated 100 trees with the bagging approach) and corrects the overfitting of the trees, which is adaptable to small dataset. Our result is the average of the 100 regression trees.

We used four validation metrics to comprehensively evaluate the three model performances in two time periods individually: Root Mean Squared Error (RMSE), Mean Absolute Percentage Error (MAPE), Mean Absolute Error (MAE), and adjusted R2. Though adjusted R2 was widely applied on many model validations, researchers [29] reported that the metrics could not have a reasonable evaluation on non-linear algorithm. We used five-fold validations to validate model performances where the training set was split into five subsets. One year in the subset was used as a validation year and the remaining 80% of data was used for training (

Table 1. Validation scheme for five-fold validation.

Order	Training Set (24)	Validation Set (6)	Validation Year
Fold 1	1987–2010	1981–1986	1985
Fold 2	1981–1986, 1993–2010	1987–1992	1990
Fold 3	1981–1993, 1999–2010	1993–1998	1995
Fold 4	1981–1998, 2005–2010	1999–2004	2000
Fold 5	1981–2004	2005–2010	2005

Note: Validation year is one of the years in the validation sets.

).

In summary, our research was a pixel-level reconstruction. We picked a vegetation pixel I, the annual NDVI values of which were from 1981 to 2010 (30 years) as our dependent variable. Within the circle, the center and radius of which were the pixel i and 1100 km, respectively, we chose all tree-ring indices within the circle as our independent variables. The plots with various available periods were classified into two groups: 1850–2010 and 1955–2010. The NDVI values and tree-ring indices were both available from 1981 to 2010 as our training set. Three approaches (RF, SVM and GRNN) were applied on pixel i and the evaluation metrics (RMSE, MAPE, MAE, and adjusted R²) of three approaches were calculated for the pixel i. We ran the same process on each vegetation pixel with three models. Then, we averaged the metrics and chose the best model from the three.

3.5. Statistics Analysis

We averaged the NDVI in the growing season (in the latter part, the annual value indicates the averaged data in every growing season) for each state and each ecoregion and mapped the reconstructed and real time series of NDVI images. The objective of making the averaged NDVI map is to display the general vegetation pattern for the whole U.S., U.S. which is more representative than any annual maps in one specific year, because the averaged value could avoid the influences of extreme climate or natural disasters. We also calculated annual temperature and precipitation for the contiguous U.S. and each state from 1895 to 2010. We intended to compare the vegetation density between the Little Ice Age and the present. To avoid the influences from one or two specific years with extreme weather, we averaged the NDVI maps for the first 10 years (1850–1859) and the last 10 years (2001–2010) using Equation (2).

Change rate = (NDVI 2001–2010−NDVI 1850–1859)/NDVI 1850–1859

(2)

We also computed the tendencies for each ecoregion and correlations between climate factors (precipitation and temperature) and the NDVI for each state with five-year-averaged NDVI values.

4. Results

4.1. Regression Model Performances

We analyzed 2 sub-periods: Sub-period 1 (1850–2010) with 187 plots and Sub-period 2 (1955–2010) with 391 plots. We conducted three regression models for each pixel in both sub-periods. Due to space constraints, we only displayed the scatter plots of Fold 1 (Year 1985) between the observed and the simulated (Figure 4). However, the evaluation metrics in the figures were the average of five-fold validation. The detailed information for each fold can be found in the Supplementary Table S2. Within Sub-period 1 (1850–2010), the GRNN was the best among all with the lowest RMSE (0.0392), lowest MAPE (11.58%), lowest MAE (0.0309), and highest adjusted R² (0.87). The SVM regression model was the worst. A similar story was also found in the Sub-period 2 (1955–2010) where the GRNN model outperformed the rest with the lowest RMSE (0.0375), lowest MAPE (11.05%), lowest MAE (0.0299), and highest adjusted R² (0.87).

Comparing the GRNN in Sub-period 1 and Sub-period 2, all metrics (RMSE, MAPE, MAE, and adjusted R2) in Sub-period 2 (1955–2010) were better than those in Sub-period 1 (1850–2010). However, the errors in GRNN Sub-period 1 were still acceptable. So, we divided our entire study period into three parts: 1850–1954, 1955–1980, and 1981–2010. From 1850 to 1954, our NDVI database used reconstructed annual NDVI maps from the Sub-period 1 GRNN model. From 1955 to 1980, our NDVI database used the reconstructed maps from the Sub-period 2 GRNN model. From 1981 to 2010, we directly used the real data from the AVHRR NDVI.

4.2. NDVI Temporal Differences in the U.S.

With the averaged NDVI maps overlaid by the annual NDVI maps from 1850 to 2010 (Figure 5A), we could divide the spatial distribution of the NDVI across the contiguous U.S. into three parts. The western part around the Pacific Ocean (including California, Oregon, and Washington) had dense vegetation around the mountains. The mountainous part between the Rocky Mountains and the Great Plains had comparatively low vegetation densities; however, the Greater Yellowstone areas, especially southeast Idaho, remained green. The eastern part had an NDVI increase as the gradient flattened from west to east. From the Little Ice Age to the present, the overall NDVI increased by 6.73% (mean NDVI 1850–1859 = 0.297, mean NDVI 2001–2010 = 0.317, Figure 5B). However, two regions had substantial NDVI decreases: southern New Mexico—northern Texas and Florida. The largest NDVI increases were in northern New Mexico and western Kansas, where the NDVI increases in some areas went beyond 20% (Figure 5B). The averaged NDVI map and change rate map only described the general temporal changes; however, annual variations in temperature and precipitation caused year-to-year NDVI pattern changes. In the year 1850 (the Little Ice Age period), some areas in Montana marked with a blue circle had lower vegetation density, compared to other selected years (Figure 6A). In the year 1933 (the Dust Bowl Drought, one of the most severe droughts in the U.S. history), the NDVI in New Mexico marked with a blue circle was lower than the rest of the years (Figure 6B); meanwhile, the overall NDVI in 1933 was also one of the lowest values in the simulated period (1850–1980, Figure 6E). In the year 1965, the temperature and precipitation in most areas were suitable for vegetation growth (PDSI = 1.27) but in Florida, some low-density vegetation areas stood out (Figure 6C). In the year 2010, sunlight and precipitation were abundant (PDSI = 3.21), while low-density vegetation (the blue circle in Figure 6D) disappeared in Florida.

4.3. NDVI Spatial Differences in the Contiguous U.S.

In addition to temporal differences, the spatial differences of the NDVI were noticeable in the contiguous U.S.. We divided our study areas into 10 ecoregions (Figure 1). The national NDVI tendency was increasing but there were five ecoregions with decreasing tendencies: Northern Forest (p = 0.015), Tropical Humid Forest (p = 0.002), Northern West Forest Mountains, Marin West Coast Forests, and Mediterranean California (Figure 7). In the 1850–1865 period, the NDVI density was relatively high in the Tropical Humid Forest ecoregion, compared to other years. The other five ecoregions showed increasing tendencies: Eastern Temperate Forests, Great Plains, North American Deserts (p = 0.035), Temperate Sierras (p = 0.048), and Southern Semi-Arid Highlands. Even though the number of increasing and decreasing tendencies were half–half, the three largest ecoregions (Eastern Temperate Forests, Great Plains, and North American Deserts, Figure 1) almost determined the overall tendency (positive) of the entire country.

4.4. Correlations Between Climate Factors and NDVI

We had two climate factors: temperature and precipitation. If the absolute value of R (NDVI and climate factors) in the Pearson correlation was 0.3 or higher [30], we assumed that vegetation in that state was sensitive to the corresponding climate factor. There were six temperature-sensitive states: Louisiana (R = 0.46, p = 0.03), Washington (R = 0.41, p = 0.05), Illinois (R = 0.39), Idaho (R =0.35), Oregon (R = 0.31), and Texas (R = 0.31), whose correlations were all positive (Figure 8A). There were also six precipitation-sensitive states: Oklahoma (R = 0.49, p = 0.02), Washington (R = −0.47, p = 0.02), California (R = −0.36), Indiana (R = 0.36), Illinois (R = 0.31), and Michigan (R = −0.30). Half of them had positive influences on the NDVI while the other half had negative influences (Figure 8B). Washington was the only state whose vegetation growing had strong correlations with both climate factors (temperature and precipitation).

5. Discussion

5.1. Evaluating Our NDVI Models

We performed pixel-by-pixel regression on each vegetation pixel. The conventional regression encountered two problems: choosing the reconstruction period for each pixel and spatial surrounding tree-ring plots for each pixel.

Classical regression typically runs one round of regression, which was suitable for the tree cores collected in the same year within the same available period. However, most of our tree-ring plots were from the ITRDB dataset, so the available periods of the plots vary from other research projects. If we set our study time as being only in the 1955–2010 period, we would lose the RWI values from plots earlier than 1955. If we set our study time as being only for the 1850–2010 period, we would lose the plots between 1850 and 1955. So, we ran a two-round regression for the periods of 1850–2010 and 1955–2010. In the 1850–2010 period, we could explore earlier NDVI distributions with fewer plots, while in the 1955–2010 period, we could have shorter NDVI reconstructed periods with more plots. The results also supported our assumption that the 1955–2010 reconstruction with 319 plots had less errors than the 1850–2010 reconstruction with 181 plots (Figure 4), so the annual simulated NDVI maps from 1955 to 1980 using 391 plots and GRNN were chosen as our primary approach. Even though the simulated NDVI from 1850 to 1954 had higher RMSEs than those from 1955 to 1980, their performance was still acceptable. The SVM model, the worst model among the three, is very sensitive to key hyperparameters whose poor performance may result from the selected kernel function.

The other issue is to choose the surrounding plots for each pixel. Previous studies preferred to choose fewer plots for each pixel, but the distances were close to the target pixel. Researchers [31] reconstructed PDSI for North America since 1300. The data sources came from weather stations which recorded early historical records, allowing long reconstruction periods. However, the NDVI images were available since the 1980s, so we need more available plots to simulate the NDVI. The extremely long search radius (1100 km) is for the pixels in southern Texas. We only chose the closest 20 plots, although there were numerous plots within the radius. Even with the extreme 1100 km radius, there was more than 42.55% of pixels where the closest 20 plots were within 550 km.

5.2. Spatial–Temporal NDVI Variations in the U.S.

The objective of making an averaged NDVI map is to display the general vegetation pattern for the whole U.S., which is more representative than the map in one specific year because the averaged value could avoid the influences of extreme climate or natural disasters. There are no significant spatial differences between the averaged NDVI from 160 annual maps and the recent NDVI map. The vegetation density was low in the Rocky Mountains, while there was high vegetation density in the Pacific west coastal areas and the US east (Figure 5A). For the overall contiguous U.S., greenness has increased by 6.73% since 1850 (Figure 6E and Figure 7). However, some extreme weather events led to some inter-annual vegetation variations during our study period.

There were two extreme national events from 1850 to 2010: the Little Ice Age and the Dust Bowl, largely affecting the vegetation biomass and density across the country, though the magnitude of the effect varied in different regions.

In the Little Ice Age period, the vegetation density in the mountainous northern U.S. was very low, as highlighted in Figure 6A. Researchers [32] had similar observations with aerial photos in 1870 and found that the tree densities around the tree line in 1870 were only half of the modern densities. We speculated that the low density and higher vegetation mortality resulted from the cooler summers and shorter growing seasons. Previous studies [33]’s research also confirmed our speculation, and they showed that tree establishment and regeneration were highly affected by cool temperatures and frost. However, the Little Ice Age had positive contributions on vegetation growth in Florida (Figure 6A). With pollen accumulations, ice cores, and other proxies in eastern Mesoamerica, researchers [34] had a similar finding: that there was a tropical forest expansion during the Little Ice Age. The best explanation for the vegetation density increase could be the increased meridional flow led by the Little Ice Age that released the local drought and brought more winter precipitation. In the Southern Short Grass Prairie (New Mexico and Texas), we observed substantial vegetation degradation (Figure 5B). The Dust Bowl drought was not only a megadrought with extreme water shortages but also brought secondary disasters to local vegetation. Previous studies [35] supported our findings, and they observed the extreme degradation in Quay County, New Mexico where there is a lack of precipitation and unstable sand dunes with spoiled vegetation. Wind erosion eroded the fine soil in the Great Plains and caused nutrient loss via soil leaching which also affected vegetation growth [36]. The results from Pompa-Garcia et. al.’s study [37] using the tree-ring index pointed out that the reconstructed NDVI in 1971 and 1972 were pretty low, which matched our results in Figure 6E. Their study also attributed the poor growth to local water shortages. However, in 2010, temperature and moisture were more suitable for vegetation growth (PDSI = 3.21) and low-density vegetation disappeared in Florida (Figure 6C,D).

Apart from the NDVI inter-annual changes, our research also pinpointed the spatial change in the 10 ecoregions. Five of them had decreasing vegetation density tendencies, especially for the Northern Forests and Tropical Humid Forests with p values less than 0.05 for both ecoregions. Researchers [38] corroborated our findings that forest cover in Northern Forests decreased and the rest of the forests fragmented with Landsat and MODIS satellite images. They asserted that the real estate market was the main cause of deforestation. The vegetation loss may not only change wildlife communities but also affect local ecosystems [39]. The largest three ecoregions had increasing tendencies, which significantly influenced the overall greenness. Among them, North American Desert had a significant increasing tendency (p = 0.038). We speculated that the general NDVI increase came from climate warming and agriculture expansion (in our study, agriculture land belongs to grassland). Present surface temperatures are around 2 degrees Celsius warmer than the Little Ice Age’s temperature [40]. Warmer and drier conditions also prompted the regeneration and recruitment of the young whitebark pine and Engelmann spruce after the Little Ice Age [33]. Previous studies [41] suggested that from 1850 to 1920 and from 1920 to 2020, there were land cover changes of 59.01 mega hectares and 11.40 mega hectares, from shrub land to cropland, around the contiguous U.S. The highest land cover change intensities were observed in the South Central, North Central, and the Great Plains.

5.3. Climate Drivers of NDVI Change

In general, the NDVI values in most states showed positive correlations with temperature in the growing season while the correlations between NDVI and precipitation varied across the country. The NDVI in Louisiana (p = 0.03) and Washington (p = 0.05) showed positive correlations with temperature in the growing season. It is highly likely that those states were situated around the coastal areas, where marine layers caused by strong ocean upwelling could develop fog or overcast weather [42]. In all seasons, cloudy days were always cooler than clear days, and there were more than 200 cloudy days in western Washington in one year [43].

The NDVI in Washington (p = 0.02), Oregon, California, and Michigan had negative correlations with the precipitations in the growing season. Those states had long coastlines and were close to the Pacific Ocean or the Great Lakes with mild temperatures and abundant precipitations. Therefore, the drought frequencies in those states were low, especially for western Washington and western Oregon [44]. Excessive moisture could reduce oxygen in the soil and destroy the plant root system, which increased vegetation mortality. Positive correlations between the NDVI and precipitation were observed in the states of Oklahoma (p = 0.02) and Indiana. A similar study from researchers [45] supported our results and ranked Oklahoma as the most drought-vulnerable state in the U.S. There were five major droughts in the state from 1901 to 2014, where prolonged dry weather reduced crop growth and livestock production [46]. Researchers [47] believed the drought in Indiana could worsen and predicted that with the increased frequency of drought, Indiana would lose 8–21% of corn and soybean yields in 2100.

5.4. Limitation and Future Studies

There are three limitations to our study that are worth mentioning, and we hope further studies can overcome these issues. Firstly, the spatial distribution of the selected plots was uneven because most of the plots were from the ITRDB, where numerous studies with various research objectives uploaded their results online. States like Kansas, Nebraska, and South Dakota did not have available plots in the ITRDB. To improve our model’s accuracy, it is recommended to collect more plots from those states. Secondly, the simulated results could not capture extreme NDVI values. When we applied three machine learning approaches to reconstruct NDVI values, the approaches heavily rely on training samples with extreme values and a large sample size. They could not capture the extreme values, even when the training set had similar data. To effectively reconstruct extreme NDVI values, we could get rid of the training-prediction process, a general workflow of a machine learning model. We recommend that further studies try multiple linear regressions, which were widely used in dendrochronology reconstruction [48,49]. Lastly, our study did not include the role of anthropogenic activities. In the last 50 years, anthropogenic activities have been one of the main drivers altering vegetation patterns in the U.S., but there are no available data on the annual human spatial distribution over our study period. A higher population usually brought urban sprawl and farmland increase, which are the main causes of land use change. New housing and more food demands switched forest into urban areas and agriculture. If we could accurately evaluate the population increase in the long-term, we may successfully assess the land cover change as a population index. Then, we may input it into our model to prompt our model. We also realized that it is meaningful to explore and deepen the analysis between the NDVI and climate factors (PDSI, temperature and precipitation) at pixel level, so we plan to simulate the annual land surface temperature and moisture maps with available remote sensing images and tree-ring plots.

6. Conclusions

We took full advantage of 391 plots and 3 machine learning approaches to reconstruct annual NDVI maps since 1850. After validations for two periods, we found out that GRNN was the best and most robust approach from 1850 to 1980, with the least error and highest adjusted R². To improve model accuracy, we plan to collect more tree cores, especially in South Dakota, Nebraska, Kansas, and Texas. With the simulated annual NDVI maps for the contiguous U.S., we observed an overall NDVI increase of 6.73%. Temporal differences in the NDVI existed, especially during the Little Ice Age and the Dust Bowl drought. The simulated NDVI values for Montana in 1895 and Florida in 1933 showed low vegetation clusters, but by 2010, those low vegetation densities had disappeared. There were also spatial differences in the contiguous U.S. Among ten ecoregions, five had NDVI increasing tendencies (Northern Forest Tropical Humid Forest, Northern West Forest Mountains, Marin West Coast Forests, and Mediterranean California), while the rest had decreasing tendencies (Eastern Temperate Forests, Great Plains, North American Deserts, Temperate Sierras, and Southern Semi-Arid Highlands). Washington and Louisiana had significantly positive correlations with temperature (both p values were less than 0.05), likely because cloudy and overcast days caused by marine layers might cool daylight temperatures, making these states prefer warmer temperatures. Oklahoma and Indiana, being inland states, showed strong positive correlations with precipitation, while coastal states like Washington, Oregon, California, and Michigan displayed negative correlations with precipitation, where abundant water could impair vegetation growth.

Our findings illustrate how the NDVI changed for over 150 years and provide insights into the relationship between the NDVI and climate drivers (temperature and precipitation) using a reliable model. With the annual simulated NDVI maps, this is the first time we can vividly depict the long-term impact of events on vegetation across the entire country.