Phenological Shifts of the Deciduous Forests and Their Responses to Climate Variations in North America

Abstract

Forests play a vital role in sequestering carbon dioxide from the atmosphere. Vegetation phenology is sensitive to climate changes and natural environments. Exploring the patterns in phenological events of the forests can provide useful insights for understanding the dynamics of vegetation growth and their responses to climate variations. Deciduous forest in North America is an important part of global forests. Here we apply time-series remote sensing imagery to map the critical dates of vegetation phenological events, including the start of season (SOS), end of season (EOS), and growth length (GL) of the deciduous forests in North America during the past two decades. The findings show that the SOS and EOS present considerable spatial and temporal variations. Earlier SOS, delayed EOS, and therefore extended GL are detected in a large part of the study area from temporal trend analysis over the years, though the magnitude of the trend varies at different locations. The phenological events are found to correlate to the environmental factors and the impact on the vegetation phenology from the factors is location-dependent. The findings confirm that the phenology of the deciduous forests in North America is updated such as advanced SOS and delayed EOS in the last two decades and the climate variations are likely among the driving forces for the updates. Considering that previous studies warn that shifts in vegetation phenology could reverse the role of forests as net emitters or net sinks, we suggest that forest management should be strengthened to forests that experience significant changes in the phenological events.

1. Introduction

Forests play a vital role in global carbon cycling. Forests sink the most carbon dioxide from the atmosphere and own the largest carbon pool compared to other terrestrial ecosystems such as grassland and cropland [1]. North America is one of the continents with a great part of the terrestrial lands covered by forests. Over the past few decades, the forest area in America has been reduced due to human activities [2,3]. The role of vegetation, including trees in forests, on the global environments is multifold. Vegetation sequesters carbon dioxide, and serves not only as a critical and irreplaceable biological pool, but also as a central hub for the interconnection between various spheres, including biosphere, hydrosphere, and atmosphere [4]. Studies for understanding the spatio-temporal vegetation dynamics have been carried out for a long time; those studies focus on a number of topics such as vegetation coverage, vegetation primary productivity, and vegetation landscape patterns, as well as the impact on the topics from climate changes and human activities [5,6]. Mapping the vegetation dynamics for forests is particularly important considering the contributions of forests to global carbon cycling [7].

Vegetation phenological events such as vegetation green-up, peak growth, leaf senescence and dormancy, are among the indicators to the ecological functions from vegetation [8]. Vegetation phenology can not only reflect switches in seasons or in a broader sense climate changes, but also open a window for understanding the mechanism for plants to adapt to the climate changes [9]. Vegetation phenology influences biomass productivity, biodiversity evolution, water balance and carbon balance at local, regional and global scales [10]. Previous studies show that the increase of temperature leads to significant advance of phenology in spring [11] and the delay of phenology in autumn [12], resulting in an extended plant growth season [13]. A few phenological events are within the study scope of vegetation phenology, including the start of growing season (SOS), the end of growing season (EOS) and the length of growing season [14]. For deciduous forests, the vegetation phenological events are often examined at an annual cycle and are subject to change under the influence of climatic and other environmental factors [15]. Vegetation in deciduous forests usually starts to turn green in spring, reaches peak greenness in summer and withers in autumn or early winter. However, vegetation phenological events are location-specific, as a result of the vegetation’s responses and adaptation to the environments [8,16]. Therefore, knowing the spatio-temporal patterns in vegetation phenology has both theoretical and practical implications for better managing ecosystems and understanding climate changes [17].

Several ways can be applied to make observations of vegetation phenology. Theoretically, the field-based observation is the most reliable method. For example, flux tower GPP time series data can be taken to examine the dynamics of vegetation phenology [18]. However, field-based observation is time consuming and labor intensive; it is thus not a cost-effective method to map vegetation phenology at multi-temporal points over a large area. Another difficulty from field-based observation is the determination of the phenological dates, which requires consistent standards or protocols. Unfortunately, such protocols are not widely accepted or even well documented, causing different phenological dates even from the same datasets when determined by different people at different times [10,19]. Such inconsistency in the vegetation phenology due to lack of agreed protocols could cause challenges in spatio-temporal pattern analysis on the phenology because the derived phenological dates may not be inter-comparable both spatially and temporally. Remote sensing has played an important role in understanding the temporal and spatial dynamic changes of plant phenology on a large scale [16]. Satellite remote sensing allows long-term observations on vegetation growth over a large area. Vegetation dynamics can be inversed from remote sensing imagery using vegetation indices such as the normalized distribution of vegetation index (NDVI). Once the identical protocols to extract the phenological events are applied to remotely sensed time-series imagery, it is possible to obtain spatio-temporally consistent phenological datasets for a large area and thus map the phenological dynamics from regional to global scales [20].

Forests are important to mitigate global warming. Out of all terrestrial ecosystems, forests absorb the largest part of carbon dioxide from the atmosphere as they grow and store the carbon in tree woods and soils [21,22]. Deciduous forests are characterized by typical vegetation growth cycles, including active growth in the spring, peak growth in the summer, leaf senescence in the autumn, and dormancy in the winter, and indicate to seasonal changes and climate variations. Thus, current work tries to understand the spatio-temporal patterns of the deciduous forests in North America in the past two decades (2001–2020) and explore the responses of the phenological events to the climate variations. We first use time-series NDVI data from remote sensing imagery to extract vegetation phenology of the deciduous forests. We then analyze the spatio-temporal patterns in the extracted dates of the phenological events and their dynamics against temperature and precipitation. The Google Earth Engine (GEE) platform allows efficient computation and mapping of the phenological events and their pattern analysis. The paper is organized as follows: the next section introduces the study area and data sources, the methodology will be presented in Section 3, the main findings are reported in Section 4, discussions and conclusion are given in Section 5 and Section 6, respectively.

2. Study Area and Data Sources

2.1. Study Area

Deciduous forest is found in middle-latitude regions with a temperate climate characterized by a winter season. Deciduous forests are composed primarily of trees that grow and shed their leaves in an annual cycle. The phenology of deciduous forests is sensitive to climate changes and phenological shifts can bring substantial impacts on ecosystem function and services due to their large carbon pools and fluxes [23]. Changes in spring and autumn phenology of temperate plants in recent decades have become iconic bio-indicators of rapid climate change [24]. The land area of North America is broadly covered by level II regions that are nested within 15 level I ecological regions (www.epa.gov/eco-research/ecoregions-north-america, accessed on 1 May 2022) (Figure 1). The ecological regions at level II are useful for national and subcontinental overviews of ecological patterns. The land area, particularly in its eastern part, is one of the few regions that are characterized by deciduous forests. The response of vegetation phenology to climate changes has regional differences. Vegetation in mid-to-high latitudes (between 30–70° N) is particularly sensitive to climate changes [25]. We are therefore interested in the boundary of the area between latitudes between 30–70° N within the ecoregions of North America and with coverage of deciduous forests. We only selected the ecoregions that are related to forestry, including northern forests, northwestern forested mountains, marine west coast forest, and eastern temperate forests. To extract the coverage of the deciduous forests only, the land cover data from MODIS MCD12Q1 land cover data are used. The MCD12Q1 dataset applies the International Geosphere-Biosphere Programme classification protocol which divides the deciduous forests in the land cover types labeled as #3 (deciduous needleleaf forests) and #4 (deciduous broadleaf forests). All other land cover types are masked out. The deciduous forests are primarily in the eastern part of the study area, though they also scatter across the north and western regions (Figure 1).

2.2. Data Sources

Remote sensing technology provides an effective way to map vegetation growth, and thus SOS and EOS over a large-scale area. Multiple remote sensed datasets are applied to extract the forest phenology and the associated factors such as climatic variables. The time window covers the years from 2001 to 2020. All below mentioned raster datasets, if not at the scale of 500 m, are resampled to a resolution of 500 m.

One of the frequently adopted methods using remote sensing data to measure vegetation growth is through the normalized distribution of vegetation index (NDVI), which reflects the abundance of vegetation greenness [26]. For deciduous forests, NDVI starts to increase during the spring green-up, peaks in the full growth period, and decreases after the autumn dormancy. When multi-temporal NDVI imagery is collected for each year, the trajectory of the intensity of vegetation greenness can be analyzed against the days of the year, which allows to extract the phenological events. This study takes use of MOD13A1 (ver. 6.1) product, which provides NDVI at a per pixel basis at 500 m in spatial resolution and 16 days in temporal resolution.

Climate variations have been proved to affect vegetation phenology. To examine the impact from climate factors, we consider various climatic variables and soil properties, including summed annual precipitation (PREC), annual average temperature (AAT), vapor pressure deficit (VPD), monthly minimum temperature (MINT), monthly maximum temperature (MAXT), surface shortwave radiation (SSR), and soil moisture (SM). The climatic variables have been included in many light-use efficiency vegetation productivity models like Miami (PREC and AAT) [27] and CASA (PREC, MINT, MAXT, and SSR) [28]. In addition, VPD and SM are important to model water availability to vegetation growth. TerraClimate, derived from climatically aided interpolation by combining the WorldClim dataset and the CRU Ts4.0 dataset, has been widely applied for mapping regional and global ecological parameters [29], and is used to derive most of the variables, including PREC, MINT, MAXT, VPD, SSR, and SM. However, TerraClimate does not include AAT. Therefore, AAT used in this study comes from the fifth generation from the European Centre for Medium-Range Weather Forecasts reanalysis for the global climate and weather (https://www.ecmwf.int/en/era5-land, accessed on 3 May 2022), which provides one month per pixel atmosphere temperature at about 11 km spatial resolution.

To validate the phenological events from the remote sensing imager, the flux data from six site-based stations within the international flux network (FLUXNET, https://fluxnet.org/data, accessed on 12 March 2022) are collected to derive field-based results. The consistency between the field-based and the remotely sensed SOS and EOS is then evaluated. Each of the stations collects gross primary productivity for at least 5 years, with details shown in

Table 1. Detail information of the flux sites.

Site Id	Site Name	Latitude	Longitude	Elevation (m)	Time Window
CA-Oas	Saskatchewan—Western Boreal, Mature Aspen	53.629	−106.198	530	2001–2010
US-MMS	Morgan Monroe State Forest	39.323	−86.413	275	2001–2010
US-UMB	Univ. of Mich. Biological Station	45.560	−84.714	234	2001–2010
US-Oho	Oak Openings	41.554	−83.844	230	2004–2006, 2008–2010
US-WCr	Willow Creek	45.806	−90.080	520	2001–2003, 2005–2009
US-Ha1	Harvard Forest EMS Tower	42.538	−72.172	340	2002, 2005, 2007–2010

.

3. Methodology

3.1. Extraction of the Phenological Events

The annual NDVI curve against the date of year (DOY) provides important information for obtaining the dates of phenological events for both SOS and EOS. NDVI in the deciduous forests starts to increase in spring and decreases in autumn. The dates that the NDVI shows the fastest increase during the spring green-up and the fastest decline during autumn are usually believed to be the points that vegetation reaches the onsets of active growth and senescence. The MOD13A1 dataset can be taken to build such NDVI curves from the time-series observations made on the same location across a year. Theoretically, a standard NDVI curve in natural forests is smooth along the date of the year because it is very unlikely that forest vegetation presents abrupt changes without direct human interferences. However, noises could be introduced in the remotely sensed imagery and thus the NDVI product in areas where clouds and aerosols are present during satellite Earth observations from outer space [29]. Such noises must be removed as much as possible so that the NDVI curves can be smoothed along the day of year (DOY). There are several smoothing methods intended to minimize the data noises in NDVI, including the Savizky-Golay filter [30], mean iterative filter [31], and spline interpolation-sliding window average [19]. In this study, we use the Savitzky-Golay filtering method to denoise the NDVI data because it is easy to apply and capable of reconstructing a high-quality NDVI time-series data set [32]. The Savitzky-Golay method is a weighted moving average filtering method, which depends on a polynomial degree of least squares fitting within a window [33], and its calculation formula is as follows

$${Y_j}^{*}=\frac{{{\displaystyle \sum }}_{i=-m}^{i=m}{C}_i{Y}_{j+i}}{N}$$

(1)

where Y is the original NDVI value, Y* is the smoothed NDVI value, C_i is the coefficient or weight of the ith sample which can be obtained by least squares fitting from a given high-order polynomial, m is a parameter that determines the size of the moving window, and N (=2m + 1) is the sample size in the filtering window, Y_j+i is the NDVI value for the i+jth sample, and j is the jth sample in the NDVI time series.

A cubic polynomial function is applied to interpolate the NDVI dataset. We first fit the function using the available NDIV points,

$$Y_t=a_0+a_1\ast t+a_2\ast t^2+a_3\ast t^3$$

(2)

where Y_t represents the NDVI for date t (t = 1, 16, 32, …) in DOY, a₀~a₃ are the polynomial coefficients which are fit using the least squares method from the available NDVI time series. After the coefficients of the function have been determined, NDVI for all dates in the year can be calculated.

The interpolated NDVI time series can be used to extract the dates of vegetation phenological events through threshold method [34], change point detection method [35], hybrid algorithm [14], or derivative analysis [10]. This study uses a dynamic threshold to extract the phenological events. The traditional threshold method predefines value to determine the phenological events. For example, when NDVI reaches a particular value in the spring season, the corresponding date is then assumed to be SOS. This method does not differentiate the vegetation status at various locations. Instead, a dynamic threshold method applies different threshold values to different locations. This dynamic threshold method considers both the maximum and minimum NDVI at a location, and is a location-dependent approach, which is more reasonable to obtain the threshold values [36]. For any given time-series NDVI curve representing the changes of vegetation growth within a year for a position, a location-dependent threshold is calculated as follows,

NDVI_sh = (NDVI_max − NDVI_min) × p + NDVI_min

(3)

where NDVI_sh indicates an NDVI threshold corresponding to the critical point for the SOS or EOS, NDVI_max and NDVI_min represent the maximum and minimum values of annual NDVI, respectively, the coefficient p is a ratio which is associated with vegetation cover types and is decided to be 0.5 [37], meaning both SOS and EOS are determined as the day in which the NDVI increases or drops at least 50% from the summer maximum. The date in DOY with NDVI reaching the location-dependent threshold is taken as SOS in spring or EOS in autumn for that particular location. Based on the SOS and EOS, the growing season length (GL) for the forests is defined as the total days between EOS and SOS. The three phenological events, including SOS, EOS and GL, are processed for each year during 2001–2020. To validate the result of the extracted SOS and EOS, the daily GPP data observed from the six flux stations (Table 1) are taken to derive the station-based phenological events, termed as photosynthesis phenology [38]. To make SOS and EOS from the two methods comparable, a similar extraction schema is adopted to compute the flux GPP, including Savizky-Golay filtering and dynamic threshold phenology extraction in which the coefficient p is set to be 0.15 [39], considering the difference of variables, i.e., NDVI from the remote sensed approach versus GPP from the flux stations. The SOS and EOS at the locations of the flux stations are extracted from the remote sensing approach and then compared with the results derived from the flux stations.

3.2. Phenological Trend Analysis

The Theil-Sen median regression is widely used to make a trend analysis. It is a nonparametric estimator for fitting a line to a set of points that chooses the median of the slopes of all lines through pairs of two-dimensional sample points [40]. The Theil-Sen regression method has many advantages, which include but are not limited to no requirement for any distribution assumption, robustness to outliers, and simplicity in computation [40]. The outcome from Theil-Sen regression is the Sen’s slope which takes the form of,

$$Sen’s slope=Median\left\{\frac{X_j-X_i}{j-i}:i<j\right\}$$

(4)

where x_j and x_i (i = 1, 2, …, n−1, and j = 2, 3, …, n) represent the ith and jth sample in the time series, and Median is a function for selecting the median value from a set. When Sen’s slope is greater than 0, an increasing trend in the time-series data samples can be made, while a value less than 0 indicates a decreasing trend. The Sen’s slope provides a method to assess the magnitude of change in sequence of measurements in time.

We then use the Mann-Kendall significance test (M-K test) to check if the Sen’s slope derived from Equation (4) is statistically significant. The M-K trend test is a nonparametric method widely applied to compute trend analysis for climate changes [40]. Several steps are involved for the test. First, an S statistic is computed using,

$$S={{\displaystyle \sum }}_{j=1}^{n-1}{\displaystyle \sum }_{i=j+1}^{n}f\left(x_i-x_j\right)$$

(5)

and

$$f\left(x_i-x_j\right)=\{\begin{array}{c}\begin{array}{c}-1 x_i-x_j<0\\ 0 x_i-x_j=0\end{array}\\ 1 x_i-x_j>0\end{array}$$

(6)

where x_i and x_j represent time series data respectively, and n is the number of samples in the time series. The next step is to calculate the variance of S statistics and Z value in the forms of,

$$VAR\left(S\right)=\frac{1}{18}[n\left(n-1\right)\left(2n+5\right)-{\displaystyle \sum }_{p=1}^g{t}_p\left(t_p-1\right)\left(2t_p+5\right)]$$

(7)

$$\mathrm{Z}=\{\begin{array}{c}\frac{S-1}{\sqrt{VAR\left(S\right)}} S>0\\ 0 S=0\\ \frac{S+1}{VAR\left(S\right)} S<0\end{array}$$

(8)

where n represents the number of samples in the time series data, g is the number of tied values, t_p is the number of ties for pth value. A positive Z indicates an upward trend in the time series, and a negative Z indicates a negative trend. A |Z| ≥ 1.96 indicates that there is a significant trend in the time-series data where 1.96 is the critical value of Z_1−α/2 for a probability value of 0.05 from the standard normal table [41].

3.3. Correlation Analysis between the Phenological Events and Climate Changes

The relationship between the extracted SOS, EOS, and GL and the explanatory climatic and soil variables, i.e., AAT, MINT, MAXT, PREC, SSR, VPD, and SM, was studied using correlation analysis. We aggregated each of the explanatory variables by averaging the monthly result into different time scales, including an annual step, an early season (from January to April) step, and a growing season step (from May to August), respectively. The partial Pearson correlation coefficient was then computed between the dates of each of the phenological events (SOS, EOS, or GL) and one of the aggregated explanatory variables with input samples from the years of 2001–2020.

4. Results

4.1. Spatio-Temporal Patterns in the Phenological Shifts

The SOS and EOS extracted from the remote sensing imagery (RS-based) and those from the flux stations (flux-based) are shown in Figure 2. The consistency between the flux-based and RS-based phenological events is shown in Figure 2a. Particularly, the coefficient of determination (R²) of the linear model for the SOS from the RS-based and flux-based methods reaches 0.697, suggesting that the SOS extracted from the remote sensing imagery is highly consistent with that from the flux observation (Figure 2b). The R² (=0.358) in the linear model for the EOS shows a fair agreement between the EOS derived from the RS-based and flux-based methods (Figure 2c). Therefore, the dates extracted from the remote sensed imagery can largely indicate the seasonal growth cycles of the deciduous forests.

The mean and standard deviation (std) of the phenological events including SOS, EOS, and GL in the North American forests from 2001 to 2020 are presented in Figure 3. The spatially averaged SOS is 101.61 ± 11.04 DOY (Figure 3a), EOS is 304.26 ± 11.24 DOY (Figure 3c), and GL is 202.17 ± 19.34 days (Figure 3e). The green-up date indicated by the day of the year (DOY) generally presents an earlier trend along the latitude belts from the north to the south (Figure 3a). The EOS shows a more complex pattern; for example, the areas presenting earlier EOS mainly locate in the middle and the most northern parts of the study area, while the regions in the south and north-east tend to show relatively later EOS (Figure 3c). The SOS and EOS can have annual variations, due to many factors including climate fluctuations in the years. Both SOS and EOS show considerable temporal variations as reflected from the spatial distributions of the standard deviation (std) in both SOS and EOS (Figure 3b and Figure 4d); the SOS seems more sensitive to climate changes than EOS, as the SOS shows a larger annual variation than the EOS. The mean growth season length (GL), calculated based on the SOS to EOS of each year, shows a longer duration in the northwestern and southern parts of the study area (Figure 3e), with higher annual variation compared to those of SOS and EOS (Figure 3f).

In order to reveal the trend of the phenological events in the North American forests from 2001 to 2020, this study applied the model of Sen’s slope and the M-K test. Figure 4 shows the results from the Sen’s slope analysis and M-K test. The Sen’s slope in SOS indicates that majority of the study area presents negative values, suggesting an advanced or earlier green-up in the forests along the years (Figure 4a), while the positive values in the Sen’s slope for EOS indicate a delayed dormancy trend in the forests over the years (Figure 4c). Furthermore, the M-K test suggests that the trend in nearly half of those areas presenting advanced SOS or delayed EOS is statistically significant (Figure 4b and Figure 5d). The advanced SOS and delayed EOS result in a longer growing season of the vegetation (GL, from Figure S1), and thus may enable the forests to sequester more carbon from the atmosphere, which can be regarded as an adaptive response to the increased concentration of atmospheric CO2 [42].

4.2. The Impact of Temperature on the Phenological Shifts

Temperature is commonly considered to be one of the most important environmental factors affecting the forest phenological events for both SOS and EOS [43,44]. The impact from the annual average temperature (AAT) is shown in Figure 5. Most parts in the study area present a negative correlation between AAT and SOS, except a small area located in the eastern side where higher AAT seems to cause delayed SOS (Figure 5a). Furthermore, the significance test suggests that the advancing impact from higher AAT on SOS was significant for a large part of the study area (Figure 5b). The impact from AAT on EOS is spatially differentiated. In the northern part of the study area, EOS is negatively correlated to AAT, meaning higher AAT might advance EOS. In the middle and southern parts of the study area, mixed effects from AAT on EOS, i.e., advancing or delaying effect, were observed (Figure 5c), though such patterns are not significant for a great part of the area (Figure 5d). The update in the SOS and EOS from AAT results in prolonged length of vegetation growing season (GL) in a great part of the study area during the study period (Figure S2a,b).

4.3. The Impact of Precipitation on the Phenological Shifts

Precipitation (PREC) was previously proved to affect plant phenology [45]. At the same time, the water availability in soils is not only attributed to precipitation but also to soil properties [46], evapotranspiration [47], and landforms [48]. Thus, the sensitivity of the phenological events in the forests is even more spatially varied than that from the temperature changes. The impact from the annual precipitation is either positively or negatively correlated to SOS (Figure 6a), though such impact is not significant for most parts of the study areas (Figure 6b). The impact from PREC on EOS is positive in the southern part but mostly negative in the northern areas (Figure 6c), suggesting that more precipitation benefits later dormancy in the forests in the southern part of the study area while less precipitation favors later dormancy in the northern part. Similarly, the significance test shows that the impact from PREC on EOS is not consistently strong for most parts of the study area (Figure 6d). Overall, the changes in the SOS and EOS from PREC lead to a shortened length of the vegetation growing season (GL) (Figure S2c,d).

5. Discussions

5.1. The Phenological Events in the Deciduous Forests of North America

The current study examined vegetation green-up date, senescence date and growing season length (GL) in the deciduous forests of North America. The deciduous forests in North America primarily distribute in the eastern part, which is close to the Atlantic Ocean (Figure 1). On average, the SOS in the deciduous forests starts on 101.61 (±11.04) DOY while EOS occurs on 304.26 (±11.24) DOY. The higher long-term temporal variation of SOS than that in EOS (Figure 3b vs. Figure 3d) during the years suggests that SOS is more sensitive to contributing factors (e.g., precipitation and temperature) and thus may vary considerably because of the fluctuations in all attributed factors. The spatial variations of both SOS and EOS across the study area affect the growing season length (GL) of the vegetation. On average, GL exhibits 202.17 (±19.34) days, which is very close to the growing season length of forests at the global scale [49]. The SOS shows graded belts along the latitude, meaning that the forest vegetation presents earlier green-up from the north to the south along the belts because of the increasing temperature, though the spatial variations in the other climatic factors (e.g., precipitation and sun radiation) may explain the difference of the SOS in the same latitude belt. The EOS shows less spatial variations, but the spatial patterns reveal later vegetation dormancy in the north-eastern and southern parts of the study area than the middle area located far from the Atlantic Ocean. The temporal variations in the SOS, EOS, and GL suggest that EOS is the most stable phenological event in terms of the annual changes among the three indicators.

5.2. Modification of the Phenology in the Deciduous Forests from the Climate Variations

Temperature and precipitation are among the most important climate factors that can update vegetation phenology [24]. Vegetation needs soil water to support its growth and suitable temperature to trigger the physiological metabolism [7,50]. For a large-scale area like the one examined in this study that covers different latitudes and longitudes of the globe, the climate variables may differ considerably across locations. In addition, the study area meets the Atlantic Ocean on the eastern side and the inland climate could also be profoundly affected by the vast sea area (Figure 1). For example, the areas in the north-east and south-east part of the area present earlier SOS and later EOS. This study analyzes the impact from the environmental variables on the phenological events. The study indicates that AAT shows a consistent and homogeneous impact on SOS given that the majority of the area presents significantly advanced SOS. Furthermore, the higher temperature in the northern part of the study area shows an intensive impact on advancing the SOS. Though global warming is widely believed to advance SOS, the current study reveals that later SOS could be resulted from higher AAT in part of the eastern area close to the Atlantic Ocean (Figure 5a), which is in agreement with the finding that continued warming strengthens the effect of later fulfillment of chilling requirements and attenuates or even reverses the advancing trend in spring phenology [51]. Although warm springs led to an advance of the growing season, warm conditions in winter caused a delay of the spring phases [51].

The impact from annual PREC showed a more complex impact on the phenological events. Vegetation growth can be updated by the availability of soil water, which is not only attributed to precipitation but also to soil properties and landforms. The current study reveals that annual precipitation correlates to the phenological events, but the spatial patterns look rather heterogeneous, which could be resulted from the multiple factors contributing to water availability in soil that is essential to vegetation growth. Though precipitation is an important factor to soil water availability, we suggest that the impact from precipitation should be analyzed with combination of landforms and soil properties. Furthermore, nearby oceans influence regional climate conditions in several ways, including absorbing solar radiation and releasing heat, and/or updating aerosols that influence cloud covers [52], which may interact with precipitation to exert a joint effect on vegetation phenology.

The vegetation growing season length (GL) seems to be closely related to the adjacency to the Atlantic Ocean. The impact from the ocean favors GL, meaning the more closeness to the ocean, the earlier the SOS and later the EOS, and thus a longer growing season length for the vegetation in the forests. This finding supports the previous study that seas or oceans tend to minimize the climate variations [53], meaning the closeness to oceans would favor vegetation growth, and thus advance vegetation green-up and postpone vegetation dormancy.

In addition, we expect that other natural conditions may show impacts on the phenological events. Furthermore, environmental factors aggregated in different seasons (e.g., early growth season and growth season) rather than at an annual time step may present varied impact patterns. To highlight the different impact from the environmental variables, the correlation coefficients between the target (SOS, EOS, or GL) and the updated versions of the explanatory variables aggregated during January-April and May-August were computed and compared (Figure 7). In general, the temperature variables (AAT, MINT, and MAXT) show closer correlation to the phenological events under all the time spans, while the variables related to water availability including soil moisture (SM), water pressure deficit (VPD), and precipitation (PREC) show much less correlation to the phenological events. Surface shortwave radiation (SSR), which is an important factor to vegetation growth, also shows weak correlation to the target variables. Moreover, the environmental factors present a consistent impact on SOS and EOS, given the correlation coefficients with temperature variables (AAT, MINT, and MAXT) are all positive while such a pattern is reversed for GL at all the examined time spans. The findings probably suggest that temperature is the most important factor affecting the phenological events, either positively or negatively, while the other factors show less influence on the phenology of the deciduous forests in the study area.

5.3. Implications of the Shifts in the Phenological Events to Global Carbon Cycling

Global warming due to excessive carbon emissions has been reported and become a severe topic among scientists and governments in many countries of the world [42]. A continuous climbing of temperature seems to have caused and will continue to induce critical consequences to the safety of the global environments. The changes in vegetation phenology over a large area can not only indicate local climate changes but also send feedback to global carbon cycling [16]. Recent development in the remote sensing technology has made it possible to model large-scale spatio-temporal vegetation phenology, thanks to the availability of historically archived remote sensing datasets. When combined with climate datasets, it is more useful to understand the feedback of vegetation phenological events to global carbon cycling.

The current study confirms that the phenological events of the deciduous forests in North America have been greatly updated in the last two decades. Particularly, vegetation green-up dates have been advanced and vegetation dormancy postponed, which has prolonged the plant growing season. One benefit from a longer vegetation growth period is that it may allow forests to photosynthesize more carbon dioxide from the atmosphere, which partially neutralizes the excessively emitted carbon from human activities. However, the net effect from the extended growing season may vary from location to location, with consideration of the heterotrophic respiration from soil [54]. While extended growing season is believed to remove more atmospheric carbon in well-maintained forests, previous studies also pointed out that the soil respiration might overtake the increased amount of carbon absorption from photosynthesis in sparsely vegetation covered forests. For example, forests can emit more carbon than they absorb due to climate change; in southeastern Brazil, the dry and warm seasonal forests combined with gross deforestation were found to progressively absorb less carbon while releasing more over time [55]. From this perspective, the extended growing season length (GL) because of the changes in the SOS and EOS may result in two opposite results, either improved carbon sink or carbon source, depending on the vegetation cover and other environmental factors. However, forest management can help to realize a positive response from the forests, i.e., sequestering more carbon from the atmosphere in the context of the GL changes [56,57]. Therefore, we suggest special attention in terms of forest managements such as forest restoration and recovery programs be given to the areas that have presented a prolonged growing season period over the years (Figure 3e).

6. Conclusions

Forests play a vital role in global carbon cycling as they are the largest carbon pool and sink the majority of the carbon dioxide from the atmosphere [58]. The climate changes may alter the patterns of vegetation phenology and thus update the capacity of carbon sequestration from forests [59]. Understanding the spatio-temporal phenological patterns in forests in the context of global climate changes can provide valuable information for making decisions on forest management targeting improved carbon sequestration. The current study applies remotely sensed imagery to map the phenological events of the deciduous forests in North America and analyzes their spatio-temporal patterns in the last two decades, including the start of green-up season (SOS), end of growing season (EOS), and vegetation growing season length (GL). This study uses threshold algorithm to extract the phenological events by the Savitzky-Golay filtered NDVI, which are consistent with the results from the long-term flux monitoring sites, suggesting that the phenological events could be accurately extracted from remote sensed imagery when combined with appropriate phenology extraction methods. Some critical findings from the current work are summarized as follows,

• The start of growing season (SOS) or vegetation green-up dates represented by day of the year (DOY) of the deciduous forests is closely related to the latitude belts. The SOS increases along the belts from south to north, meaning vegetation green-up date becomes later in the more northern latitudes. Forests showing later end of growing season (EOS) and longer vegetation growth period are found to cluster in the north-eastern and southern parts of the study area.
• The SOS in the deciduous forests generally shows a significant advancing trend and EOS tends to delay to later DOY over the last two decades, resulting in substantial extension of the vegetation growing length (GL), which is likely to sink more carbon as a result of an adaptive response to the continuously rising concentration of carbon dioxide in the atmosphere.
• Vegetation phenology in the forests is closely correlated to climate variations in the study area. Higher annual average temperature (AAT) can lead to an earlier SOS. On the contrary, vegetation will wither earlier in lower latitude belts showing higher AAT, while delayed EOS is mainly observed in areas at higher latitude. The impact of precipitation on the phenological events shows as more spatially heterogeneous than that from AAT.

Though accumulated annual precipitation and annual average temperature are important climatic variables, temporal variation, e.g., monthly difference, in precipitation or temperature could discover a varied pattern in the relationships between climate fluctuation and vegetation phenology, which will be further examined. In addition to the climate factors such as AAT and precipitation, vegetation phenology in forests could be jointly affected by many other factors such as soil properties, landforms, and the closeness to oceans. The current study applies the Pearson’s correlation to examine the impact from each of the climate factors. In the future, advanced analytical models, such as partial least squares structural equation modeling [60], that incorporate all the above important factors in more detail in temporal granularity may be more effective to reveal each direct or their interaction impact on the vegetation phenological events (SOS, EOS, and/or GL) in the deciduous forests.