Linking Spaceborne and Ground Observations of Autumn Foliage Senescence in Southern Québec, Canada

Abstract

Autumn senescence progresses over several weeks during which leaves change their colors. The onset of leaf coloring and its progression have environmental and economic consequences, however, very few efforts have been devoted to monitoring regional foliage color change in autumn using remote sensing imagery. This study aimed to monitor the progression of autumn phenology using satellite remote sensing across a region in Southern Québec, Canada, where phenological observations are frequently performed in autumn across a large number of sites, and to evaluate the satellite retrievals against these in-situ observations. We used a temporally-normalized time-series of Normalized Difference Vegetation Index (NDVI) extracted from Moderate Resolution Imaging Spectroradiometer (MODIS) imagery to monitor the different phases of autumn foliage during 2011–2015, and compared the results with ground observations from 38 locations. Since the NDVI time-series is separately normalized per pixel, the outcome is a time-series of foliage coloration status that is independent of the land cover. The results show a significant correlation between the timing of peak autumn coloration to elevation and latitude, but not to longitude, and suggest that temperature is likely a main driver of variation in autumn foliage progression. The interannual coloration phase differences for MODIS retrievals are larger than for ground observations, but most ground site observations correlate significantly with MODIS retrievals. The mean absolute error for the timing of all foliage phases is smaller than the frequency of both ground observation reports and the frequency of the MODIS NDVI time-series, and therefore considered acceptable. Despite this, the observations at four of the ground sites did not correspond well with the MODIS retrievals, and therefore we conclude that further methodological refinements to improve the quality of the time series are required for MODIS spatial monitoring of autumn phenology over Québec to be operationally employed in a reliable manner.

1. Introduction

During the autumn season, as the day light hours shorten, and the temperatures drops, the leaves of many deciduous trees and shrubs change their color from green to various shades of yellow, orange, red and brown. This phenomenon of senescence progresses over several weeks, as the transport of water and nutrients to the leaves is gradually impeded, causing the chlorophyll in the leaves to decrease [1,2]. The high chlorophyll content during the growing season often masks the color of other pigments, but as the chlorophyll degrades, the orange-yellow color of carotenoids becomes dominant [3]. In addition, towards the end of summer, red-purple anthocyanins develops in the leaves of some plants. Accordingly, the color of many leaves during senescence is determined by a combination of these two types of pigments and the brown color of the leaf cell walls to produce a collage of autumn colors [3].

Autumn foliage color changes are of interest for several reasons. The gradual change of the landscape color has significant effects on the earth’s surface energy balance [4]. Senescence, which marks the end of the growing season, has a central role in biogeochemical cycles, and strongly affects both the surface albedo and primary productivity of terrestrial systems [5,6]. Plant phenology, and specifically the timing of autumn coloration, is an effective indicator for climate change [7,8]. It is also found to be closely synchronized with other natural phenomena such as trout spawning events [9]. In addition, the timing of autumn coloration is of major interest to the tourism industry in certain areas. Many tourists are drawn to the varied mosaic of autumn foliage colors [10]. For example, this has been estimated to be a US $400 million tourism industry in parts of north-eastern USA and south-eastern Canada (in 2008) [11]. As such, the timing and progression of the autumn foliage phenomena is of significant interest for environmental, as well as economic, reasons [12].

Vegetation phenology, and more specifically autumn foliage, is monitored in several distinct ways. Traditionally, vegetation phenology was observed from the ground. Frequent ground observations of phenology in some places extends back for many years, although this is labor intensive and observer bias is known to occur [13]. Recently, digital cameras have also been used to observe and monitor vegetation phenology from the ground [14,15]. Ground based cameras have the advantage of providing a continuous and permanent record that can be inspected at any time. In theory, their maintenance is not as labor intensive as conducting ground surveys at a high temporal frequency. On the other hand, they cannot be used to obtain a synoptic perspective, as is the case with satellite remote sensing. In the last few decades, satellite remote sensing has emerged as an important method to monitor the phenology of vegetation over large areas [8]. So far, the majority of remote sensing studies have focused on spring phenology, while autumn phenology has received far less attention [16], although some studies have considered the full length of the growing season [17,18,19].

Remote sensing from both satellite and ground based sensors has been shown to be more reliable for spring compared to autumn phenology when correlated with visual ground observations [14,17]. The reason for this is that the rate of change in vegetation indices tends to be more gradual in the autumn, which makes estimates of phenological metrics more challenging relative to spring [20]. Nevertheless, in large areas where ground observations are sparse, satellite imagery provides a means to track the phenomenon in real-time, not only at discrete locations, but in a spatially continuous manner. The gradual change of autumn senescence calls for special attention to the different phases of its progression. However, most remote sensing studies that monitored the autumn’s onset reduced it to a single end-of-season date [17,18,19], and only a minority of the studies follow the progression of the coloration phenomenon throughout its development [21]. As such, there is significant interest in developing remote sensing approaches to estimate the progression of autumn foliage [22].

Previous remote sensing studies have followed the progression of autumn foliage by analyzing a time series of vegetation indices, namely the Normalized Difference Vegetation Index (NDVI) [8,20,21,22,23,24], the Enhanced Vegetation Index (EVI) [20,22,25], Leaf Area Index (LAI) [26], Normalized Difference Water Index (NDWI) [27] and the Normalized Difference Infrared Index (NDII) [28]. Most commonly, the data used to compose the time-series is derived from the Moderate Resolution Imaging Spectroradiometer (MODIS), Advanced Very High Resolution Radiometer (AVHRR) and Satellite Pour l’Observation de la Terre-Vegetation (SPOT-VGT) [29]. While many different algorithms are used to analyze these time-series, the most prevalent approaches use curve fitting to model the vegetation index change using a logistic function, or a Sigmoid [26,30]. Some reviews already cover the most commonly employed datasets and methodologies [31,32].

The most serious difficulty in autumn phenology studies using remote sensing today is validating the satellite retrievals with ground data. Both the local nature of ground observations, covering only a few selected plant species, and the shortage in ground validation data, pose a limit to validation efforts [26]. Usually, less than ten validation sites are used for regional assessments [20,24,26], or otherwise very small areas covering one or two pixels are finely covered by ground observations [7,22]. A minority of the studies does not perform any ground validation of the retrieved vegetation phenology [25], while others use dozens of ground stations to validate satellite observations across a region. For example, one study used 46 stations in China’s temperate zone during the period 1986–2005 to validate the beginning and ending of the growing season [23]. However, to the best of our knowledge, no study has used a large number of in-situ observations to validate the progression of autumn foliage across a region with regard to the different phases of foliage color progression, and not just one date that marks the end of the growing season. Therefore, the aim of this study is to monitor the progression of autumn foliage color using satellite remote sensing across a region in Southern Québec, Canada, where autumn phenological observations are frequently performed across a large number of sites, and to assess the error of the satellite retrievals against these in-situ observations.

2. Materials and Methods

2.1. Phenology Ground Observations

Ground observations used to determine foliage coloration are crucial for assessing and validating the phenology monitoring from satellite imagery. In Québec, autumn foliage color status is tracked on a weekly basis by the staff at regional tourist offices and parks in 38 locations throughout the province (Figure 1). These observations are reported in the form of an assignment to one of six classes: beginning soon (0–10% colored leaves), early (11–25% colored leaves), mid-point (26–50% colored leaves), nearly peak (51–75% colored leaves), peak (76–100% colored leaves), and past peak (end of season, most leaves have fallen to the ground). This categorical assignment is expected to overcome some of the bias that may occur as a result of different observers that estimate the leaf coloration. Weekly reports from all sites are collected by the Québec Ministère du Tourisme, and updates are available online at: https://www.quebecoriginal.com/en/discover/seasons-in-quebec#fall (accessible during the autumn season). In this study, we used available reports from five autumn seasons in 2011–2015. The reports in 2011–2013 covered 35 sites. Two more sites were covered in 2014 (Portneuf and Mont Lac-Vert), and an additional one in 2015 (Montebello).

2.2. Satellite Imagery and Pre-Processing

In this study, a time-series of satellite imagery was analyzed to obtain the temporally-normalized brownness index [21]. This index, which is described in this section, was used to model the autumn foliage coloration phase and link it with the fraction of colored leaves. The advantage of using the brownness index for this purpose is its independence of the surface background, vegetation abundance, and species composition. Therefore, it is suitable for the diverse forests in the study area, which cover over 200,000 km², ranging from boreal to deciduous and mixed forests.

Several MODIS collection 5 products [33] for two land tiles (h12v04, h13v04, each covering ~10-by-10 degrees latitude/longitude) were used in this study. The primary time series was assembled for day of year 233–305 of 2011–2015, such that it covers the entire autumn season of every year. Since MODIS has a large swath, bi-directional reflectance distribution function (BRDF) effects must be corrected to allow for temporally and spatially consistent and comparable data. Therefore, we used MCD43A4, the Nadir BRDF Adjusted Reflectance (NBAR) 16-Day L3 Global 500 m product [34]. This product is operationally produced at an interval of 8 days using cloud-free, and atmospherically-corrected surface reflectance from both Terra and Aqua satellites to calculate the surface reflectance anisotropy over a 16-day period and to retrieve a reflectance anisotropy model for each pixel. A temporal resolution of 8 days may be longer than optimal, yet this compromise in temporal accuracy is a necessary trade-off to derive vegetation indices based on NBAR. This is the reason that previous studies also used MODIS NBAR products to track autumn foliage [21,22], while studies that did not use NBAR data sometimes failed in monitoring autumn phenology [27]. Two different products were used to ensure that snow covered pixels were identified and removed from the analysis. The MCD43A2 quality product, which is complementary to the NBAR reflectance product, was used since it contains a snow albedo flag that can be used to identify snow covered pixels. A time series of the MODIS land surface temperature (LST) product (MOD11A2) that provides the average values of clear-sky LSTs during an 8-day period with a spatial resolution of 1 km [35] was used to further identify and eliminate partially snow covered pixels. Lastly, the MODIS Land Cover Type product (MCD12Q1) was used to exclude water pixels from the analysis [33].

The pre-processing chain followed the methodology of Zhang & Goldberg [8] that proved useful in tracking the progression of autumn foliage phases. In short, the NDVI was derived from NBAR bands 1 and 2 for every pixel in the time-series. Missing values were filled in by linear interpolation while ignoring trailing and leading gaps in the time series, and instead using the two nearest good quality neighbours. NDVI values in pixels that were either identified by the snow flag or corresponded to LST under 5 °C were replaced using the nearest snow free NDVI value. The NDVI value in pixels identified as water by the Land Cover Type product were set to zero. Outliers (namely unusually high NDVI values) were taken out of the analysis by a moving window median filter with a window size of 3 values that was run iteratively until the time series stopped changing.

The pre-processed MODIS NDVI values for every autumn season were fitted with a logistic equation to model the autumn foliage color progression. This approach for modeling both spring and autumn progression has been gaining popularity in recent years since it proved to be suitable for phenology monitoring in many sites around the globe [8,20,26,31,32,36,37]. This logistic function is a common type of sigmoid model that resembles an s-shaped curve in spring, and the mirror image of an s-shape in autumn. It is useful for simulating a natural process beginning with an exponential growth/decline, followed by a tapering growth/decline rate as saturation begins, indicating maturity of development (in spring) or senescence (in autumn). The equation of this model can be written as [30]:

$$\mathrm{NDVI}\left(t\right) = \frac{c}{1 + e^{\left(a + bt\right)}} + d$$

(1)

where t is time (in days), a and b are coefficients that determine the pace of greening up or senescence, c and d control the lower and upper limits of the function, such that c + d is the maximum autumn NDVI, and d is the initial background NDVI value. Parameters a and b were calculated by fitting the time series of MODIS NDVI in every pixel for every autumn season with Equation (1). Subsequently, these parameters were used to derive the temporally-normalized brownness index [21]:

$$Brownness\left(t\right) = 1 - \frac{1}{1 + e^{\left(a + bt\right)}}$$

(2)

The temporally-normalized brownness index describes the relative variation in foliage color over time in each pixel independently. Therefore, it makes the comparison of autumn phenology dynamics over large areas indifferent to ecosystem heterogeneity, surface background, plant cover type, abundance, and species composition. The relative variation in foliage color is later assigned into the same six classes as the ground observations (described in Section 2.1).

2.3. Validating the Remote Sensing Retrievals of Land Surface Phenology Using Ground Based Observations

For the purpose of comparison, both the temporally-normalized brownness time series (with an 8-day resolution), and the ground observation time series per site (with a 7-day resolution), were linearly interpolated to a daily resolution. The pixel brownness values at each observation site were extracted. The temporally-normalized brownness values were classified into the same categories used for ground observations (described in Section 2.1). The starting day of each foliage phase was determined for each site according to both satellite and ground observations in every year between 2011 and 2015. The Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) [38] over the research area was resampled to 500 m resolution using the nearest neighbour method to derive the elevation of each MODIS pixel. The linear correlation between elevation, latitude and longitude and the onset of autumn foliage coloration phases at the 38 sites was determined for both ground observations and MODIS retrievals of foliage color progression. Subsequently, the relations between elevation, latitude and longitude to MODIS foliage color progression retrievals over the entire research area were analyzed. To account for variance in autumn foliage color progression across years, the interannual phase difference in each site was determined by calculating the average and standard deviation of the start of each foliage phase for both satellite and ground observations. In addition, the Mean Absolute Error (MAE) between ground observations and MODIS retrievals was used to determine their agreement. MAEs were calculated per color progression phase for every site, as well as for all the sites cumulatively. Finally, the linear correlation between the start day of each foliage phase, as measured by satellite and ground observations, was derived.

3. Results

The temporally-normalized brownness in 2015 over the study region is visualized in Figure 2. It clearly shows that senescence begins in the northern latitudes and shifts toward the south. The ground observations of autumn foliage color progression correlate significantly with this pattern as well (

Table 1. The linear correlation between the average day of peak coloration onset and elevation, latitude and longitude. NS = not significant, (p > 0.05).

	Ground Observation		Satellite Observations
	r	p-Value	r	p-Value
Elevation	0.61	<0.001	0.35	0.029
Latitude	0.34	0.036	0.53	<0.001
Longitude	0.2	NS	0.22	NS

). This pattern was consistent during all years of the study, and is in agreement with previous studies [17,18,19,21], and with Hopkin’s Law of Bioclimatics [39]. Hopkin’s Law correlates phenological phenomena such as autumn foliage with elevation, latitude, and longitude. These variables are related to daylight hours and temperature, which are known to control leaf senescence. Table 1 shows the correlation between elevation, latitude, longitude and the average day of peak coloration onset during 2011–2015 according to both ground observations and MODIS retrievals at the pixels corresponding to the ground observation sites. While significant correlations were found between elevation and latitude to the average day of peak coloration onset, longitude was not found to be significantly correlated to it.

Moreover, the timing of the onset of each foliage coloration phase and its duration are strongly related to elevation, and latitude, but display a weak relation to longitude (Figure 3). The analysis of MODIS observations across the study area shows that high elevation pixels are associated with poor data quality, possibly because of high cloud coverage (Figure 3A,B), which results in noisy retrievals of foliage phase onset and duration. In elevations up to 500 m, as elevation increases, autumn foliage starts sooner, and the duration of each phase is shortened. In elevations above 500 m, however, the data becomes noisy. Due to this noise, the correlation between satellite data and elevation is lower than for ground data. A similar trend is shown for latitude (Figure 3C,D): at northern latitudes the onset of each autumn foliage phase is sooner, and the duration of the peak coloration phase is shortened. On the other hand, longitude does not seem to have an effect on the onset and duration of autumn foliage phases (Figure 3E,F), except for the most eastern region in the study area (the Maritimes and the Gaspé Peninsula), where autumn foliage begins later, and lasts longer compared to the rest of the study area. This is probably the result of the proximity of this region to the ocean, which moderates temperature changes, as well as differences in species composition. These findings provide additional support to Hopkin’s Law of Bioclimatics, and suggest that temperature is likely a main driver of variation in autumn foliage progression.

The average time of occurrence of foliage phases differed considerably across ground observations sites, and also between ground and satellite retrievals of these sites (

Table 2. The average start of each autumn foliage phase (day of year, rounded to the nearest integer), and the standard deviation (in brackets) at each site for 2011–2015.

Site	Ground Observations					MODIS Observations
	Early	Midpoint	Near Peak	Peak	Past Peak	Early	Midpoint	Near Peak	Peak	Past Peak
Parc national d’Aiguebelle	261 (1.53)	265 (3.3)	269 (4.16)	275 (5.31)	283 (7.04)	260 (7.3)	265 (12.3)	271 (17.4)	268 (2.89)	273 (5.66)
Témiscamingue	262 (1.15)	266 (2.28)	272 (1.41)	279 (3.36)	289 (2.65)	260 (5.46)	266 (5.46)	273 (5.41)	280 (6.06)	292 (7.22)
Témiscouata	265 (2.52)	268 (5.48)	272 (6.11)	277 (3.39)	285 (6.98)	260 (4.65)	267 (7.46)	264 (13.8)	274 (7.18)	284 (8.85)
Rivière-du-Loup/Littoral	264 (2.58)	269 (4.22)	276 (3.27)	281 (5.17)	288 (4.57)	265 (8.87)	273 (10.2)	281 (11.8)	289 (13.7)	283 (12)
Parc national du Mont-Orford	263 (3.86)	268 (5.22)	274 (5.41)	280 (5.17)	285 (4.03)	264 (6.72)	271 (6.61)	278 (8.62)	284 (11.4)	291 (8.74)
Sutton et Bromont	261 (2.94)	267 (3.36)	275 (5.36)	279 (4.83)	286 (5.48)	263 (4.72)	271 (5.4)	278 (6.35)	286 (7.26)	296 (6.18)
Mont-Mégantic	259 (1.15)	263 (3.3)	268 (5.54)	273 (5.07)	284 (6.82)	258 (14)	267 (10.9)	277 (11.2)	286 (15.5)	286 (3.21)
Bois-Francs	266 (0.707)	271 (4.18)	275 (5.29)	278 (6.06)	286 (1.83)	261 (7.27)	266 (5.93)	272 (6.73)	277 (8.93)	287 (14.8)
Nicolet/Bécancour	267 (0.837)	272 (4.6)	276 (5.39)	280 (5.68)	288 (3.83)	263 (5.03)	268 (5.24)	274 (5.74)	279 (6.38)	289 (8.64)
Baie-Saint-Paul	260 (2)	269 (4.36)	276 (3.63)	284 (3.65)	292 (3.77)	260 (5.08)	265 (4.66)	271 (4.34)	276 (4.77)	287 (5.4)
La Malbaie	260 (2)	270 (4.92)	279 (2.99)	284 (3.96)	292 (3.77)	262 (4.28)	267 (4.64)	272 (5.07)	277 (5.5)	286 (6.78)
Parc national des Grands-Jardins	260 (2)	270 (4.92)	277 (3.46)	284 (3.65)	292 (3.77)	258 (3.21)	264 (5.4)	271 (7.62)	278 (9.86)	285 (4.65)
Le littoral	260 (2.12)	264 (1.52)	271 (3.05)	275 (3.54)	284 (4.04)	263 (4.72)	271 (6.14)	278 (8.96)	285 (11.4)	291 (7.8)
Saint-Georges/Thetford Mines	260 (2.12)	266 (1.95)	271 (3.35)	276 (4.44)	283 (3.5)	262 (3.81)	266 (4.47)	270 (5.63)	275 (6.54)	282 (7.33)
Les Appalaches	260 (0)	260 (2.06)	265 (1.52)	268 (1.3)	276 (4.47)	257 (11.9)	262 (8.28)	271 (9.82)	272 (8.43)	279 (12)
Parc national de la Gaspésie	264 (3.61)	265 (5.13)	270 (4.39)	277 (5.89)	283 (6)	267 (21.1)	264 (20.4)	269 (18.7)	260 (22.9)	257 (27.7)
Parc national du Canada Forillon	264 (3.96)	270 (3.77)	274 (2.61)	278 (2.39)	285 (6.34)	260 (4.47)	265 (5.81)	270 (6.95)	276 (8.91)	286 (12.7)
Rawdon	263 (2.61)	268 (2.12)	275 (3.7)	284 (5.63)	292 (5.06)	256 (8.32)	263 (6.24)	270 (6.04)	277 (7.95)	289 (14)
Saint-Côme	261 (2.7)	266 (2.61)	272 (3.35)	280 (7.85)	290 (8.06)	262 (5.72)	267 (5.31)	271 (4.66)	277 (4.22)	285 (4.21)
Saint-Donat	260 (2.52)	264 (3.77)	270 (4.57)	274 (6.8)	288 (7.87)	260 (14.5)	262 (8.62)	269 (8.06)	276 (8.14)	288 (11)
Sud: Basses-Laurentides et Argenteuil	264 (2.97)	273 (4.98)	285 (2.79)	289 (2.99)	291 (1.41)	265 (6.04)	275 (4.55)	285 (4.09)	295 (4.32)	309 (4.58)
Centre: Cœur des Laurentides	261 (0.577)	264 (2.24)	269 (4.15)	273 (4.98)	280 (4.83)	262 (5.34)	266 (6.3)	270 (6.83)	275 (7.98)	282 (10.9)
Nord: Hautes-Laurentides	260 (1.41)	265 (2.07)	271 (0.837)	278 (4.39)	285 (5.64)	261 (7.38)	267 (7.23)	273 (7.12)	279 (7.75)	290 (8.38)
Centre de la nature	266 (4.53)	276 (6.43)	283 (3.95)	287 (4.11)	295 (4.86)	263 (6.72)	273 (5.22)	283 (4.44)	293 (4.39)	308 (4)
Parc national du Canada de la Mauricie	261 (2.52)	269 (2.38)	273 (1.95)	278 (2.19)	293 (1.71)	261 (7.76)	266 (6.84)	270 (5.32)	275 (4.82)	284 (4.15)
Trois-Rivières	264 (3.06)	270 (2.38)	273 (2.77)	279 (3.77)	293 (1.63)	264 (7.73)	273 (5.93)	282 (4.82)	292 (4.77)	305 (7.23)
Mont Saint-Hilaire	262 (2.22)	272 (4.51)	283 (3.21)	287 (1.5)	295 (1.5)	267 (7.5)	275 (5.93)	283 (4.66)	292 (3.58)	308 (3.36)
Rougemont	260 (1.82)	269 (3.83)	279 (5.1)	286 (5.19)	296 (1.91)	267 (6.69)	276 (5.59)	284 (4.51)	293 (3.58)	308 (3.44)
Mont Rigaud	264 (4.55)	272 (4.44)	281 (3.97)	287 (3)	294 (1.63)	274 (6.44)	281 (4.97)	287 (4.1)	294 (3.42)	306 (3.56)
Mont Sainte-Marie	264 (3.06)	271 (3.81)	276 (3.9)	281 (3.61)	288 (3.5)	263 (3.63)	269 (4.18)	276 (5.39)	283 (6.02)	295 (8.41)
Parc de la Gatineau	266 (3.77)	271 (3.77)	276 (3.27)	281 (3.27)	288 (3.5)	262 (7.66)	270 (5.32)	277 (3.83)	284 (3.29)	298 (6.53)
Lac Beauport/Stoneham	263 (0)	263 (3.63)	270 (2.28)	274 (3.05)	282 (4.93)	260 (3.81)	268 (6.07)	277 (8.53)	286 (11.3)	295 (8.04)
Mont Sainte-Anne	263 (0)	263 (3.63)	270 (2.35)	273 (3)	282 (4.93)	261 (5.18)	265 (6.02)	271 (6.78)	276 (7.4)	285 (9.82)
Parc national du Saguenay (secteur Rivière-Éternité)	262 (1.15)	265 (1.82)	270 (2.92)	274 (3.08)	283 (6.85)	263 (8.14)	267 (7.11)	271 (6.73)	276 (6.31)	280 (6.23)
Parc national des Monts-Valin	260 (2.45)	263 (2.99)	267 (3.3)	271 (3.96)	279 (3.94)	273 (16.3)	273 (20.1)	275 (25)	292 (9.29)	302 (4.36)
Portneuf	264 (4.24)	269 (7.07)	273 (7.07)	277 (6.36)	276 (0)	258 (9.28)	264 (8.43)	271 (7.98)	277 (8)	289 (8.83)
Mont Lac-Vert	262 (1.41)	266 (2.12)	275 (9.19)	284 (2.12)	289 (0)	262 (5.06)	262 (12.5)	273 (7.89)	276 (12.4)	283 (20.5)
Montebello	264 (0)	276 (0)	283 (0)	286 (0)	n/a	266 (5.59)	273 (4.83)	281 (5.18)	288 (5.15)	301 (5.5)

, Figure 4). The interannual phase difference at each site during 2011–2015 (as represented by the standard deviation of the starting day) is much larger for MODIS retrievals compared to ground observations. In the early, midpoint, near peak, peak and past peak coloration phases the standard deviations range between 0–5, 0–7, 0–9, 0–8, and 0–8 days, respectively, for ground observations, and 3–21, 4–20, 4–25, 3–23, 3–28 days, respectively, for MODIS retrievals. The interannual standard deviation of occurrence of the different foliage color phases as observed on the ground and reported by Zhang & Goldberg (2011) is 1.7–3.6 for Harvard Forest in 2001–2004, and 3–6.6 for Hubbard Brook Experimental Forest in 2001–2003, which is within the range found for the sites in Québec. Zhang & Goldberg (2011) also found that the interannual variation in the timing of peak and near-peak colorations is less than 10 days in more than 90% of North-Eastern America for 2001–2004. However, our results for 2011–2015 show greater interannual variation: only 10% of the research area shows an interannual variation of less than 10 days in the timing of each foliage stage, and the interannual variation in the timing of peak and near-peak colorations is up to 39 days in more than 90% of the research area.

In most sites, a significant correlation was found between ground observations and MODIS retrievals (

Table 3. Mean absolute error (MAE) and the linear correlation between ground observations and MODIS retrievals of the start day of each autumn foliage phase. Significant R² (p ≤ 0.01) are marked in bold.

Site	Early	Midpoint	Near Peak	Peak	Past Peak	R2
Parc national d’Aiguebelle	5.33	9.50	9.40	6.75	12.25	0.11
Témiscamingue	5.00	4.60	3.80	4.80	7.00	0.74
Témiscouata	4.67	3.50	8.00	3.75	4.33	0.78
Rivière-du-Loup/Littoral	7.75	7.00	8.60	11.80	4.50	0.50
Parc national du Mont-Orford	6.25	4.00	4.40	6.80	7.75	0.74
Sutton et Bromont	5.00	5.00	3.20	6.80	10.25	0.84
Mont-Mégantic	14.00	7.75	9.00	13.00	7.33	0.59
Bois-Francs	7.80	4.80	5.40	7.00	14.50	0.54
Nicolet/Bécancour	5.80	5.00	3.60	2.80	2.50	0.83
Baie-Saint-Paul	5.33	7.00	7.20	8.40	7.00	0.57
La Malbaie	6.00	4.00	7.75	8.20	8.25	0.59
Parc national des Grands-Jardins	2.33	6.25	8.40	11.00	8.25	0.57
Le littoral	5.00	6.40	7.20	9.80	7.25	0.78
Saint-Georges/Thetford Mines	4.40	2.00	2.20	4.00	6.75	0.68
Les Appalaches	14.00	5.50	7.60	6.50	7.75	0.35
Parc national de la Gaspésie	12.00	20.40	18.80	19.75	27.67	0.01
Parc national du Canada Forillon	4.60	6.60	6.40	7.00	9.50	0.55
Rawdon	8.40	6.60	7.20	11.00	13.00	0.54
Saint-Côme	4.80	3.20	3.80	8.00	7.00	0.66
Saint-Donat	12.67	4.67	5.33	8.50	12.33	0.37
Sud: Basses-Laurentides et Argenteuil	4.60	2.20	3.00	4.75	12.00	0.84
Centre: Cœur des Laurentides	4.00	3.00	2.60	3.80	5.00	0.68
Nord: Hautes-Laurentides	5.25	4.40	4.40	6.20	6.80	0.62
Centre de la nature	5.60	4.60	2.75	4.75	12.33	0.81
Parc national du Canada de la Mauricie	7.67	6.25	3.80	3.20	10.50	0.62
Trois-Rivières	4.00	5.00	8.60	12.80	11.00	0.88
Mont Saint-Hilaire	5.75	6.80	2.20	3.50	11.25	0.79
Rougemont	6.80	7.20	5.40	5.50	12.25	0.82
Mont Rigaud	9.50	10.60	6.80	6.75	10.75	0.73
Mont Sainte-Marie	2.33	2.00	1.40	3.60	3.75	0.92
Parc de la Gatineau	7.25	3.60	2.60	4.00	9.75	0.81
Lac Beauport/Stoneham	3.50	5.60	7.20	12.20	13.50	0.76
Mont Sainte-Anne	5.00	4.60	3.80	4.60	5.50	0.61
Parc national du Saguenay (secteur Rivière-Éternité)	7.25	4.80	4.20	4.00	5.00	0.50
Parc national des Monts-Valin	15.00	16.67	18.00	18.67	20.67	0.24
Portneuf	8.50	5.00	4.00	6.50	13.00	0.92
Mont Lac-Vert	5.00	5.00	6.00	12.50	24.00	0.26
Montebello	10.00	3.00	3.00	7.00	n/a	0.89
All sites	6.52	5.87	5.95	7.48	9.50	0.48

). Figure 5 shows an example of a relatively good correlation between the two types of phenological monitoring at the Mont Sainte-Marie site. However, the temporally-normalized brownness index mismatches the ground observations of foliage progression in four locations (Gaspésie National Park, Monts-Valin National Park, Mont Lac-Vert, d’Aiguebelle National Park). The linear correlation between ground observations and MODIS retrievals of the start day of each autumn foliage phase is not significant in these four sites. The MAE of each autumn foliage phase in these sites is also relatively high (Table 3), and therefore the MODIS retrievals at these sites are considered as outliers. This mismatch might be the result of a highly clouded environment over those sites that degraded the quality of the MODIS time series. Additionally, both Gaspésie National Park and Monts-Valin National Park are high elevation sites which, as shown in Figure 3A,B, are prone to data gaps due to cloudiness and poor quality imagery. Figure 6 shows the scatterplots summarizing the onset day of autumn foliage coloration phases between in-situ measurements and MODIS retrievals for 2011–2015 at all sites. There is a highly significant positive linear correlation between ground observations and MODIS retrievals (R² = 0.48, p < 0.001), and this correlation is further improved when excluding the four above mentioned outliers from the analysis (R² = 0.58, p < 0.001). The MAE across all sites and all coloration phases is 7 days, and improves to 6.32 days by removing the outliers (6.23, 5.24, 5.21, 6.9, 8.55 days for the early, midpoint, near peak, peak and past peak coloration phases). Accordingly, the MAE for the early, midpoint, near peak, and peak coloration phase are all lower than the temporal resolution of the ground observations (7 days) and therefore acceptable, while the MAE at the past peak coloration phase is slightly higher. Figure 7 shows the MAE per site at each coloration phase as a function of the site elevation, and demonstrates that high elevation sites have the highest retrieval errors.

4. Discussion

The results presented here are the first attempt to validate the MODIS based temporally-normalized brownness index for autumn foliage progression monitoring over a period of five years using in-situ observations at a large number of sites across Quebec, Canada. This validation is necessary in order to reliably track this phenomenon in a spatially continuous manner across the region, thus extending monitoring currently performed in discrete locations. While previous studies demonstrated a very good correspondence between a few ground validation sites and MODIS retrievals of autumn foliage (

Table 4. Comparison between different studies that used ground observations and MODIS retrievals of autumn foliage progression.

Study	Number of Sites	Number of Years	MAE
Hmimina et al. (2013)	7	9	4–15
Zhang & Goldberg (2011)	2	8	3.9–8.3
Liu et al. (2015)	2	2	5.25
Che et al. (2014)	10	10	5–16
Hufkens et al. (2012)	27	5-10	14–36
Kobayashi et al. (2016)	17	1-3	16–26
This study	38	5	5.24–8.55

), by validating the MODIS retrievals at a large number of sites, we show that this is not always the case. Most sites did show a good agreement between ground-based foliage progression reports and MODIS retrievals. However, a few sites where the NDVI time series was degraded exhibit errors in foliage progression estimations during some or all years of the study. Therefore, it is still possible to improve our ability to automatically track autumn foliage changes, and further research should be performed to improve the detection accuracy of autumn foliage phases. Overcoming this challenge would allow for near-real-time tracking of autumn foliage progression, and the exclusion of erroneous predictions.

The considerable improvement in the linear correlation between ground observations and MODIS retrievals (from R² = 0.48 to R² = 0.58), as well as the improvement in MAE (from 7 to 6.32), that were achieved by removing outliers from the analysis, show that there is potential for improvement in the MODIS retrieval of autumn foliage status. Our reliance on 16-days MODIS products to derive the NDVI time series could cause some uncertainties in the classification of autumn foliage phases. Moreover, cloud contamination influences the representativeness of the NBAR reflectance used to derive the NDVI time series. Accordingly, in the future, protocols for assessing the classification confidence in each pixel should be developed. Such protocols could aid in the elimination of pixels with poor data quality, and thus exclude some noise from the analysis.

Another unavoidable caveat of this study is the use of the ground observations as ground truth. Unlike previous studies that used few sites with strong control over field data collection protocols [9,20,21,22], our study used data collected in multiple sites with little control over observation protocols and accuracy. The classification of foliage progression into six coloration phases (beginning soon, early, midpoint, near peak, peak and past peak) may reduce observation errors of the fraction of colored foliage, however, because the observations are performed by different personnel, at different locations, and where the vegetation composition is different, this can result in some bias in those ground estimations. As a result, our ability to confidently extrapolate these subjective observations over large areas is limited and varies between sites. On the other hand, it is quite remarkable that a significant correlation exists between MODIS pixels that are 25 ha in area, and in-situ observations that probably cover only a fraction of that area. Therefore, bias in ground estimations of autumn foliage progression might affect the apparent error in satellite based estimations of foliage progression. Other studies that used ground red-green-blue (RGB) cameras [14,15,20] to validate satellite phenology retrievals, encountered similar issues with scale and representativeness of images acquired by ground cameras. To a certain extent, this inhibited the correlation between ground and satellite remote sensing of forest phenology. Therefore, it seems that some degree of inconsistency between ground and spaceborne observations of autumn foliage progression is unavoidable, but should be minimized as much as possible.

On a regional scale, the timing of autumn foliage onset was found to be related to elevation and latitude. This is in agreement with previous studies that noted latitudinal gradients of the autumn onset date [17,18,19]. Additional factors that are likely associated with the autumn foliage progression pattern, such as microclimate and plant species composition, were not studied here. However, it has been established that the timing of change in leaf color varies considerably between different species [21]. Spaceborne observations are often characterized by a tradeoff between spatial coverage, spectral resolution, and temporal resolution. In this study, the use of MODIS NBAR products allowed for the coverage of a large area with a relatively high temporal resolution at the expense of spatial resolution, and thus the analysis of individual trees was not possible. Should a smaller area be the focus of investigation, a time series of high (i.e., Sentinel-2, SPOT, RapidEye), and very-high (i.e., IKONOS, WorldView satellites, airborne sensors), spatial resolution images could be used to classify tree species [40] and thus track autumn foliage progression for individual tree species.

Since the interannual variation found by Zhang & Goldberg [8] for 2001–2004 is considerably smaller than what we have found for 2011–2015, an analysis of a longer time series would be useful to explore. MODIS records extend from 2000 until the present, and therefore in a future study, this data could be used to examine trends in autumn foliage progression over 17 years. It would be particularly interesting to examine the autumn foliage onset and progression with respect to climate variation since it is known to affect green leaf phenology [20]. Phenological changes are also related to climate change since changes in the length of the growing season can contribute to changes to climate drivers such as albedo, and also atmospheric CO₂ concentrations due to carbon sequestration in terrestrial plants [41]. In addition, significant changes in the peak of autumn foliage can also influence the tourism industry with regards to the best timing for foliage trips and their promotion. In China, climate changes resulted in an average delay of the autumn foliage vacation season by 4–5 days per decade over a span of 30 years [12]. Climate change effects on autumn foliage could impose serious challenges to the tourism industry in North America as well. For example, changes in species composition from maple-beech-birch forests that are responsible for some of the spectacular autumn colored landscapes, into oak-hickory as a result of warming, is expected to reduce the bright red coloration that is the attraction for many visitors [11]. In order to study these effects in Québec, a longer time-series is needed, but assuming that Québec is not immune to regional and global warming trends, the regional timing of autumn foliage is likely to be affected.’

5. Conclusions

This paper presents the first regional scale validation of the MODIS temporally normalized brownness index for tracking autumn foliage phases using 38 validation sites. We showed that spaceborne data significantly correlates with ground observations in most validation sites (R² = 0.58), and the MAE is smaller than the frequency of these observations (MAE = 6.32). However, this was not the case in four of the validation sites, where the MODIS time series was degraded by clouds, and as a result the foliage coloration status was misestimated. Therefore, future research emphasis should be placed on estimating the MODIS time-series quality or the foliage status classification confidence level, in order to exclude noisy pixels from the analysis and reliably track the foliage color progression through autumn. In spite of the uncertainties in regional tracking of this phenomena from space, this type of study is important for the local tourism economy, as well as for stimulating further examination of phenology shifts in climate models and to improve their forecast abilities by introducing local scale findings.