Abstract
The timing of fall foliage coloration, especially peak coloration, is of great importance to the climate change research community as it has implications for carbon storage in forests. However, its long-term variation and response to climate change are poorly understood. To address this issue, we examined the long-term trends and breakpoints in satellite derived peak coloration onset from 1982 to 2014 using an innovative approach that combines Singular Spectrum Analysis (SSA) with Breaks for Additive Seasonal and Trend (BFAST). The peak coloration trend was then evaluated using both field foliage coloration observations and flux tower measurements. Finally, interannual changes in peak coloration onset were correlated with temperature and precipitation variation. Results showed that temporal trends in satellite-derived peak coloration onset were comparable with both field observations and flux tower measurements of gross primary productivity. Specifically, a breakpoint in long-term peak coloration onset was detected in 25% of pixels which were mainly distributed at latitudes north of 37°N. The breakpoint tended to occur between 1998 and 2004. Peak coloration onset was delayed before the breakpoint while it was transformed to an early trend after the breakpoint in nearly all pixels. The remaining 75% of pixels exhibited monotonic trends, 35% of which revealed a late trend and 40% an early trend. The results indicate that the onset of peak coloration experienced a late trend during the 1980s and 1990s in most deciduous and mixed forests. However, the trend was reversed during the most recent decade when the timing of peak coloration became earlier. The onset of peak coloration was significantly correlated with late summer and autumn temperature in 55.5% of pixels from 1982 to 2014. This pattern of temperature impacts was also verified using field observations and flux tower measurements. In the remaining 44.5% of pixels, 12.2% of pixels showed significantly positive correlation between the onset of peak coloration and cumulative precipitation during late summer and autumn period from 1982 to 2014. Our findings can improve understanding of the impact of changes in autumn phenology on carbon uptake in forests, which in turn facilitate more reliable measures of carbon dynamics in vegetation–climate interactions models.
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Introduction
Vegetation phenology plays a critical role in characterizing land-surface carbon fluxes, parameterizing numerical weather prediction models (Gutman and Ignatov 1998) and land surface process models (Thornes 1985, Zhang et al 2002), and monitoring ecosystem dynamics and climate impacts (Menzel et al 2006, Zhang et al 2014, Liu et al 2016). Recent warming has been reported to lead to earlier spring and later autumn events (Myneni et al 1997, Keenan et al 2014), and increased vegetation activity (Graven et al 2013). Although early greening trends and lengthening of the growing season are expected to increase carbon sequestration and extend the net carbon uptake period in deciduous forests (Black et al 2000, Churkina et al 2005, Dragoni et al 2011), autumn warming suggests net carbon dioxide losses may occur in northern ecosystems (Piao et al 2008).
Spring phenology, in particular, has received much attention from the global change community during the past several decades because it is one of the most widely observed and recorded seasonal events and has proven to be an important bio-indicator of climate change (Beaubien and Freeland 2000, Badeck et al 2004, Richardson et al 2009). However, less effort has been dedicated to studying autumn phenology due to its complex environmental drivers (Cleland et al 2007, Yang et al 2012). Nonetheless, quantifying changes in autumn phenology is important as it alters the reproductive capacity of individual species, regulates growing season duration, and affects net productivity and carbon uptake of ecosystems (Garonna et al 2014, Gallinat et al 2015, Liu et al 2015).
Previous investigations of autumn phenophases mainly focused on intra- or inter- annual variations in senescence and dormancy onset (Jeong et al 2011, Zhu et al 2012, Yang et al 2015, Gill et al 2015). In fact, autumn phenology, ranging from the date when leaves first change color to the date when 100% abscission occurs, is considered to be a multi-event progression (Denny et al 2014). The United States Weather Chanel (www.weather.com/) and field foliage network (www.foliagenetwork.net/) define the fall foliage color status as little/no change, low coloration, moderate coloration, peak coloration and post-peak coloration. Among all fall foliage stages, peak coloration has the broadest range of applications for the tourism industry related to the timing of 'leaf peeping' trips (Ge et al 2013) and significant implications for climate change and ecosystem management activities. Various fall foliage stages at a large scale can be successfully monitored from satellite-derived temporally-normalized brownness indices (Zhang and Goldberg 2011, Zhang et al 2012). However, examination of the temporal trend in peak coloration and its response to climate change at a large scale has not yet been undertaken.
The temporal trend in phenological events has been frequently detected using a simple linear regression (SLR) approach to monitor recent climate change (Menzel and Fabian 1999, Schwartz and Reiter 2000, Garonna et al 2014). The SLR method provides a slope representing the average rate of change in phenological events over a study period, but this approach has a number of major limitations when searching for signals in short and noisy phenological times series (Sparks and Menzel 2002, Sparks and Tryjanowski 2005). Specifically, the SLR method requires the time series to have a comparatively linear trend and it is very sensitive to outliers and boundary values (Schlittgen and Streitberg 1999). Moreover, the SLR slope (magnitude) is largely influenced by the length of the time series and its start and end dates, especially for highly variable phenological time series over a few decades. In order to overcome some of these shortcomings, Keatley and Hudson (2004) introduced the use of singular spectrum analysis (SSA) into the examination of long-term phenological time series. SSA has no specific distributional assumptions. In addition, it uses a data-adaptive basis set, instead of fitting an assumed model to the available series, to overcome the problems of finite sample length and noisiness of sampled time series (Hudson and Keatley 2010).
Analyzing long-term phenological time series reveals temporal variations associated with the study period. This is evident in North America where trends in spring green up differ dramatically between different time periods, such as, 1982–2005 (Zhang et al 2007), 1982–2006 (White et al 2009), 1982–2003 (Reed 2006), and 1982–2008 (Jeong et al 2011). The difference is likely associated with the research approaches in these studies but the length of the time series, start and end dates, and extreme values are also critical in determining whether the trends can be detected at a significant level (Sparks and Menzel 2002). Thus, it is necessary to incorporate mathematical methods capable of automatically identifying the change point in a long-term phenological time series. Such change points have been detected using several approaches, such as Breaks for Additive Seasonal and Trend (BFAST) (Verbesselt et al 2010), nonparametric methods (Itoh and Kurths 2010), and Bayesian-based approaches (Dose and Menzel 2004, Mohammad-Djafari and Féron 2006).
This study aimed to investigate the spatio-temporal change in the onset of peak coloration in both deciduous broadleaf and mixed forests in the central and eastern United States from 1982 to 2014 and its responses to climate change. To reach this goal, a time series of daily EVI2 (two band enhanced vegetation index) was derived from the Advanced Very High Resolution Radiometer (AVHRR, 1982–1999) and Moderate Resolution Imaging Spectroradiometer (MODIS, 2000–2014) records at a spatial resolution of 0.05°. A Hybrid Piecewise Logistic Model (HPLM) was then used to reconstruct the EVI2 time series. The reconstructed EVI2 data were converted to a temporally-normalized brownness index to determine the onset of peak coloration. Subsequently, the temporal trends and change points in the time series of long-term peak coloration onset were both identified at a pixel level by combining SSA with BFAST methods. In addition, the trends were evaluated using field observations including leaf coloration and gross primary productivity data. Finally, the response of peak coloration onset to interannual variation in temperature and precipitation was examined at each pixel and field observation site.
Materials and methods
Phenological detection from satellite data
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).
Figure 1 Spatial distribution of the flux sites (solid circles) in deciduous broadleaf and mixed forests (green) in the central and eastern United States. The solid red circles represent the sites with gap-filled GPP datasets over 10 yr, which were used in this study.
Download figure:
Standard image High-resolution imageHPLM 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.
Field observations of fall foliage coloration at Harvard Forest
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).
Phenological detections from GPP data
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).
Figure 2 An example of raw GPP data with a reconstructed curve from Park Falls in 2001 and the diagram of the dynamic threshold approach for extracting the onset of peak coloration (PC).
Download figure:
Standard image High-resolution imageDaily temperature and precipitation
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.
Trend analysis for phenological time series
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.
Investigation of peak coloration response to temperature and precipitation
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:
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.
Results
Climatology of 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)).
Figure 3 Spatial distribution of peak coloration onset in deciduous and mixed forests in the central and eastern United States from 1982 to 2014. (a) Mean onset of peak coloration (DOY) and (b) standard deviation. The solid line is the state boundary in the United States.
Download figure:
Standard image High-resolution imageEvaluation of trends in satellite-derived peak coloration onset at three observation sites
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.
Figure 4 Long- and short-term trends in peak coloration onset derived from satellite (a) and (b), field observation (c) and (d), and GPP data (e) and (f) at Harvard Forest from 1991 to 2014. The straight dashed line represents linear regression line (a), (c) and (e). The vertical dashed line represents the breakpoint in the time series (b), (d) and (f).
Download figure:
Standard image High-resolution imageAn 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−1 (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.
Figure 5 Long- and short-term trends in peak coloration onset derived from satellite data (a) and (b) and GPP data (c) and (d) at Park Falls from 1997 to 2014 using linear regression (a) and (c) and SSA (b) and (d). The straight dashed line represents linear regression line (a) and (c).
Download figure:
Standard image High-resolution imageThe 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)).
Figure 6 Long-term and short-term trends in peak coloration onset derived from satellite (a) and (b) and GPP data (c) and (d) in Morgan Monroe State forest from 1999 to 2014 using linear regression (a) and (c) and SSA (b) and (d) methods. The straight dashed line represents linear regression line (a) and (c). The vertical dashed line represents the breakpoint in the time series (b) and (d).
Download figure:
Standard image High-resolution imageSpatial pattern in peak coloration trends
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)).
Figure 7 Spatial pattern of long-term trends in peak coloration onset (DOY) from 1982 to 2014 in deciduous and mixed forests in the central and eastern United States. (a) The year of breakpoint occurrence; (b) trend (d y−1) in pixels before the breakpoint (first time period, from 1982 to the year of breakpoint occurrence); (c) trend (d y−1) in pixels after the breakpoint (second time period, from the year of breakpoint occurrence to 2014); and (d) trend (d y−1) in pixels without a breakpoint. The dark grey color in (a), (b) and (c) represents the pixels without a breakpoint, but in (d) represents the pixels with a breakpoint.
Download figure:
Standard image High-resolution imageCorrelation of peak coloration onset with temperature and precipitation at three observation sites
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.
Figure 8 Correlation coefficients between peak coloration and mean temperature during the different PLs at Harvard Forest from 1991 to 2014 (a), Park Falls from 1997 to 2014 (d) and Morgan Monroe State Forest from 1999 to 2014 (g); long-term trend in mean temperature within the optimal PL using linear regression (b), (e), (h) and SSA (c), (f), (i) at Harvard Forest, Park Falls and Morgan Monroe State Forest, respectively.
Download figure:
Standard image High-resolution imageCorrelation of peak coloration onset with temperature and precipitation across the central and eastern United States
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.
Figure 9 Spatial pattern of correlation between peak coloration onset and temperature or precipitation. (a) Optimal period length of temperature, (b) correlation coefficient between peak coloration onset and mean temperature within the optimal PL (P < 0.1), (c) optimal period length of cumulative precipitation, and (d) correlation coefficient between peak coloration onset and cumulative precipitation within the optimal PL (P < 0.1). The dark grey color indicates the deciduous and mixed forests with non-significant correlation between the onset of peak coloration and optimal PL temperature or precipitation.
Download figure:
Standard image High-resolution image
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).
Table 1. Pixel proportions (%) with significant and non-significant correlations between the onset of peak coloration and mean temperature within the optimal PL, which were then divided by the correlations between the onset of peak coloration and cumulative precipitation within the optimal PL. Note no pixels with significantly negative correlations between the onset of peak coloration and mean temperature. Sig_Positive—significantly positive correlation, Sig_Negative—significantly negative correlation, and Sig_Non—non-significant correlation.
| Temperature | Sig_Positive | Sig_Non | ||||
| 55.5% | 44.5% | |||||
| Precipitation | Sig_Negative | Sig_Positive | Sig_Non | Sig_Negative | Sig_Positive | Sig_Non |
| 18.8% | 5.4% | 31.3% | 5.8% | 12.2% | 26.5% | |
Discussion and Conclusions
Statistically significant temporal trends in long-term phenological time series are not always detected using a linear regression model because phenological changes may not exhibit a monotonic trend (Menzel et al 2006, Shen et al 2014, Zhou et al 2014). However, unlike linear regression SSA does not impose linearity on the time series and has been successfully used to describe the complex variation in spring phenology of three perennial plant species from 1983 to 2002 (Keatley and Hudson 2004). The current study demonstrated the advantage of using SSA in the analyses of a long-term time series of peak coloration onset derived from satellite data. For example, a simple linear regression found no significant trends in the onset of peak coloration observed from either satellite, field, or flux tower measurements at Harvard Forest from 1991 to 2014 but SSA revealed a delay in peak coloration onset from 1991 to 2002 and an early trend from 2003 to 2014. Although SSA produces a set of components, the first component of the time series of peak coloration onset described variation in the long-term trend because it accounted for more than 99.0% of signal (Vautard et al 1992). The trends identified by SSA were considered to be 'significant' because the first component was usually above the noise floor (Shun and Duffy 1999, D'Odorico et al 2002), however, further work is required to apply statistical significance to SSA-based trends (Elsner and Tsonis 1996).
Breakpoints were identified from the first component of peak coloration time series derived from SSA because it contained very high signal-to-noise ratio (Verbesselt et al 2010). In particular, we determined breakpoints caused by reversed slope directions for each pixel using BFAST although breakpoints could also be a result of an abrupt change in amplitude. Abrupt changes in slope amplitude may occur with a change in satellite sensor (De Beurs and Henebry 2005) such as inconsistencies between AVHHR on different satellites or between AVHRR and MODIS. However, the reversed slope directions are unlikely caused by sensor artifacts. Given that the break point mainly shifts between 1998 and 2004 and that observations from field sites and satellite tend to align, the break points identified in this study were unlikely associated with the transition from AVHRR to MODIS. This implies that the breakpoints derived in the current study were mainly associated with climatic variation. It should be noted, however, that the SSA-BFAST method in detecting break points requires the choice of a length of window (or segment) based on the structural changes of the time series (Keatley and Hudson 2004, Verbesselt et al 2010). Thus, the change of window length could have potential influences on the results.
This study is the first to investigate and validate temporal trends in satellite derived fall foliage coloration. Validation of satellite-based peak coloration is important for accurate determination of the onset of the autumn season but direct validation is challenging and impractical for a coarse resolution pixel. This is because the samples of field observations that are spatially compatible with satellite coarse pixels are very limited although landscape phenology could be upscaled from plot-level observations (Liang et al 2011). However, long-term trends should be consistent if interannual variation between a small plot and a large pixel is similarly controlled by climate change. Therefore, satellite-derived long-term trends in this study were validated using trends in field foliage color observations and GPP measurements over a period of more than 10 years. The validation revealed that the temporal trends in peak coloration onset from the three datasets were consistent whether the breakpoint existed (such as Harvard Forest and Morgan Monroe State Forest) or not (such as Park Falls). This suggests that the satellite measurements could effectively monitor the temporal trend in peak coloration onset of deciduous and mixed forests.
The trend analysis used here demonstrated that a breakpoint was detected in 25% of the pixels in which the onset of peak coloration was delayed over the first time period and advanced over the second time period. These breakpoints mainly occurred during 1998–2004 and distributed at latitudes north of 37°N. In the remaining 75% of pixels, peak coloration onset was delayed in 35% and advanced in 40% of pixels from 1982 to 2014. These results suggested that the onset of peak coloration in most deciduous and mixed forests across the northeastern United States switched from a delayed trend in the 1980s and 1990s to an early trend during the past decade. Compared to previous work simply showing delayed autumn phenology across a range of ecosystems over the past several decades (Jeong et al 2011, Dragoni and Rahman 2012, Gill et al 2015), our study revealed the dynamic trends in autumn phenology for all pixels using SSA. Earlier autumn phenology was found to also contribute to a shorter growing season, which may shorten the period of net carbon dioxide uptake and thus decrease carbon storage (Churkina et al 2005).
Our results showed strong evidence of a control of temperature within an optimal PL on interannual variation in peak coloration onset of deciduous and mixed forests in the central and eastern United States. Examinations of the validation sites revealed that the long-term trend in peak coloration onset basically tracked the trend in optimal PL temperature although the pattern varied greatly with variation in local climate. The results across the central and eastern United States demonstrated that peak coloration onset was significantly correlated with the mean temperature within an optimal period length in about 55.5% of pixels. Spatial patterns indicated that the optimal period length was about two months shorter in northern than in southern regions. This latitudinal gradient was comparable to the variation in the onset of peak coloration. This coincident spatial variation between peak coloration onset and optimal period length likely suggests that summer temperature may influence autumn foliage coloration. The spatial variation in peak coloration onset and optimal period length also implies that temperature impacts vary geographically instead of consistently across the region, which is likely controlled by the geographic adaption relationships between vegetation phenology and climate influences through both contemporary and historical variability (Liang et al 2016). The plant optimal adaption or temperature memory may determine to a large extent the geographic distributions and vegetation seasonality (Rohde and Junttila 2008). This has been demonstrated in a recent study that spatial pattern of long-term spring phenology across the United States was associated with vegetation climate adaptation (Liang et al 2016). However, mapping vegetation climate adaptation patterns at a large scale is challenging (Liang et al 2016) and the mechanism of plant memory to dynamic environments is still poorly understood (Crisp et al 2016). As a result, further research is needed to investigate the impacts of temperature memory on the autumn phenology and its spatial pattern.
Precipitation could be another important factor influencing peak coloration onset (Leuzinger et al 2005, Shen et al 2015), but the mechanism of effects of precipitation on autumn phenology is much more complex than temperature. The correlation between peak coloration onset and cumulative precipitation within the optimal PL was significantly positive in 12.2% of pixels where non-significant correlations between the onset of peak coloration and mean temperature within the optimal PL were found, which indicated that less precipitation during late summer and autumn period could advance the occurrence of peak coloration if temperature controls were weak. This is because more precipitation, especially in summer, should promote vegetation growth, which, in turn, could delay the timing of peak coloration onset (Shen et al 2015). The combination of temperature and precipitation could explain the interannual variations in the onset of peak coloration in 67.6% (55.5% plus 12.2%) of pixels across the central and eastern United States. However, 24.6% of pixels showed significantly negative correlations between the onset of peak coloration and cumulative precipitation within the optimal PL, which did not necessarily have a particular biophysical meaning. This is because the onset of peak coloration in most of them (18.8%) was significantly controlled by temperature. In other words, the influence of precipitation was likely overridden by temperature.
In conclusion, satellite data over the past 33 years revealed a delayed trend in the onset of peak coloration during the 1980s and 1990s. However, this trend was reversed during the most recent decade when the onset of peak coloration became earlier in most deciduous and mixed forests of the central and eastern United States. The change in trend was strongly associated with variation in late summer and autumn temperature in 55.5% of pixels. In 12.2% of the remaining 44.5% of pixel, precipitation could well explain the interannual variations in the onset of peak coloration. These results are important for determining the impact of future changes in the onset of peak coloration on carbon uptake in forests across the US and elsewhere.
Acknowledgments
This work was supported by NASA contracts NNX15AB96A and NNX14AJ32G, and NOAA contract JPSS_PGRR2_14.