Detecting spatiotemporal changes of peak foliage coloration in deciduous and mixedforests across the Central and Eastern United States

The long-term daily EVI2 time series over the last 33 yr (1982–2014) was used to retrieve the onset of peak coloration in deciduous and mixed forests across the central and eastern United States (figure 1). Specifically, EVI2 was calculated from the daily red and near-infrared reflectance in the AVHRR long term data record (LTDR, 1982–1999) and the MODIS climate modeling grid (CMG, 2000–2014) data at a spatial resolution of 0.05° (∽5 km). This dataset was obtained from the University of Arizona (http://vip.arizona.edu/viplab_data_explorer). EVI2 has advantages over NDVI (normalized difference vegetation index) in quantifying vegetation activity but also remains functionally equivalent to EVI (enhanced vegetation index) (Huete et al 2002, Jiang et al 2008). AVHRR LTDR quality was improved by preprocessing radiometric corrections, viewing and illumination adjustment, and cloud screen and water-vapor correction (Pedelty et al 2007). Furthermore, AVHRR LTDR was bridged to MODIS CMG using a set of linear polynomials (Tsend-Ayush et al 2012).

Zoom In Zoom Out Reset image size

HPLM was used to reconstruct the EVI2 temporal trajectory, which has been demonstrated to effectively minimize noise, such as the presence of snow cover and cloud, remove irregular values, and fill in missing observations in the EVI2 time series (Zhang 2015). The reconstructed EVI2 time series from the mid maturity phase to mid dormancy phase was then used to calculate the temporally-normalized brownness index during the senescence phase at the pixel level (Zhang and Goldberg 2011). The temporally-normalized brownness index is capable of measuring the relative change in colored foliage within a pixel. It is independent of surface background, vegetation abundance, and species composition. In particular, the onset of peak coloration was defined as the date at which the temporally-normalized brownness index was equal to 0.6, which presented the maximal colored foliage cover in tree canopies (Zhang and Goldberg 2011).

The onset of peak coloration in the deciduous and mixed forests was extracted by overlaying the extracted peak coloration map on the MODIS Land Cover Type CMG product (MCD12C1) at a 0.05° spatial resolution. The climatology of peak coloration onset was further calculated using average and standard deviation from 1982 to 2014.

To evaluate the temporal trend in satellite-derived peak coloration, field observations of fall foliage coloration were collected from Harvard Forest from 1991 to 2014 (http://harvardforest.fas.harvard.edu:8080/exist/apps/datasets/showData.html?id=hf003).This site is located at 42°32^'N and 72°11^'W in central Massachusetts and is the only one with a long-term record of woody plant fall phenology in the United States. The mean annual temperature is 6.62 °C and mean annual precipitation is 1071 mm. The fall foliage progress was observed for 33 dominant species representing both canopy and understory individuals. The observations included the percentage of leaves that had changed color on a given plant and the percentage of leaves that had fallen at 3- to 7 d intervals during the senescence phase. The species based observations could be accurately aggregated by weighting their abundance to be compatible with remotely sensed phenology (Liang et al 2011, Liu et al 2015), however, the limitation of an unknown percent area of each species in this site makes it impossible to obtain area-weighted foliage coloration data. Thus, to match field measurements to satellite pixels, we calculated the average fraction of foliage coloration at each observation date. The temporal observations of foliage coloration were subsequently fitted using a logistic model for each year. The start date of peak coloration was estimated when the percentage of colored foliage reached 60% (Zhang and Goldberg 2011).

We further collected flux-tower measurements of gross primary productivity (GPP) as an alternative dataset to further evaluate the trend in satellite-based peak coloration onset. Across the central and eastern United States, there are a total of 21 flux tower sites with deciduous broadleaf and mixed forests (figure 1). However, there are only three sites with gap-filled GPP data longer than 10 yr, which are Harvard Forest in Massachusetts (US-Ha1; 1991–2014), Morgan Monroe State Forest in Indiana (US-MMS; 1999–2014), and Park Falls in Wisconsin (US-PFa; 1997–2014). At Park Falls, the mean annual temperature is 4.33 °C and mean annual precipitation is 823 mm. At Morgan Monroe State Forest, the mean annual temperature is 10.85 °C and mean annual precipitation is 1032 mm. The level2 flux data in these three sites were downloaded from the American FLUXNET (http://ameriflux.ornl.gov).

The level2 data were recorded hourly, so that we first calculated daily mean GPP using the arithmetic average. Due to the noise inherent in original GPP time series, we smoothed the daily GPP time series using the HPLM method (Zhang 2015). An example of smoothing daily GPP time series is presented in figure 2. Based on daily smoothed time series, a dynamic threshold method could be applied to detect fall phenological events (Richardson et al 2010). In this study, we calculated the GPP-based onset date of peak foliage coloration according to the dynamic threshold that was the same as the threshold for the temporally-normalized brownness index (figure 2).

Zoom In Zoom Out Reset image size

We obtained 3-hourly land surface temperature and precipitation data, from 1982 to 2014 at a spatial resolution of 32 km for the United States from the NCEP (National Centers for Environmental Prediction) North America Regional Reanalysis (NARR) (Mesinger et al 2006). This dataset was produced using a fixed assimilation/forecast model and is the most accurate and consistent long-term dataset that covers the entire North American continent (Mesinger et al 2006). The daily mean temperature was generated by averaging the 3-hourly land surface temperature values. The daily total precipitation was generated by summing the 3-hourly precipitation values.

Field measurements of temperature and precipitation were obtained for the field sites. The daily air temperature and precipitation at Harvard Forest was downloaded from ClimDB/HydroDB Data Retrieval Program (http://climhy.lternet.edu/plot.pl). The air temperature and precipitation at Morgan Monroe State Forest and Park Falls were included in the level2 flux dataset. The daily mean temperature and daily total precipitation at these two flux tower sites were computed by averaging hourly values and summing hourly values, respectively.

The satellite-based onset of peak coloration from 1982 to 2014 in the central and eastern United States was examined using SSA in the R package ('Rssa') to identify temporal trends. The SSA method is a time series analysis method which is capable of decomposing and forecasting time series. The original time series is decomposed by singular value decomposition (SVD) as a trend and oscillatory components that could be related to seasonality and noise (Golyandina and Korobeynikov 2014). The principal components that were considered to describe the trend were determined using the scree plot which is a plot of the eigenvalues of a correlation matrix in descending order of their magnitude (Vautard and Ghil 1989). Given that SSA requires that the datasets are complete, we only analyzed phenological time series without missing values. The window length in SSA was set to 10 years due to the climate variability in 1980s, 1990s, and 2000s.

Breakpoints in the first reconstructed series from SSA were examined using the BFAST method in the R package ('BFAST'). This method is capable of capturing long-term phenological changes, abrupt changes and gradual trends in satellite image time series (Verbesselt et al 2010). Specifically, the ordinary least squares residuals-based moving sum (MOSUM) test was first used to test whether one or more breakpoints occurred in the time series. If the MOSUM test indicated significant change (p < 0.05), the number of breakpoints was determined by the Bayesian Information Criterion (BIC), and the date and associated confidence interval for each breakpoint were estimated (Zeileis et al 2002, Bai and Perron 2003, Zeileis 2005). In this study, the breakpoint was defined as the year that divided the phenological time series into two stages of an early and a late trend. Although the BFAST was able to detect all breakpoints that were caused by changes in either amplitude or direction of the slope in the long-term phenological time series, we only selected the breakpoint in which the slope direction changed. For pixels in which a breakpoint was detected, the entire time period was separated into the first time period and the second time period.

Similarly, temporal trends and breakpoints were examined for the long-term field-observed phenological time series at Harvard Forest and GPP-based peak coloration time series at Harvard Forest, Morgan Monroe State Forest, and Park Falls flux tower sites using the SSA-BFAST method.

The timing of phenological events have been found to be mainly influenced by temperature within a specified period length (PL, days) prior to the phenological event occurrence (Matsumoto et al 2003, Chen and Xu 2012). This specified period can be quantified using a critical climate period analysis (Craine et al 2012). To determine the optimal PL (critical climate period), we first calculated the correlation coefficients between the onset of peak coloration in a time series and a set of mean temperatures calculated with different PL values. The PL was defined as follows:

$$\begin{equation} {\rm PL} = {\rm bPL} + {\rm mPL} \end{equation}$$

where bPL is the basic PL that refers to the time period (days) between the earliest and latest date in the time series of peak coloration onset in a pixel, and mPL is a lag changing from 1 to 90 d at a one-day interval.

For a given mPL, the mean temperature during the PL was calculated for each year. A correlation was further established from the time series (from 1982 to 2014) of peak coloration onset and the mean temperature for a given pixel. With the variation in mPL (1–90 d) and the corresponding PL, 90 correlation coefficients were generated. Subsequently, we obtained the optimal PL by identifying the largest correlation coefficient for a given pixel.

This approach was also applied to calculate the optimal PL between the onset of peak coloration and cumulative precipitation in order to investigate the effects of precipitation. The optimal PL of cumulative precipitation differs from that of mean temperature because these two factors influence the onset of peak coloration with different mechanisms.

Similarly, we also determined the optimal PL temperature and precipitation for the field observations at Harvard Forest, Morgan Monroe State Forest, and Park Falls. Furthermore, we examined the temporal trend and breakpoint in optimal PL mean temperature and cumulative precipitation time series using both simple linear regression and SSA-BFAST to explain the long- and short-term trends in peak coloration onset.

Figure 3 shows the spatial pattern of mean peak coloration onset from 1982 to 2014 together with the standard deviation for deciduous and mixed forests across the central and eastern United States. In general, the timing of peak coloration onset in day of year (DOY) was gradually delayed with decreasing latitude, ranging from DOY 260 in northern Minnesota, Wisconsin and Maine to DOY 320 in the central region of Louisiana, Alabama and Georgia (figure 3(a)). The spatial shift took about two months from northern to southern areas. The standard deviation of peak coloration dates was generally less than 10 d in the northern region, whereas it was up to 20 d in the southern area (figure 3(b)).

Zoom In Zoom Out Reset image size

While landscape phenology can be scaled up from field observations of individual species to effectively evaluate satellite-derived phenology (Liang et al 2011, Liu et al 2015), site-specific comparison of long-term trends provides an alternative way to evaluate the reliability of the satellite-derived trends. Figure 4 displays the interannual variation in the timing of peak coloration onset derived from satellite data, field foliage color observations and GPP data at Harvard Forest. Linear least-squares regression, which has been commonly used to determine temporal trends in phenological events, showed no significant trend in any of the three datasets (figures 4(a), (c) and (e)). In contrast, SSA produced the first reconstructed series of peak coloration onset, which accounted for more than 99% of the variation and closely represented the dynamic trend. The trend indicated that satellite-based peak coloration was delayed by up to 3 d from 1991 to 2001 and advanced by up to 5 d from 2002 to 2014 (figure 4(b)). Similarly, the field-based peak coloration was also delayed by up to 3 d from 1991 to 2002 and advanced by up to 4 d from 2003 to 2014 (figure 4(d)). The GPP-based peak coloration showed a late trend from 1991 to 2001 an early trend from 2002 to 2006 and little change from 2007 to 2014 (figure 4(f)). The breakpoints calculated using the BFAST method were found to be located in 2001 or 2002 for all three datasets. Evidently, the trend in satellite-based peak coloration onset was consistent with field observations and GPP-based measurements during two periods, 1991–2001 and 2002–2014, at Harvard Forest.

Zoom In Zoom Out Reset image size

An early trend in the timing of peak coloration was revealed from both satellite observations and GPP measurements at Park Falls with no distinguishable breakpoints (figure 5). The satellite-derived trend became earlier at a rate of 0.5 d yr⁻¹ (P < 0.05) from 1997 to 2014 as revealed by linear regression (figure 5(a)), which was similar to the temporal pattern derived from SSA (figure 5(b)). Similar temporal patterns were found for GPP-based peak coloration onset although the GPP data presented a somewhat larger interannual variation than the satellite data (figure 5(c)). The GPP-based peak coloration onset was advanced by up to 20 d from 1997 to 2014 (figure 5(d)). Overall, the temporal trends in peak coloration onset as determined from satellite and GPP showed close agreement.

Zoom In Zoom Out Reset image size

The time series of peak coloration onset from both satellite and GPP data at Morgan Monroe State forests exhibited no significant trends from 1999 to 2014 using simple linear regression (figures 6(a) and (c)). However, using SSA, we found that satellite-based peak coloration onset was delayed by up to 3 d from 1999 to 2006 and advanced by up to 4 d from 2007 to 2014 (figure 6(b)). This temporal pattern was in close agreement with the trends in GPP-based peak coloration onset from 1999 to 2014 (figure 6(d)).

Zoom In Zoom Out Reset image size

Figure 7 presents the spatial pattern of the long-term trend in satellite-based peak coloration onset in both deciduous and mixed forests using SSA-BFAST. Among all pixels, 25% were determined to have one breakpoint that mainly occurred from 1998 to 2004 although exceptions were found in a few pixels (figure 7(a)). These pixels were mainly distributed at latitudes north of 37°N, including northern Minnesota, Wisconsin and Michigan, western Maine, eastern New Hampshire, and northeast New York (figure 7(a)). In 96% of the 25% pixels with one breakpoint, the onset of peak coloration was delayed during the first time period, but it transitioned to an advance in timing following the break point (figures 7(b) and (c)). 75% of the pixels showed monotonic trends, 35% of which exhibited a late trend and 40% an early trend (figure 7(d)). The pixels with a monotonic late trend were mainly distributed in northern Minnesota, Wisconsin and Michigan, West Virginia and western Louisiana and Arkansas. In contrast, pixels exhibiting a monotonic early trend were mainly distributed in eastern Maine, southeastern New York, Pennsylvania, and western Alabama (figure 7(d)).

Zoom In Zoom Out Reset image size

Figure 8 illustrates the correlation between peak coloration onset and mean temperature within the optimal PL at three field sites. The correlation coefficient indicated that the onset of peak coloration was significantly correlated with mean temperature within the optimal PL (P < 0.05) although the optimal PL was 76 d, 25 d, and 68 d at Harvard Forest, Park Falls and Morgan Monroe State Forest, respectively (figures 8(a), (d) and (g)). Overall, the trends in mean temperature within the optimal PL varied greatly with sites. At Harvard Forest, the mean temperature significantly increased by 0.04 °C per year (P < 0.05) from 1991 to 2014 (figure 8(b)). This pattern could not fully explain the non-significant changes in peak coloration (figure 4(c)). However, SSA revealed that the first reconstructed series of optimal PL mean temperature accounted for more than 98% of the signal. The first reconstructed temperature series increased from 1991 to 2004 then decreased from 2005 to 2009 and slightly increased again from 2010 to 2014 (figure 8(c)). This temporal pattern could better explain the delay in peak coloration from 1991 to 2002 and the advance from 2003 to 2014 (figure 4(d)). At Park Falls, Wisconsin, the mean temperature displayed with a large amount of interannual variation (SD = 1.67 °C). Therefore, no significant decrease was found using linear regression (figure 8(e)). However, SSA revealed that optimal PL temperature clearly decreased from 1997 to 2014 (figure 8(f)), which could explain the advanced peak coloration onset at this site (figure 5). In Morgan Monroe State Forest, simple linear regression revealed no significant temporal trend in mean temperature over the 16 yr of the study (figure 8(h)). In contrast, SSA indicated that the mean temperature increased by 0.5 °C from 1999 to 2006 and decreased by 0.5 °C from 2007 to 2014 (figure 8(i)), which was in close agreement with temporal changes in the onset of peak coloration at this site (figure 6). In addition, there were no significant correlation between the onset of peak coloration and cumulative precipitation within the optimal PL at these three sites, so we did not present these results.

Zoom In Zoom Out Reset image size

Figure 9 presents the spatial pattern of optimal period length and the correlation coefficients of peak coloration onset with mean temperature and cumulative precipitation during the optimal PL in deciduous and mixed forests, respectively. In general, the optimal PL followed a latitudinal gradient, which was longer in southern than in northern regions (figure 9(a)). Specifically, autumn optimal PL was around 30 d in northern Minnesota, around 60 d in Maine, Connecticut, Massachusetts and West Virginia, and around 90 d in Western Arkansas, Louisiana, central Alabama and Georgia. The daily mean temperature within the optimal PL was significantly correlated with peak coloration onset (P < 0.1) in approximately 55.5% of pixels (table 1). The pixels with significant positive correlations were mainly distributed in northern Minnesota and Michigan, Maine, Arkansas and Louisiana, Connecticut, and West Virginia.

Zoom In Zoom Out Reset image size

The optimal PL of cumulative precipitation varied across the central and eastern United States, ranging from 30 d to 150 d, but regular spatial pattern was not revealed (figure 9(c)). The correlation between the onset of peak coloration and cumulative precipitation within the optimal PL was significantly positive in 17.6% pixels and negative in 24.6% pixels (figure 9(d) and table 1). Compared to figures 9(b) and (d) indicates that the positive correlations were mainly distributed in the region where the peak coloration onset was non-significantly correlated with mean temperature, such as northwestern Pennsylvania, central Virginia, western Alabama, Mississippi, eastern South Carolina and Georgia. In contrast, the negative correlation was generally distributed in Minnesota, Michigan, Maine, Arkansas, Louisiana (figure 9(d)), where the peak coloration onset was significantly positively correlated with mean temperature (table 1).