Articles | Volume 18, issue 11
Research article
07 Jun 2021
Research article |  | 07 Jun 2021

A survey of proximal methods for monitoring leaf phenology in temperate deciduous forests

Kamel Soudani, Nicolas Delpierre, Daniel Berveiller, Gabriel Hmimina, Jean-Yves Pontailler, Lou Seureau, Gaëlle Vincent, and Éric Dufrêne

Tree phenology is a major driver of forest–atmosphere mass and energy exchanges. Yet, tree phenology has rarely been monitored in a consistent way throughout the life of a flux-tower site. Here, we used seasonal time series of ground-based NDVI (Normalized Difference Vegetation Index), RGB camera GCC (greenness chromatic coordinate), broadband NDVI, LAI (leaf area index), fAPAR (fraction of absorbed photosynthetic active radiation), CC (canopy closure), fRvis (fraction of reflected radiation) and GPP (gross primary productivity) to predict six phenological markers detecting the start, middle and end of budburst and of leaf senescence in a temperate deciduous forest using an asymmetric double sigmoid function (ADS) fitted to the time series. We compared them to observations of budburst and leaf senescence achieved by field phenologists over a 13-year period. GCC, NDVI and CC captured the interannual variability of spring phenology very well (R2>0.80) and provided the best estimates of the observed budburst dates, with a mean absolute deviation (MAD) of less than 4 d. For the CC and GCC methods, mid-amplitude (50 %) threshold dates during spring phenological transition agreed well with the observed phenological dates. For the NDVI-based method, on average, the mean observed date coincides with the date when NDVI reaches 25 % of its amplitude of annual variation. For the other methods, MAD ranges from 6 to 17 d. The ADS method used to derive the phenological markers provides the most biased estimates for the GPP and GCC. During the leaf senescence stage, NDVI- and CC-derived dates correlated significantly with observed dates (R2=0.63 and 0.80 for NDVI and CC, respectively), with an MAD of less than 7 d. Our results show that proximal-sensing methods can be used to derive robust phenological metrics. They can be used to retrieve long-term phenological series at eddy covariance (EC) flux measurement sites and help interpret the interannual variability and trends of mass and energy exchanges.

1 Introduction

In the temperate and boreal climate zone, the timing of phenological events is strongly controlled by temperature and is thus responsive to the ongoing climate change (Menzel et al., 2006; Badeck et al., 2004; Piao et al., 2019). The opening of buds (“budburst”) in spring and the coloration and fall of leaves (“leaf senescence”) in autumn are the key steps in the phenological cycle of forest trees. These stages mark the start and end of the photosynthetically active period and as such strongly influence the carbon and water exchanges between the ecosystem and the atmosphere (Goulden et al., 1996; Delpierre et al., 2009a; Richardson et al., 2010; Dragoni et al., 2011). Historically, the timing of these events has been monitored through direct and periodic human-eye observations of the state of buds and leaves in the field (Sparks and Carey, 1995). However, this method is time-consuming, laborious and subject to an observer effect (Roetzer et al., 2000; Schaber and Badeck, 2002; Klosterman et al., 2014). Alternative, ground-based indirect methods have been tested for monitoring the phenology of different ecosystems. Proximal-sensing methods based on measuring the radiation reflected, transmitted or absorbed by the canopy (henceforth “radiation-based methods”) are increasingly being used. The broadband NDVI (Normalized Difference Vegetation Index) calculated from measurements of the fraction of reflected radiation in the photosynthetically active radiation (PAR) spectral domain and shortwave bands, proposed by Huemmrich et al. (1999), has been successfully used in order to monitor vegetation phenology in many studies (Huemmrich et al., 1999; Richardson et al., 2007; Liu et al., 2019). Forest phenology was also described from measurements of the fraction of transmitted PAR through the canopy (Toda and Richardson, 2018; Perot et al., 2019) and leaf area index (LAI) (Keenan et al., 2014). Spectral vegetation indices derived from tower-mounted hyperspectral spectroradiometers (Kobayashi et al., 2018; Lu et al., 2018), RGB/IR (red, green, blue and infrared bands) cameras (Richardson et al., 2007; Klosterman et al., 2014; Richardson et al., 2018a; Richardson, 2019; Milliman et al., 2019), or from two-band red and near-infrared proximal sensors (Ryu et al., 2010; Eklundh et al., 2011; Soudani et al., 2012; Hmimina et al., 2013) have also been assessed. More recently, passive sun-induced fluorescence has been used (Lu et al., 2018). At vegetation sites where continuous measurements of carbon flux are available, phenology has also been estimated from the dynamics of GPP (gross primary productivity) and net ecosystem exchange (NEE) (Gonsamo et al., 2013; Wu et al., 2017; Garrity et al., 2011).

Over the past 2 decades, hundreds of experimental sites measuring CO2, water and energy exchanges between ecosystems and the atmosphere have been set up worldwide. These sites are organized in networks (FLUXNET, ICOS (Integrated Carbon Observation System), etc.) and aim to record long-term data according to standardized protocols (Baldocchi et al., 2001; Franz et al., 2018). These sites acquire high temporal-resolution time series combining both mass (CO2 and water) flux data with ancillary data, which include incident, reflected and transmitted radiation measurements in different spectral ranges, and also LAI, NDVI and RGB images of the canopy. Yet, the phenology of the vegetation cover is not routinely monitored over all sites, precluding the assessment of its influence on carbon and water exchanges. These sites provide data which allow the comparison of various radiation-based methods for monitoring forest phenology. However, the comparative studies cited above and those carried out at some of the carbon flux measurement sites did not cover all the methods at the same site and were also limited to a few and short periods of time. Also, most of these studies suffered from a lack of direct and independent phenological observations. As underlined in Klosterman et al. (2014), this is a key challenge in interpreting estimates from the various approaches. Indeed, most of the radiation-based methods use optical signals at different wavelengths and at different spectral resolutions. Depending on species and sensor specifications (spectral, radiometric and geometric responses), this could lead to possible mismatches between observed and estimated phenology due to the well-known selective absorption properties of plant components (Sims and Gamon, 2002). The measurement conditions (sun-view geometry, field of view) may also differ (Sonnentag et al., 2012). Also, some mainly observe the top of the canopy (down-looking sensors mounted above the canopy), while others are more integrative of the whole canopy (indirect methods that use transmitted or absorbed radiation). Therefore, there is a need to conduct comparative studies to establish rigorously the correspondence between phenological dates recorded by field phenologists and phenological metrics predicted by indirect proximal methods.

In this study, we present an exhaustive comparative survey of various proximal methods to estimate both spring and autumn phenology in a mature deciduous forest ecosystem surrounding the Fontainebleau–Barbeau carbon flux tower. The main objective is to evaluate the performance of each of the methods in reproducing interannual variation in spring and autumn phenology directly observed by field phenologists over a 13-year period.

2 Materials and methods

2.1 Site description

Data were mainly acquired at the eddy covariance (EC) flux measurement site at the Fontainebleau–Barbeau forest (482826′′ N, 24657′′ E), 53 km southeast of Paris, France. Fontainebleau–Barbeau is a deciduous forest mainly composed of mature sessile oak (Quercus petraea (Matt.) Liebl) and an understory of hornbeam (Carpinus betulus L.). The average stand LAI, based on measurements using the litter collection method over the 2012–2018 period, is 5.8 m2 m−2, ranging from 4.6 to 6.8 m2 m−2 (unpublished data). Hornbeam contribution to stand LAI accounts for 30 %, ranging from 24 % to 39 % from year to year.

At the Fontainebleau–Barbeau EC flux measurement site, which belongs to the European ICOS-RI Ecosystem network (Integrated Carbon Observation System – Research Infrastructure, FR-Fon code), a 35 m high tower was installed in 2005 in order to measure energy and CO2 exchanges between the forest and the atmosphere with the eddy covariance technique. More details about the study site and flux calculation are given in Delpierre et al. (2016). The tower has been equipped with various proximal sensors that we used here to estimate the timings of phenological events (Table 1). More details about the instrumentation and measurements achieved at this site are available at, last access: 27 May 2021).

Table 1Methods and variables used in the calculation of phenology metrics in the Fontainebleau–Barbeau forest. NDVI: narrow-band Normalized Difference Vegetation Index; NDVIbr: broadband NDVI; fRvis: fraction of reflected radiation by the canopy in the PAR spectral domain; GCC: greenness chromatic coordinate from RGB camera images; fAPAR: fraction of absorbed radiation in the PAR spectral domain; CC: canopy closure; LAI: leaf area index; GPP: gross primary productivity. These vegetation variables are named Vv hereafter.

* See text for details.

Download Print Version | Download XLSX

2.2 Extraction of phenological markers

The data and methods used in the calculation of phenology metrics are summarized in Table 1. The general principle of the phenological-metrics extraction method consists in building time series at daily resolution that describe the canopy foliage dynamics during the whole seasonal cycle of vegetation (Fig. 1). This method applies to all the variables (“vegetation variable”, Vv) listed in Table 1. Then, to compare the different vegetation proxies without possible methodological biases, we opted for the same method using an asymmetric double sigmoid (ADS) similar to Zhang et al. (2003), Soudani et al. (2008) and Klosterman et al. (2014).

Figure 1Illustration of phenological markers extracted from ADS (asymmetric double sigmoid) functions fitted to NDVI data acquired in 2015 (empty circles and red curve). Vertical lines in blue: SOS, MOS and EOS are dates of start, middle and end of leaf onset in spring. SOF, MOF and EOF are dates of start, middle and end of leaf senescence (colored and fallen leaves) in autumn. The third derivative of the ADS function shows peaks and holes corresponding to the six phenological dates (black dotted line).


Briefly, an asymmetric double sigmoidal function was fitted on Vv time series according to the following equation:

(1) V v t = w 1 + w 2 + 1 2 w 1 - w 2 [ tanh w 3 t - u - tanh w 4 t - v ] .

Vv (t) is the considered vegetation variable (percentage of open buds and percentage of non-senescent leaves, NDVI, NDVIbr, fRvis, fAPAR, CC, LAI, GCC (greenness chromatic coordinate) or GPP). t is the time (day of year). tanh is the hyperbolic tangent and w1, w2, w3, w4, u and v are the fitting parameters. (w1+w2) is the Vv minimum in the unleafy season. (w1w2) is the total amplitude of variation in Vv over the year. The two phenological markers u and v are the dates of the two inflection points when Vv increases during the spring (u) and decreases during the autumn (v). For these two dates u and v, Vv (t) is very close to 50 % of its total amplitude of variation in spring and autumn, respectively. Four other phenological markers are determined numerically from the extrema of the third derivative of the ADS function according to Zhang et al. (2003) (Fig. 1). The six phenological markers are named as follows according to Klosterman et al. (2014): SOS, MOS and EOS for the start, middle and end of leaf onset (budburst) in spring and SOF, MOF and EOF for the start, middle and end of leaf senescence in autumn, corresponding approximately to 10 %, 50 % and 90 % of total amplitude during the increase and the decline in canopy greenness in spring and autumn, respectively.

Fitting was done by minimizing the sum of squares of differences between fitted (Eq. 1) and measured Vv. In order to better constrain the fitting at the end of the leafy season, each year of data was extended to the end of January of the following year. Thus, potentially, each time series is composed of 396 d instead of 365 d.

2.3 Data

2.3.1 Field phenological observations (OBS)

We collected spring and autumn phenological field observations at the Fontainebleau–Barbeau forest over 13 years (2006–2018; see Delpierre et al., 2020, Denéchère et al., 2019) through complementary sampling schemes. Over the 2006–2018 period, we implemented an “extensive” survey in which, biweekly over March–April, we monitored the bud development of >100 randomly chosen dominant sessile oak trees and recorded the date at which 50 % of the individuals displayed at least 50 % open buds in their crowns (corresponding to stage 7 of the BBCH (Biologische Bundesanstalt für Land und Forstwirtschaft) scale). Observations were done with binoculars by three intercalibrated observers. This date is referred to as BB-OBS (BB for budburst) in the following. In years 2015–2017 we complemented this protocol with an “intensive” survey. In the plot, 27 to 66 individual trees (depending on years) were tagged and monitored for budburst from 0 % budburst to 100 % budburst. This survey yielded the progress of budburst for each tree crown, which we averaged to get the progress of budburst for the tree population (Fig. 2a). We further monitored the progress of leaf senescence (percentage of colored or fallen leaves) in each individual tree crown weekly in autumn and averaged the individual values to get the progress of leaf senescence at the tree population scale (Fig. 2a). We fitted the ADS (Eq. 1) function to these continuous data and retrieved the MOS-OBS (in spring) and MOF-OBS (in autumn) metrics. The MOS-OBS (obtained from the intensive survey) and BB-OBS (obtained from the extensive survey) dates compare very well, their maximum absolute difference being 1 d (Delpierre et al., 2020). Hence in the following we will use the BB-OBS as the observed date of budburst over the whole (2006–2018) study period. All spring phenological observations were conducted on a biweekly basis. Hence the uncertainty of BB-OBS is 3.5 d.

We completed the MOF-OBS (autumn) metrics obtained at Fontainebleau–Barbeau through the intensive survey over 2015–2017 with leaf senescence data obtained over 2011–2014 from a phenological survey site 50 km away from Fontainebleau–Barbeau (Orsay site). At this site, we deployed an intensive-monitoring protocol of leaf senescence (30 to 60 tagged sessile oaks monitored weekly for the percentage of colored or fallen leaves during autumn) from which we obtained the LS-OBS metrics, that is the date at which 50 % trees had 50 % leaves colored or fallen. In 2015, autumn phenological observations were conducted simultaneously in Fontainebleau–Barbeau and Orsay: the MOF-OBS (Fontainebleau–Barbeau, DOY 300) and LS-OBS (Orsay, DOY 295) dates compared well. Considering that leaf senescence dates are comparable at a scale of tens of kilometers (Delpierre et al., 2009b), we used the 2011–2014 Orsay LS-OBS data to complement the 2015–2017 Fontainebleau–Barbeau MOF-OBS data. All spring phenological observations were conducted on a weekly basis. Hence the uncertainty of MOF-OBS and LS-OBS is 7 d.

2.3.2 Narrow-band NDVI

The NDVI is calculated as follows:

(2) NDVI = ( NIR - R ) / ( NIR + R ) .

R and NIR are radiances in the red (640–660 nm) and the near-infrared (780–920 nm) bands, respectively. Radiances are measured using a laboratory-made NDVI sensor (Pontailler et al., 2003). A description of this sensor and its use for estimating phenological metrics in various biomes is given in Soudani et al. (2012) and Hmimina et al. (2013). Briefly, the sensor is positioned at the top of the EC flux tower in the Fontainebleau–Barbeau forest, about 7 m above the canopy, directed downwards and inclined about 20–30 to the vertical and facing south to avoid the hot-spot effects in canopy reflectance when the viewing direction is collinear with the solar direction. The field of view of the sensor was 100 and the area observed is a few tens of square meters. Measurements are acquired continuously every half hour. Noisy data, mainly due to rainfall and very low radiation conditions, were removed according the procedure described in Soudani et al. (2012). This procedure consists in keeping only NDVI measurements recorded when the ratio between global radiation (RGin) measured above the canopy and the exo-atmospheric radiation (Rex) at the top of atmosphere exceeds the threshold of 0.6, considered to be the threshold for distinguishing between clear and overcast sky conditions (Soudani et al., 2012). Then, the daily average of filtered NDVI data acquired between 10:00 and 14:00 UT is considered to minimize the effects of daily variations in solar angle. Finally, filtered and daily averaged NDVI data were used in Eq. (1).

2.3.3 RGB camera

Digital pictures (resolution of 2590×1920 pixels) of the forest canopy are acquired continuously every hour between 08:00 to 17:00 UT with an AXIS P1347 camera installed next to and according to the same geometric configuration of the NDVI sensor. In order to minimize effects of changing illumination conditions, a white PVC panel is installed in the camera field of view (FOV) and used as a reference. Pictures (10 d−1) were processed automatically under MATLAB. At first, three regions of interest (ROIs) were delineated on a spring picture. Two ROIs, having an area of 3000 pixels and 1140 pixels, respectively, are located on the reference panel. The third ROI is located over the vegetation area that covers the central region of the picture (2 MP pixels). To convert RGB data measured by the camera to pseudo-reflectance (ρR, ρG, ρB), digital counts in red, green and blue bands of the vegetation ROI were averaged and divided by the averages of R, G and B measured on the two white ROIs on the reference PVC panel. These pseudo-reflectances were averaged on a daily basis (10 values per day, corresponding to the hourly sampling) and used to determine the daily GCC as follows:

(3) GCC = ρ G / ρ R + ρ G + ρ B .

Phenological markers are then extracted from GCC time series according Eq. (1).

2.3.4 Broadband NDVIbr and fraction of reflected radiation fRvis

Broadband NDVI (NDVIbr), named according to Huemmrich et al. (1999), was calculated from incoming and reflected radiation in the visible spectral region (400–700 nm) corresponding to the spectral range of PAR measured using PAR sensors (PQS1, Kipp and Zonen, Finland) and in the shortwave spectral regions (200 to 3600 nm) using a CMP22 pyranometer (Kipp and Zonen, Finland). A conversion factor of 4.57 µmol J−1 (McCree, 1972, in Wang et al., 2006) was used to convert PAR units (µmol m−2 s−1) to energy units (J m−2 s−1). As in Wohlfahrt et al. (2010), NDVIbr is calculated as below:

(4) NDVI b r = NIR out NIR in - PAR out PAR in NIR out NIR in + PAR out PAR in .

NIRin= RGin PARin, NIRout= RGout PARout, RGin, RGout, PARin and PARout are incoming and outgoing reflected radiation in shortwave and PAR spectral regions.

The fraction of reflected radiation fRvis was calculated as

(5) f R vis = PAR out PAR in .

NDVIbr and fRvis were filtered by applying the same ratio of 0.6 between RGin and Rex and limiting the period of acquisition to between 10:00 to 14:00 h UT. Finally, filtered and daily averaged fRvis and NDVIbr data were used to in Eq. (1) to extract the six phenological markers. Because fRvis was lower during the leafy season than in winter (unleafy season), Eq. (1) was applied to (1−fRvis) allowing it to have the same temporal pattern as the other variables. For simplicity, the fRvis term will be used hereafter when referring to the method.

2.3.5 Fraction of absorbed PAR (fAPAR), canopy closure (CC) and leaf area index (LAI)

Fifteen quantum PAR sensors (PQS1, Kipp and Zonen, Finland), directed towards the sky, are installed below the canopy on the ground area surrounding the EC flux tower to ensure a robust spatial sampling of the radiation transmitted through the canopy. Measurements are achieved at a half-hour time step, simultaneously with measurements of incoming and reflected PAR radiation above the canopy. The filtering of transmitted, reflected and incoming radiation measurements is carried out according to the same procedure used for NDVI, NDVIbr and fRvis. Consequently, only measurements taken between 10:00 and 14:00 h UT after filtering are used in the calculation of fAPAR, CC and LAI.

fAPAR is calculated according the following expression:

(6) f APAR = PAR in - PAR out - PAR t PAR in .

CC is calculated using a new formulation as follows:

(7) CC = 1 - PAR t cos θ / PAR in ,

where PARin and PARout are defined above in Eq. (3). PARt is the averaged over 15 sensors of transmitted radiation measured beneath the canopy. θ is the sun zenith angle calculated using the standard astronomical formula. Unlike Eq. (6) and the previous studies (Richardson et al., 2007; Garrity et al., 2011; Toda and Richardson, 2018), the division of PARt by the cosine of the sun zenith angle (Eq. 7) allows us to consider variation in PARt due solely to the variation in the path length of incident radiation passing through the forest canopy before reaching the ground according to the seasonal variation in the solar angle. In order to assess the performance of this new formulation proposed in this study, we also calculated CC without cosine correction.

Another possible alternative to this correction/normalization in order to take into account sun angle effects on transmitted PAR (Eq. 7) is to estimate leaf area index from the canopy gap fractions since the estimation of LAI using the Beer–Lambert law corrects for the effects of solar angle and considers leaf angle distribution through the extinction coefficient K. The LAI was calculated as follows:

(8) LAI = - log PAR t / PAR in / K .

Log is the natural logarithm. K is the coefficient of extinction, calculated following the expression given in Campbell and Norman (1998):

(9) K θ = x 2 + tan θ 2 x + 1.774 ( x + 1.182 ) - 0.733 .

The parameter x describes an ellipsoidal leaf angle distribution function (x=1 for spherical distribution, x>1 for planophile and x<1 for erectophile leaves). In this study and in order to let K vary according to the seasonal variations in the solar angle, we only fixed the parameter x in Eq. (9). In order to estimate an average value of the x parameter in the Fontainebleau–Barbeau forest, Eq. (8) was inverted, based on direct LAI measurements around the EC flux tower using a litter collection technique according to the ICOS protocol (Gielen et al., 2018) and the radiation measurements over the 2012–2018 period. The parameter x was about 1.4 which corresponds to an average value of K of about 0.67 during the leafy season (DOY 150–240). This value agrees with previous studies (Baldocchi et al., 1984; Holst et al., 2004). Thus, we note that K is calibrated from the “true” average green LAI measured by the litter collection method, and thus it corrects for clumping effects and woody components. The term LAI is used in the present study instead of the term PAI (plant area index, including leaf and woody components) usually used when it is estimated from canopy transmittance and using assumptions about leaf angle distribution in order to estimate the extinction coefficient (Campbell, 1986).

Similarly, to the other vegetation variables, phenological metrics were extracted from time series of fAPAR, CC and LAI according Eq. (1).

2.3.6 GPP data

Half-hourly GPP data were estimated for the ecosystem from net-carbon flux measurements acquired by an eddy covariance system. Details of instrumentation and processing are provided in Delpierre et al. (2016) and are found at (last access: 17 May 2021). GPP was aggregated daily and used to create continuous time series from 2006–2018. The extraction of phenological markers was done according the same procedure (using Eq. 1).

2.4 Statistical analysis

The performance of each of the indirect methods presented above was evaluated with respect to the field phenological observations using three criteria, which are (1) the coefficient of determination (R2) calculated from a simple linear regression between estimated (Pi) and observed dates (Oi) for the different years (N), (2) the mean bias error (MBE), and (3) the mean absolute deviation (MAD) calculated as follows:


Figure 2Illustration of 1-year (2015) time series of OBS (a), NDVI (b), NDVIbr (c), GCC (d), 1−fRvis (e), fAPAR (f), CC (g), LAI (h) and GPP (i) in the Fontainebleau–Barbeau forest. Data are shown as empty circles. The red bold continuous curve is the ADS function (Eq. 1) fitted to the time series. For visual observations, data shown are as a percentage of open buds in spring and as a percentage of non-senescent leaves (100 % – observed percentage of senescent leaves) in autumn. The percentage of open buds is forced to 100 % for the summer growing season and to 0 % during the winter dormancy season. Vertical lines: spring and autumn phenology estimates using MOS and MOF (black) and observed dates (BB-OBS and LS-OBS) (blue).


3 Results

An illustration of the time series of the vegetation variables used (OBS, NDVI, NDVIbr, GCC, (1−fRvis), fAPAR, CC, LAI and GPP) is provided in Fig. 2. Time series of all years (2006–2018) are given in the Supplement Fig. S1.

Time series in Fig. 2 for the year 2015 and in Fig. S1 for all years show that the general patterns of phenological transitions corresponding to the onset of leaves in the spring and to leaf senescence in the autumn are reproduced by all indirect methods but with a variable bias in comparison with the field observation. However, in the autumn, GPP time series show a decline that appears very early in the year, practically from the beginning of summer. GCC time series may also show atypical interannual patterns, with some years during which a GCC decline is also observed very early in the year (2014, 2016–2018 in Fig. S1), although it is slower than the one observed on GPP.

Average phenological dates observed (BB-OBS and LS-OBS) and estimated from the different methods using MOS and MOF markers are given in Fig. 3. All phenological dates, using the six phenological markers (SOS, MOS, EOS, SOF, MOF, EOF), are given Table S2.

Figure 3Average phenological dates in spring (a) and autumn (b) using MOS and MOF phenological markers, respectively, and for the different years.


In spring, field phenological observations (BB-OBS) are earlier than the estimates provided by the majority of the indirect methods (Fig. 3a). However, whatever the method used, the interannual phenological variations are well reproduced. During the autumn, phenological observations (LS-OBS) are later than the indirect methods, except for CC and fAPAR (Fig. 3b), and the performance of the different methods seems more limited compared to spring phenology. Figure 4 shows R2 , MBE and MAD between observed and estimated phenological dates using MOS (Fig. 4a) and MOF (Fig. 4b) markers during spring and autumnal phenological transitions, respectively.

Figure 4Coefficient of determination (R2) (a, b), mean bias error (MBE) (c, d) and mean absolute deviation (MAD) (e, f) in days between observed and estimated phenological dates using MOS and MOF markers during spring (a, c, e) and autumnal (b, d, f) phenological stages. The significance levels of R2 are given by stars: * P<0.05, ** P<0.01, *** P<0.001 and **** P<0.0001. The height of grey boxes marks the average of the statistics across study years (individual years are represented by the black dots). Red horizontal lines represent temporal-resolution-related uncertainties associated with field phenological observations of 3.5 d during the spring and of 7 d during the autumn.


In the spring, R2 values between observed (BB-OBS) and estimated phenological dates (Fig. 4a) based on the MOS marker are all statistically significant (at a significance level of 0.05) and range from about 0.99 to 0.34. All indirect methods are also consistent with each other as shown by the high correlation coefficients in Fig. S3, which confirms the good reproducibility of interannual phenological variability by the different indirect methods. In comparison to BB-OBS, the best correlation is found with GCC over the period 2012–2018 during which RGB images are available (R2=0.99). NDVI and CC are also highly correlated with BB-OBS (R2∼0.89 and 0.80, respectively). Lower but significant correlations are found between BB-OBS and fAPAR, LAI, NDVIbr and 1−fRvis (R2 between 0.6 and 0.7), and the lowest correlation is found between BB-OBS and GPP (R2∼0.34).

Between the different indirect methods and during the spring, R2 between MOS estimates ranges from 0.26 to 0.96 (see correlation matrix in Fig. S3). Best correlations are found between fAPAR, NDVI, NDVIbr, LAI and fRvis(R2>0.89). Good correlations are also found between GCC, NDVI and CC (R2=0.8). Finally, we may also note good consistency between derived dates from GPP- and radiation-based methods (NDVIbr, fAPAR, LAI and fRvis; R2>0.64). The lowest correlation is found between GCC and GPP.

For budburst phenological timings, MBE between BB-OBS and MOS (Fig. 4c) is negative for GCC and CC (estimated date is earlier than observed date). MBE is about 1 d with GCC (MAD  1 d) over 2012–2018 and is also about 1 d with CC over 2006–2018 (MAD  2 d). We note that MBE and MAD (Fig. 4c and e) for these two methods are slightly less than the observation uncertainty of 3.5 d. For the other methods (NDVI, NDVIbr, LAI, fRvis, fAPAR and GPP) MBE and MAD are equal, meaning that MOS estimates from these methods always overestimate the observed phenological dates (BB-OBS). MBE (or MAD) is 3.5 d with NDVI, 6 d with fAPAR and 8 d with NDVIbr. MBE is high with LAI (10 d), fRvis (14 d) and GPP (17 d). Note that for CC, an MBE of about 1 d was obtained after cosine correction of the transmitted PAR according to Eq. (7). Without this correction, MBE increases from 1 d (MAD  2 d) to 6 d (MAD  6 d) and R2 decreases from 0.80 to 0.71. Comparison of the phenological patterns of CC time series obtained with and without cosine correction shows that the cosine correction has the effect of causing an earlier spring phenological start, thus advancing the date of the inflection point (Fig. S4).

During the autumn (Fig. 4b), interannual variation in LS-OBS is well reproduced by CC and NDVI time series which provide estimates that are significantly correlated with the observations (R2=0.80 and 0.63 for CC and NDVI, respectively). Between the indirect methods (Fig. S3), best correlations are found between NDVIbr, fAPAR, NDVI and fRvis (R2∼0.7), LAI and fRvis (R2= 0.58), NDVI and fAPAR (R2= 0.56), NDVI and CC (R2= 0.55), fRvis and fAPAR (R2= 0.55), and CC and fAPAR (R2= 0.42). Surprisingly, correlations between estimated dates from LAI and from CC during the autumn (R2= 0.1), both using the fraction of the transmitted radiation as the unique input, are low compared to what might be expected. Note that it only concerns the senescence stage since the correlation between estimates from LAI and CC during the spring is high (R2 0.74).

During the senescence phase, for the NDVI and CC methods, for which the R2 between estimates and observations are significant, MBE is about 2 d with NDVI (MAD  5 d) and about 14 d with CC (MAD  14 d) (Fig. 4d and f). For CC, MBE decreases from about 37 d without cosine correction to 14 d after correction. The cosine correction yields a faster decrease in CC during the senescence stage (Fig. S4). For CC, LS-OBS are better predicted using thresholds at SOF instead of MOF with an MBE of about 1 d (and MAD of 7 d). MOFs from LAI, fRvis, GCC and GPP provide early estimates compared to LS-OBS. MBE is about 14 d with LAI, 23 d with fRvis, 36 d with GCC and 50 d with GPP. fAPAR leads to estimates that are on average about 30 d later than LS-OBS. Note that for GCC, biases are highly variables between years. For years (2012/2013/2015) for which ADS function do not show the early decline in the autumn estimated dates are very close to OBS (MBE 7 d).

For the phenological markers estimated at the beginning and end of budburst (SOS and EOS) or autumn (SOF and EOF) (Table S2), and considering the period 2015–2017 for which the six phenological markers are available from the intensive sampling, it may be noted that SOS dates are close to observed date (DOY 97) for all methods (between DOY 94–101) except for CC. CC starts to increase earlier, at DOY 82, i.e., 15 d before SOS from OBS. Phenological field observations achieved for understory hornbeam trees over the period 2006–2016 (data not shown) show that, on average, the hornbeam budburst date (i.e., BB-OBS for hornbeam) is around DOY 96 [range 85–107]. MBE between the BB-OBS of hornbeam and SOS estimates is about 1 d (MAD  5 d) for GPP, 5 d (MAD  5 d) for NDVI, 8 d (MAD 8 d) for CC, and between 6–8 d for LAI, fAPAR, NDVIbr and fRvis. For GCC and over 2012–2016, MBE is 2 d. Significant correlations were also obtained between observed hornbeam budburst dates and SOS estimates derived from NDVI, LAI, NDVIbr, CC and fAPAR. R2 ranges between 0.73 and 0.49, and the best correlation is obtained with NDVI-based SOS estimates. Note also that there is a significant correlation between the observed budburst dates of oak and hornbeam (R2∼0.6), but on average hornbeam trees break buds about 10 d earlier than oaks.

For the end of spring, EOSs based on GCC are quite close to EOSs determined from field phenological observations (3 d earlier for GCC). For the other methods, estimated EOSs are later than observed EOS dates. MBEs are 3 d for NDVI, 8 d for fAPAR, 10 d for CC, 14 d for NDVIbr, 20 d for LAI, 28 d for fRvis and 41 d for GPP. During the senescence phase, SOF from NDVI and CC gives the best agreement with the observed SOF date (3 d on average over 2015–2017), followed by fAPAR (6 d). The observed EOF is better predicted using fRvis, CC, NDVI and GPP. MBE is about 3 d for fRvis, 6 d for CC and NDVI, and 9 d for GPP.

As an illustration of the above, Fig. 5 shows average phenological patterns of vegetation variables derived from average parameters of modeled time series through an ADS function fitted to data over the period 2012–2017, common to all vegetation variables, for the spring (Fig. 5a) and the autumn (Fig. 5b) phenological stages. The correspondence between field-observed dates and phenological metrics derived from indirect methods is also shown.

Figure 5Average phenological patterns during budburst (a) and senescence (b) during the period 2012–2017 using modeled time series through an ADS function fitted to the measured time series of NDVI (Normalized Difference Vegetation Index), GCC (greenness chromatic coordinate), broadband NDVI (NDVIbr), LAI (leaf area index), fAPAR, CC (canopy closure), fRvis (fraction of reflected radiation) and GPP (gross primary production). Amplitudes of variations are normalized to 1. Horizontal dotted lines: for each variable, the proportion of the average amplitude that equals the average of the BB-OBS (Fig. 5a) and LS-OBS (Fig. 5b) dates. Horizontal bold red line (y axis = 0.5): mid-amplitude (50 %) corresponding to mid-onset of spring (MOS) and mid-onset of senescence (MOF). Vertical black line: averages of observed phenological dates during 2012–2017 for budburst (BB-OBS) and for senescence (LS-OBS).


Figure 5 illustrates what is described above by showing average temporal patterns during budburst and senescence over the period 2012–2017, common to all eight methods and for which field phenological observations are available in both spring and autumn. Figure 5a shows the good correspondence between the observed dates and the estimates derived from CC and GCC using the mid-amplitude (50 %) MOS threshold. For CC and GCC, MOS clearly marks the budburst date as characterized in the field using the observation protocol used in our study (50 % of trees with at least 50 % open buds per tree crown, BB-OBS). For the NDVI-based method, on average, the mean observed BB-OBS date coincides with the date when NDVI reaches 25 % of its amplitude of variation between the NDVI minimum in winter and NDVI maximum at the end of spring. For the other methods including fAPAR, NDVIbr, LAI, fRvis and GPP, estimated dates at the mid-amplitude threshold are later than BB-OBS with an MAD ranging from 6 to 17 d. A threshold at 20 % of the spring amplitude for GPP, fRvis and NDVIbr and at 10 % for LAI and fAPAR provide estimates with a bias <2 d. During the leaf senescence phase (Fig. 5b), NDVI at mid-amplitude and CC time series and CC when it begins to decrease (∼95 % of its amplitude) provide estimates consistent with the observations. For the other methods, the thresholds shown in Fig. 5b are only valid on average over the period 2012–2017 since the relationships between observations and estimates are not statistically significant as shown in Fig. 4b.

Figure 5a and b also show that the different methods perform relatively well in the spring but deviate from each other in the autumn. Figure S5 shows that the relationships between the different variables are dependent on the considered phenological stage. This is clearly the case in the relationships between fAPAR, NDVI, GCC, GPP and 1−fRvis. It may be noted that the same NDVI value corresponds to a lower fAPAR in spring than in autumn. In other words, NDVI and fAPAR responses to changes in canopy properties follow two different trajectories depending on the season. This “hysteresis” phenomenon may explain the shift between NDVI- and fAPAR-based estimates during the senescence phase (overestimation of the senescence date by the fAPAR), while both predict very close dates during the spring. This phenomenon of hysteresis is also observed in the same way between fAPAR and GCC and fAPAR and GPP. A given GPP or GCC value corresponds to a lower fAPAR in spring than in autumn. We may also note that the relationships between NDVI and GCC are different depending on the season, but a higher GCC corresponds to the same NDVI in spring than in autumn.

4 Discussion

4.1 Ability of GCC to detect phenological transitions

Using RGB-based GCC time series, the MAD with BB-OBS is about 1 d over the 7 years of comparison (2012–2018). This result is in line with previous studies, particularly the study of Richardson et al. (2018b), who compared RGB-camera-based estimates to independent human-eye observations achieved over deciduous forests. They observed average biases of less than 7 d depending on the site, and the best agreement was obtained using GCC at 25 % of its amplitude as a threshold. Many other studies comparing GCC and indirect visual phenological estimates from the same photographs (Klosterman et al., 2014; Wingate et al., 2015) have also concluded that the GCC method yields estimations of the spring phenological date with an average bias of around 7–8 d. In our study, we show that over the 7-year period (Fig. 5a), GCC at the inflection point (MOS) in spring, which corresponds to about 50 % of its annual amplitude derived from modeled time series, is the best predictor of the human-eye-observed BB-OBS dates, which correspond to 50 % of sampled oak trees having at least 50 % open buds (in fact corresponding to about 50 % open buds at the population scale, Nicolas Delpierre unpublished results). This result supports the fact that the camera accurately reports what is observed by the human eye in the field during the spring and that the GCC index is a very good indicator of the timing of budburst. It may also be noted that the phenological field observations have been carried out by the same (three) intercalibrated observers over the study period and according to a constant protocol. This may also participate in explaining the good agreement between field observations and estimated dates from RGB-based GCC index time series. Indeed, several studies have highlighted the importance of uncertainties associated with observations due to various sources, especially the observer effect (Schaber, 2002), and the availability of good-quality data is a prerequisite for a rigorous evaluation of the various indirect methods.

On the other hand, the ability of GCC to estimate the senescence date is variable. For some years, the decline in GCC may start earlier than expected, and therefore estimated dates are strongly biased. When the senescence phase causes pronounced contrasts on RGB images between the summer growth and senescence phases, estimated dates agree with field observations, as for the years 2012, 2013 and 2015. For these years, estimated dates are very close to OBS with an MAD of about 7 d, of the same order of magnitude as the field observation uncertainty. Therefore, during autumn, data quality and data processing appear crucial to obtain reliable estimates, and the extracting of senescence dates based on ADS model may not be the right approach. Other approaches, particularly the spline-based method used for PhenoCam data that has shown good performance (Richardson et al., 2018b), deserve to be employed. Other RGB-based spectral indices using the red band, designed specifically to monitor the autumn phenological transition, such as RCC (red chromatic coordinate) (Klosterman et al., 2014; Liu et al., 2020) or GRVI (green–red vegetation index) (Motohka et al., 2010; Nagai et al., 2012), should also be evaluated. This is beyond the scope of this study and further methodological development is therefore needed to rigorously assess the real potential of this technique for estimating phenological dates during the senescence stage.

Another point to note, as shown in this study (Fig. 2d) and previously pointed out in several other studies (Sonnentag et al., 2012; Keenan et al., 2014; Klosterman et al., 2014; Petach et al., 2014) is that GCC shows annual spikes during the spring followed by a rapid decline. The annual amplitude of GCC determined from the modeled time series is generally smaller than the actual amplitude. In our study, GCC spikes are reached on day 121 on average over 2012–2018. They are not well captured by the ADS model because they are delayed by about 10 d compared to the end of the spring green-up stage determined from GCC-based EOS phenological marker. GCC spikes are also reached 10 d before LAI reaches its maximum. This result is consistent with Keenan et al. (2014). Based on intensive field measurements at canopy and leaf scales, they observed a time lag of about 2 weeks between the canopy maximum LAI measured by an LAI-2000 plant canopy analyzer and GCC spikes. They concluded that GCC depends on leaf color and saturates faster than measured canopy LAI, which was explained by the oblique viewing angle of the camera, which leads to a higher effective LAI. In the same study, they showed that GCC peaks were reached while the main leaf traits (maximum leaf area, chlorophyll content, leaf mass area) continue their development. Similar results were also reported in Yang et al. (2014) and Liu et al. (2015), who showed that GCC peaks in spring were approximately 20 d earlier than the peak of the total chlorophyll concentration. In our study, on average, GCC spikes almost coincide with maximum fAPAR and CC (EOS), whereas these two variables are based on incoming, reflected and transmitted PAR measurements using hemispherical sensors and therefore are integrative of the whole canopy. This result supports the hypothesis of a combined effect of canopy coloring and closure on GCC spikes. However, and contrary to LAI, which is estimated, in this study, only from incident and transmitted radiation, fAPAR and CC also additionally use reflected radiation. Therefore, they are also sensitive to changes in leaf color and other leaf traits during the spring. This may explain the good correspondence between the timings of GCC spikes and the timings of the maximum of fAPAR and CC.

4.2 Ability of NDVI to detect phenological transitions

Results also show that MOS and MOF of NDVI are good proxies of observed dates with an MAD of about 3–4 d in spring over the whole period 2006–2018 and 5 d in autumn over 2011–2017 period. Estimates based on NDVI are also highly correlated with spring and autumn field phenological observations with an R2 of 0.88 and of 0.62, respectively. This reflects the ability of ground-based NDVI time series to reproduce the interannual variability of phenology at this site (Figs. 3b and 4b). This potential has also been shown in previous studies, in evergreen and deciduous forest ecosystems in France, an evergreen tropical rain forest in French Guyana, an herbaceous savanna in Congo and a succession of three annual crops in Belgium (Soudani et al., 2012; Hmimina et al., 2013).

Good agreement between RGB camera indices and proximal NDVI-based measurements has also been shown in crops (Sakamoto et al., 2012) and in herbaceous species (Anderson et al., 2016). However, NDVI measurements does not show the spikes observed in GCC in late spring and our study shows that NDVI is more stable, less scattered and better representative of the LAI plateau throughout the summer growth phase observed in deciduous forests. Similar conclusions were drawn in Petach et al. (2014). In conclusion, the NDVI sensor using MOS and MOF criteria can be regarded as the best option since it provides reliable estimates for monitoring both spring and autumn phenology. In addition, and as highlighted in Hmimina et al. (2013), in situ NDVI measurements using proximal sensors are done a few meters above the top of the canopy, and because NDVI is a normalized index, the effects of the sky conditions produce little noise. Thus, measurements can be carried out under diffuse sky conditions, allowing for the monitoring of vegetation phenology at a high temporal frequency. Nevertheless, proximal NDVI sensors have the disadvantage that measurements remain limited to a narrow field of view and do not allow us to extract key phenological metrics at the individual tree level when this may be possible using RGB camera (Delpierre et al., 2020). The use of multispectral cameras with RGB+NIR bands, which are increasingly used at many sites, may allow us to overcome this inconvenience and should therefore be encouraged.

4.3 Ability of CC to detect phenological transitions

During the spring, a good performance of CC-based method was obtained after cosine correction of the transmitted PAR according to Eq. (7) (Fig. 4a and Table S2). Without this correction, MAD between estimated and observed MOS dates is 3 times larger (6 d vs. 2 d) and R2 slightly lower (0.71 vs. 0.80). It may be noted that uncorrected CC, which corresponds to the complement to 1 of the canopy transmittance, and fAPAR provide similar estimated MOS dates, which are on average about 1 week later than observed dates (Table S2). This result is in line with the study of Perot et al. (2020), conducted in a mature oak forest, which showed that on average estimated MOS dates from canopy transmittance time series are about 7 d later than the observed budburst dates.

A comparison of the phenological patterns of CC time series obtained with and without cosine correction (Fig. S4) shows that the cosine correction has the effect of causing an earlier spring phenological start, thus advancing the date of the inflection point. While the estimated date at the inflection point after cosine correction (CC-MOS) is very close to BB-OBS, SOS appears earlier than the observed SOS of oak trees. This can be explained by the budburst of the first trees of the hornbeam understory, which on average has an earlier budburst date:s about 10 d before the overstory oak trees. During the senescence phase, the cosine correction significantly improved the estimates, but the bias remains high (14 d on average). Despite this bias, autumn CC-MOF dates are the most correlated with observations LS-OBS (R2= 0.8) (Fig. 4b and Table S2). We notice that CC time series are sensitive to the intercepted radiation, which mostly depends on canopy structure and not so much on pigmental (color) properties. Here we derived LS-OBS from the monitoring of the percent of senescent (i.e., colored or fallen leaves) in the canopy, which we build from independent observations of percent colored and percent fallen leaves in the tree crowns. For those years when we continued canopy observations until complete leaf fall, we observed that 50 % leaf fall is typically attained 2–3 weeks after 50 %-senescence, at a date comparable to CC-MOF.

In summary, the cosine correction significantly improves estimated dates based on CC both in the spring and senescence seasons. The new formulation of CC calculation proposed in this study (Eq. 7), which takes into account the effects of seasonal variations in sun angle on the transmitted PAR, merits being tested at other sites in order to assess accurately its performance as it is likely to be dependent on both the canopy structure and the latitude of the site.

4.4 Ability of NDVIbr to detect phenological transitions

The phenological pattern of NDVIbr is comparable to the one obtained from NDVI time series but with greater amplitudes during the spring and autumn phenological transitions for the latter (Figs. 2 and S1). This result is also consistent with Liu et al. (2019), who compared broadband and narrowband NDVI in a temperate broadleaved deciduous forest. Like NDVI, NDVIbr is measured directly above the canopy and seems not to be very sensitive to cloud conditions as also underlined in Wang et al. (2004) and Wilson and Meyers (2007). On average, the deviation between estimated MOS dates from NDVI and NDVIbr is 5 d in spring and 1 d in autumn, respectively. However, while in spring the estimated MOS dates from NDVI and NDVIbr are highly correlated (R2= 0.87), the correlation is lower in autumn (R2= 0.49) and is non-significant between autumn NDVIbr estimates and observed dates LS-OBS. As a result, NDVI and NDVIbr seem to be decorrelated in autumn, and the performance of NDVIbr time series to describe the interannual variability of phenology is only limited to spring.

4.5 Ability of GPP to detect phenological transitions

On average over an 11-year period (2006–2016), GPP starts increasing (GPP-SOS) on DOY 96, 10 d earlier than overstory oak trees (DOY 106, Fig. 3 and Table S2). The starting date of GPP coincides exactly with the date of hornbeam budburst (DOY 96) and of the earliest oaks (Delpierre et al., 2020). However, GPP reaches a maximum in a time interval close to the summer solstice (Figs. 2 and 5a) and then starts to decline immediately after. Consequently, GPP-MOS overestimates BB-OBS by about 17 d. This result is in line with other previous studies which have shown that GPP peaks several weeks later than the other variables. Toomey et al. (2015) showed that the start of GPP in spring coincides with the onset of GCC, but GPP peaks 2–4 weeks later. They also noted an immediate decline in GPP once its peak is reached. Similar conclusions between GCC and GPP can also be drawn from Richardson et al. (2009).

During the autumn phase, based on ADS function, the GPP time series fails to produce plausible estimates of LS-OBS, either using SOF, MOF or EOF criteria.

As underlined above, among all the indirect methods evaluated in this study, estimates of budburst dates derived from GPP time series using the MOS criterion based on ADS are the most biased estimates and are also the least correlated with the observed phenological dates of oak trees (MBE 17 d, R2= 0.34, Fig. 4a). This weak correlation can be explained both by a starting of the GPP simultaneously with the budburst of the hornbeam understory and the high dependency of GPP on the solar radiation level, in addition to the effects of the increase in the LAI and the leaf maturation (Delpierre et al., 2009a). Figures 2 and 5a show that GPP reaches a short-lived plateau around the summer solstice in June, when both maximum LAI is reached and solar irradiance is at its maximum. On the other hand, MOF dates during the autumn are earlier than LS-OBS (Figs. 2, 5 and Table S2). Consequently, the length of the period of budburst and leaf development in spring between GPP-derived SOS and EOS dates, is about 57 d over the 13 years of measurements, while it is only about 17 d from in situ NDVI. The length of the growing season, between estimated dates of MOS and MOF, is also greatly reduced, and it is only 130 d based on GPP, whereas it is 192 d from NDVI and 199 d from field phenological observations. Similar results are shown in the studies of Lu et al. (2018) and Keenan et al. (2014). In conclusion, the extraction of phenology from GPP time series using inflection points of transitions in the spring and autumn is therefore not representative of the canopy leaf display, and other approaches based on absolute or relative thresholds of GPP as in Richardson et al. (2010) and in Wu et al. (2017) may be more representative. Nevertheless, GPP remains a composite signal driven by changes in vegetation phenology and physiological processes that are under the control of the fluctuations of abiotic factors, and its use to derive the timings of phenological events must be carried out with great care, as strongly emphasized in Gonsamo et al. (2013).

4.6 Hysteresis phenomena between vegetation variables according to the spring and senescence seasons

As shown in Fig. 5, the performance of the different methods for estimating key phenological dates differs between spring and autumn. While the correlations between estimates and observations are all significant during spring (Fig. 4a), only NDVI and CC provide estimates consistent with autumn observations (Fig. 4b). The hysteresis phenomenon that characterizes some relationships between the vegetation variables used in the different methods reflects their different biophysical meanings (Fig. S5). This is particularly the case for the relationships between NDVI and fAPAR and between GCC and fAPAR. In spring, the performances of NDVI and fAPAR are similar, whereas in autumn the fPAR provides very late estimates. This can be explained by a high sensitivity of NDVI and GCC to pigment changes during senescence, whereas fAPAR responds mainly to leaf fall and canopy opening.

4.7 Linking phenological dates recorded by field phenologists and phenological metrics predicted by indirect proximal methods

The analysis of the link between phenological dates based on field observation and those derived from modeled time series (Fig. 5a and b) shows that, on average over 13 years, BB-OBSs (corresponding approximately to 50 % open buds in the canopy) are better predicted by MOS (50 % of the annual amplitude of variation) for methods based on GCC and CC. For the NDVI-based method, a threshold of 25 % of its amplitude coincides with the average observed date. However, due to the rapid increase in NDVI during the spring, a 50 % threshold also provides estimates with a bias of the same order of magnitude as the uncertainty in the phenological observations (3.5 d). For the other methods (GPP, fRvis, NDVIbr, fAPAR and LAI), a threshold at 20 % of the annual amplitude appears more appropriate to estimate the average observed date of budburst. During the senescence phase, and for NDVI- and CC-based methods, for which observations and estimates are significantly correlated, the MOF of NDVI is very close to the observed LS-OBS date (50 % of trees having at least 50 % of senescent or fallen leaves per tree crown) and SOF of CC is more in line with the observed date but less stable than MOF.

Although they are based on data acquired over a long period covering 13 years of measurements and observations, these thresholds may be specific to our study site and their stability and genericity merit further study in other forest ecosystems.

4.8 Summary remarks on deriving phenological metrics from radiation-based methods at EC flux-tower sites

Many flux-tower sites that use the eddy covariance technique routinely acquire the biometeorological variables used in the calculation of GPP, LAI, fRvis, NDVIbr, fAPAR and CC. During the spring stage, LAI, fRvis and GPP-based estimates are biased by about 10 to 17 d. fRvis and GPP are the worst-performing predictors, especially GPP. On the other hand, this study shows that NDVIbr, fAPAR and CC are able to reproduce interannual variation in spring budburst with a bias of about 1 week when MOS is considered (Figs. 3 and 4, Table S2). In the same vein, the use of the CC-based method is also another robust alternative for monitoring spring and autumn phenological transitions at EC flux-tower sites. However, CC and fAPAR require additional measurements of transmitted radiation below the canopy. Indeed, such measurements are not commonly achieved at EC flux measurement sites and should be deployed as transmitted radiation data time series, in addition to phenology, can also be used to estimate leaf area index and to characterize its seasonal dynamics (Keenan et al., 2014). These measurements must be performed using an appropriate number of below-canopy radiation sensors to take the heterogeneity of the canopy structure into account (Pontailler, 1990; Link et al., 2004; Garrity et al., 2011; Webster et al., 2016). When such data are available, derived phenological metrics can be used to conduct retrospective studies in order to interpret the interannual variability of carbon fluxes and are preferable to those derived from the fluxes themselves such GPP or NEE as already pointed out in Gonsamo et al. (2013).

5 Conclusions

We used various methods to characterize the temporal dynamics of forest canopy in a temperate deciduous forest. Field phenological observations provided exhaustive multi-year samples allowing us to accurately assess the potential of each method. However, we emphasize that this potential remains relative because it was evaluated using the ADS method applied to all vegetation proxies and regarded in this study as the only method of extracting phenological dates in order not to bias their comparison. Using the ADS-based phenology extraction method led to results showing that this potential is different depending on the method and the season. During the spring phase, GCC, NDVI and CC, using the inflection point MOS criterion, provide estimates closest to observed dates with an absolute bias of less than 4 d and of the same order as the temporal resolution of phenological observations (3.5 d). For CC, this is obtained only after a cosine correction of the transmitted PAR, a correction that takes the variation in the optical path in the canopy due to the seasonal variation in the solar angle into account. Without this correction, the prediction bias increases from about 2 to 6 d. Using the MOS criterion, NDVIbr and fAPAR also give satisfactory estimates with a bias of around 1 week that corresponds to the temporal resolution generally used in phenological observations. However, for these variables as well as for fRvis, LAI and GPP, a threshold of 20 % of their transition amplitude in spring allows us to obtain more precise estimates in agreement with observed dates. During the senescence phase, only MOF of NDVI and CC can reproduce the interannual variability of leaf senescence. However, these findings are specific to the ADS-based method used to derive phenological markers from time series data. More appropriate methods, especially for GPP and GCC time series, could have provided better estimates of senescence date.

This study validated the estimates provided by the different methods by comparing them with phenological observations carried out according to the same protocol by intercalibrated observers and over 13 years of field observations for budburst and 7 years for leaf senescence. But more particularly, this study demonstrated the good performance of methods based on broadband NDVI (NDVIbr), the fraction of absorbed PAR (fAPAR) and canopy closure (CC) that use solar radiation data routinely recorded at several EC flux-tower sites. This opens real perspectives to conduct retrospective studies for a better interpretation of the interannual variation in carbon fluxes. fAPAR and CC use transmitted radiation measurements below the canopy which are less common but merit being widely deployed at EC flux measurement sites.

Data availability

The data used for this study and the computer codes (in MATLAB and R) are available on request from the corresponding author.


The supplement related to this article is available online at:

Author contributions

KS designed the study. He did the data analysis and wrote the paper with contributions from all co-authors. ND, DB, JYP, LS, GV and ED participated in the data collection and data formatting.

Competing interests

The authors declare that they have no conflict of interest.


Many thanks to all colleagues who participated in the installation of the various instruments at the Fontainebleau–Barbeau site and all those involved in the data collection used in this study. The FR-Fon study site has been funded through several French and European research framework programs (GIP Ecofor, Allenvi, CarboEurope, FP6; CarboExtreme, FP7). It is part of the Integrated Carbon Observation System (ICOS, FP7) European research infrastructure and of the SOERE-Ecofor French research network.

Review statement

This paper was edited by Trevor Keenan and reviewed by two anonymous referees.


Anderson, H. B., Nilsen, L., Tømmervik, H., Karlsen, S. R., Nagai, S., and Cooper, E. J.: Using Ordinary Digital Cameras in Place of Near-Infrared Sensors to Derive Vegetation Indices for Phenology Studies of High Arctic Vegetation, Remote Sens., 8, 847,, 2016. 

Badeck, F. W., Bondeau, A., Böttcher, K., Doktor, D., Lucht, W., Schaber, J., and Sitch, S.: Responses of spring phenology to climate change, New Phytol., 162, 295–309,, 2004. 

Baldocchi, D. D., Matt, D. R., Hutchison, B. A., and McMillen, R. T.: Solar radiation within an oak-hickory forest: an evaluation of extinction coefficients for several radiation components during fully leafed and leafless periods, Agr. Forest Meteorol., 32, 307–322,, 1984. 

Baldocchi, D. D., Falge, E., Gu, L., Olson, R., Hollinger, D., Running, S., Anthoni, P., Bernhofer, C., Davis, K., Fuentes, J., Goldstein, A., Katul, G., Law, B., Lee, X., Malhi, Y., Meyers, T., Munger, J. W., Oechel, W., Pilegaard, K., Schmid, H. P., Valentini, R., Verma, S., Vesala, T., Wilson, K., and Wofsy, S.: FLUXNET: a new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor and energy flux densities, B. Am. Meteorol. Soc., 82, 2415–2435,<2415:FANTTS>2.3.CO;2, 2001. 

Campbell, G. S.: Extinction coefficients for radiation in plant canopies calculated using an ellipsoidal inclination angle distribution, Agr. Forest Meteorol., 36, 317–321,, 1986. 

Campbell, G. S. and Norman, J. M.: The Light Environment of Plant Canopies, in: An Introduction to Environmental Biophysics, Springer, New York, USA, 1998. 

Delpierre, N., Berveiller, D., Granda, E., and Dufrêne, E.: Wood phenology, not carbon input, controls the interannual variability of wood growth in a temperate oak forest, New Phytol., 210, 459–470,, 2016. 

Delpierre, N., Soudani, K., François, C., Köstner, B., Pontailler, J., Nikinmaa, E., Misson, L., Aubinet, M., Bernhofer, C., Granier, A., Grünwald, T., Heinesch, B., Longdoz, B., Ourcival, J., Rambal, S., Vesala, T., and Dufrene, E.: Exceptional carbon uptake in European forests during the warm spring of 2007: a data-model analysis, Global Change Biol., 15, 1455–1474,, 2009a. 

Delpierre, N., Dufrêne, E., Soudani, K., Ulrich, E., Cecchini, S., Boé, J., and François, C.: Modelling interannual and spatial variability of leaf senescence for three deciduous tree species in France, Agr. Forest Meteorol., 149, 938–948,, 2009b. 

Delpierre, N., Soudani, K., Berveiller, D., Dufrêne, E., Hmimina, G., and Vincent, G.: “Green pointillism”: detecting the within-population variability of budburst in temperate deciduous trees with phenological cameras, Int. J. Biometeorol., 64, 663–670,, 2020. 

Denéchère, R., Delpierre, N., Apostol, E., Berveiller, D., Bonne, F., Cole, E., Delzon, S., Dufrêne, E., Gressler, E., Jean, F., Lebourgeois, F., Liu, G., Louvet, J., Parmentier, J., Soudani, K., and Vincent, G.: The within-population variability of leaf spring and autumn phenology is influenced by temperature in temperate deciduous trees, Int. J. Biometeorol., 65, 369–379,, 2019. 

Dragoni, D., Schmid, H. P., Wayson, C. A., Potter, H., Grimmond C. S. B., and Randolph, J. C.: Evidence of increased net ecosystem productivity associated with a longer vegetated season in a deciduous forest in south-central Indiana, USA, Global Change Biol., 17, 886–897,, 2011. 

Eklundh, L., Jin, H., Schubert, P., Guzinski, R., and Heliasz, M.: An optical sensor network for vegetation phenology monitoring and satellite data calibration, Sensors, 11, 7678–7709,, 2011. 

Franz, D., Acosta, M., Altimir, N., Arriga, N., Arrouays, D., Aubinet, M., Aurela, M., Ayres, E., López-Ballesteros, A., Barbaste, M., Berveiller, D., Biraud, S., Boukir, H., Brown, T., Brümmer, C., Buchmann, N., Burba, G., Carrara, A., Cescatti, A., Ceschia, E., Clement, R., Cremonese, E., Crill, P., Darenova, E., Dengel, S., D'Odorico, P., Gianluca, F., Fleck, S., Fratini, G., Fuß, R., Gielen, B., Gogo, S., Grace, J., Graf, A., Grelle, A., Gross, P., Grünwald, T., Haapanala, S., Hehn, M., Heinesch, B., Heiskanen, J., Herbst, M., Herschlein, C., Hörtnagl, L., Hufkens, K., Ibrom, A., Jolivet, C., Joly, L., Jones, M., Kiese, R., Klemedtsson, L., Kljun, N., Klumpp, K., Kolari, P., Kolle, O., Kowalski, A., Kutsch, W., Laurila, T., De Ligne, A., Linder, S., Lindroth, A., Lohila, A., Longdoz, B., Mammarella, I., Manise, T., Marañon-Jimenez, S., Matteucci, G., Mauder, M., Meier, P., Merbold, L., Mereu, S., Metzger, S., Migliavacca, M., Mölder, M., Montagnani, L., Moureaux, C., Nelson, D., Nemitz, E., Nicolini, G., Nilsson, M. B., Op de Beeck, M., Osborne, B., Ottosson Löfvenius, M., Pavelka, M., Peichl, M., Peltola, O., Pihlatie, M., Pitacco, A., Pokorny, R., Pumpanen, J., Ratié, C., Schrumpf, M., Sedlák, P., Serrano Ortiz, P., Siebicke, L., Šigut, L., Silvennoinen, H., Simioni, G., Skiba, U., Sonnentag, O., Soudani, K., Soulé, P., Steinbrecher, R., Tallec, T., Thimonier, A., Tuittila, E., Tuovinen, J., Vestin, P., Vincent, G., Vincke, C., Vitale, D., Waldner, P., Weslien, P., Wingate, L., Wohlfahrt, G., Zahniser, M., and Vesala, T.: Towards long-term standardised carbon and greenhouse gas observations for monitoring Europe's terrestrial ecosystems: a review, Int. Agrophys., 32, 439–455,, 2018. 

Garrity, S. R., Bohrer, G., Maurer, K. D., Mueller, K. L., Vogel, C. S., and Curtis, P. S.: A comparison of multiple phenology data sources for estimating seasonal transitions in deciduous forest carbon exchange, Agr. Forest Meteorol., 151, 1741–1752,, 2011. 

Gielen, B., Acosta, M., Altimir, N., Buchmann, N., Cescatti, A., Ceschia, E., Fleck, S., Hörtnagl, L., Klumpp, K., Kolari, P., Lohila, A., Loustau, D., Marañon-Jimenez, S., Manise, T., Matteucci, G., Merbold, L., Metzger, C., Moureaux, C., Montagnani, L., Nilsson, M. B., Osborne, B., Papale, D., Pavelka, M., Saunders, M., Simioni, G., Soudani, K., Sonnentag, O., Tallec, T., Tuittila, E., Peichl, M., Pokorny, R., Vincke, C., and Wohlfahrt, G.: Ancillary vegetation measurements at ICOS ecosystem stations, Int. Agrophys., 10, 645–664,, 2018. 

Gonsamo, A., Chen, J. M., and D'Odorico, P.: Deriving land surface phenology indicators from CO2 eddy covariance measurements, Ecol. Indic. 29, 203–207,, 2013. 

Goulden, M. L., Munger, J. W., Fan, S. M., Daube, B. C., and Wofsy, S. C.: Exchange of carbon dioxide by a deciduous forest: response to interannual climate variability, Science, 271, 1576–1578,, 1996. 

Hmimina, G., Dufrene, E., Pontailler, J.-Y., Delpierre, N., Aubinet, M., Caquet, B., de Grandcourt, A., Burban, B., Flechard, C., Granier, A., Gross, P., Heinesch, B., Longdoz, B., Moureaux, C., Ourcival, J.-M., Rambal, S., Saint-André, L., and Soudani, K.: Evaluation of the potential of MODIS satellite data to predict vegetation phenology in different biomes: An investigation using ground based NDVI measurements, Remote Sens. Environ., 132, 145–158,, 2013. 

Holst, T., Hauser, S., Kirchgässner, A., Matzarakis, A., Mayer, H., and Schindler, D.: Measuring and modelling plant area index in beech stands, Int. J. Biometeorol., 48, 192–201,, 2004. 

Huemmrich, K. F., Black, T. A., Jarvis, P. G., McCaughey, J. H., and Hall, F. G.: High temporal resolution NDVI phenology from micrometeorological radiation sensors, J. Geophys. Res.-Atmos., 104, 27935–27944,, 1999. 

Jenkins, J. P., Richardson, A. D., Braswell, B. H., Ollinger, S. V., Hollinger, D. Y., and Smith, M.-L.: Refining light-use efficiency calculations for a deciduous forest canopy using simultaneous tower-based carbon flux and radiometric measurements, Agr. Forest Meteorol., 143, 64–79,, 2006. 

Keenan, T. F., Darby, B., Felts, E., Sonnentag, O., Friedl, M. A., Hufkens, K., O'Keefe, J., Klosterman, S., Munger, J. W., Toomey, M., and Richardson, A. D.: Tracking Forest Phenology and Seasonal Physiology Using Digital Repeat Photography: A Critical Assessment, Ecol. Appl., 24, 1478–1489,, 2014. 

Klosterman, S. T., Hufkens, K., Gray, J. M., Melaas, E., Sonnentag, O., Lavine, I., Mitchell, L., Norman, R., Friedl, M. A., and Richardson, A. D.: Evaluating remote sensing of deciduous forest phenology at multiple spatial scales using PhenoCam imagery, Biogeosciences, 11, 4305–4320,, 2014. 

Kobayashi, H., Nagai, S., Kim, Y., Yang, W., Ikeda, K., Ikawa, H., Nagano, H., and Suzuki, R.: In situ observations reveal how spectral reflectance responds to growing season phenology of an open evergreen forest in Alaska, Remote Sens.-Basel, 10, 1071,, 2018. 

Link, T. E., Marks, D., and Hardy, J.: A deterministic method to characterize canopy radiative transfer properties, Hydrol. Process., 18, 3583–3594,, 2004. 

Liu, F., Wang, X., and Wang, C.: Autumn phenology of a temperate deciduous forest: Validation of remote sensing approach with decadal leaf-litterfall measurements, Agr. Forest Meteorol., 15, 107758,, 2019. 

Liu, Y., Wu, C., Sonnentag, O., Desai, A. R., and Wang, J.: Using the red chromatic coordinate to characterize the phenology of forest canopy photosynthesis, Agr. Forest Meteorol., 285–286, 107910,, 2020. 

Liu, Z., Hu, H., Yu, H., Yang, X., Yang, H., Ruan, C., Wang, Y., and Tang, J.: Relationship between leaf physiologic traits and canopy color indices during the leaf expansion period in an oak forest, Ecosphere, 6, 259,, 2015. 

Lu, X., Liu, Z., Zhou, Y., Liu, Y., An, S., and Tang, J.: Comparison of Phenology Estimated from Reflectance-Based Indices and Solar-Induced Chlorophyll Fluorescence (SIF) Observations in a Temperate Forest Using GPP-Based Phenology as the Standard, Remote Sens.-Basel, 10, 932,, 2018. 

McCree, K. J.: Test of current definitions of photosynthetically active radiation against leaf photosynthesis data, Agr. Meteorol., 10, 443–453,, 1972. 

Menzel, A., Sparks, T. H., Estrella, N., Koch, E., Aasa, A., Ahas, R., Alm-Kübler, K., Bissolli, P., Braslavská, O., Briede, A., Chmielewski, F. M., Crepinsek, Z., Curnel, Y., Dahl, Å., Defila, C., Donnelly, A., Filella, Y., Jatczak, K., Måge, F., Mestre, A., Nordli, Ø., Peñuelas, J., Pirinen, P., Remišová, V., Scheifinger, H., Striz, M., Susnik, A., Van Vliet, A. J. H., Wielgolaski, F.-E., Zach, S., and Zust, A.: European phenological response to climate change matches the warming pattern, Global Change Biol., 12, 1969–1976,, 2006. 

Milliman, T., Seyednasrollah, B., Young, A. M., Hufkens, K., Friedl, M. A., Frolking, S., Richardson, A. D., Abraha, M., Allen, D. W., Apple, M., Arain, M. A., Baker, J., Baker, J. M., Bernacchi, C. J., Bhattacharjee, J., Blanken, P., Bosch, D. D., Boughton, R., Boughton, E. H., Brown, R. F., Browning, D. M., Brunsell, N., Burns, S. P., Cavagna, M., Chu, H., Clark, P. E., Conrad, B. J., Cremonese, E., Debinski, D., Desai, A. R., Diaz-Delgado, R., Duchesne, L., Dunn, A. L., Eissenstat, D. M., El-Madany, T., Ellum, D. S. S., Ernest, S. M., Esposito, A., Fenstermaker, L., Flanagan, L. B., Forsythe, B., Gallagher, J., Gianelle, D., Griffis, T., Groffman, P., Gu, L., Guillemot, J., Halpin, M., Hanson, P. J., Hemming, D., Hove, A. A., Humphreys, E. R., Jaimes-Hernandez, A., Jaradat, A. A., Johnson, J., Keel, E., Kelly, V. R., Kirchner, J. W., Kirchner, P. B., Knapp, M., Krassovski, M., Langvall, O., Lanthier, G., Maire, G. I., Magliulo, E., Martin, T. A., McNeil, B., Meyer, G. A., Migliavacca, M., Mohanty, B. P., Moore, C. E., Mudd, R., Munger, J. W., Murrell, Z. E., Nesic, Z., Neufeld, H. S., Oechel, W., Oishi, A. C., Oswald, W. W., Perkins, T. D., Reba, M. L., Rundquist, B., Runkle, B. R., Russell, E. S., Sadler, E. J., Saha, A., Saliendra, N. Z., Schmalbeck, L., Schwartz, M. D., Scott, R. L., Smith, E. M., Sonnentag, O., Stoy, P., Strachan, S., Suvocarev, K., Thom, J. E., Thomas, R. Q., Van den Berg, A. K., Vargas, R., Vogel, C. S., Walker, J. J., Webb, N., Wetzel, P., Weyers, S., Whipple, A. V., Whitham, T. G., Wohlfahrt, G., Wood, J. D., Yang, J., Yang, X., Yenni, G., Zhang, Y., Zhang, Q., Zona, D., Baldocchi, D., and Verfaillie, J.: PhenoCam Dataset v2.0: Digital Camera Imagery from the PhenoCam Network, 2000–2018, ORNL DAAC, Oak Ridge, Tennessee, USA,, 2009. 

Motohka, T., Nasahara, K. N., Oguma, H., and Tsuchida, S.: Applicability of green-red vegetation index for remote sensing of vegetation phenology, Remote Sens.-Basel, 2, 2369–2387,, 2010. 

Nagai, S., Saitoh, T. M., Kobayashi, H., Ishihara, M., Motohka, T., Suzuki, R., Nasahara, K. N., and Muraoka, H.: In situ examination for the relationship between various vegetation indices and tree phenology in an evergreen coniferous forest, Japan, Int. J. Remote Sens., 33, 6202–6214,, 2012. 

Perot, T., Balandier, P., Couteau, C., Perret, S., Seigner, V., and Korboulewsky, N.: Transmitted light as a tool to monitor tree leaf phenology and development applied to Quercus petraea, Agr. Forest Meteorol., 275, 37–46,, 2019. 

Petach, A. R., Toomey, M., Aubrecht, D. M., and Richardson, A. D.: Monitoring vegetation phenology using an infrared-enabled security camera, Agr. Forest Meteorol., 195–196, 143–151,, 2014. 

Piao, S., Liu, Q., Chen, A., Janssens, I., Fu, Y., Dai, J., Liu, L., Lian, X., Shen, M., and Zhu, X.: Plant phenology and global climate change: current progresses and challenges, Global Change Biol., 25, 1922–1940,, 2019. 

Pontailler, J.-Y.: A Cheap Quantum Sensor Using a Gallium Arsenide Photodiode, Funct. Ecol., 4, 591–596,, 1990. 

Pontailler, J.-Y., Hymus, G. J., and Drake, B. G.: Estimation of leaf area index using ground-based remote sensed NDVI measurements: Validation and comparison with two indirect techniques, Can. J. Remote Sens., 29, 381–387,, 2003. 

Richardson, A. D.: Tracking seasonal rhythms of plants in diverse ecosystems with digital camera imagery, New Phytol., 222, 1742–1750,, 2019. 

Richardson, A. D., Jenkins, J. P., Braswell, B. H., Hollinger, D. Y., Ollinger, S. V., and Smith, M. L.: Use of digital webcam images to track spring green-up in a deciduous broadleaf forest, Oecologia, 152, 323–334,, 2007. 

Richardson, A. D., Black, T. A., Ciais, P., Delbart, N., Friedl, M. A., Gobron, N., Hollinger, D. Y., Kutsch, W. L., Longdoz, B., Luyssaert, S., Migliavacca, M., Montagnani, L., Munger, J. W., Moors, E., Piao, S., Rebmann, C., Reichstein, M., Saigusa, N., Tomelleri, E., Vargas, R., and Varlagin, A.: Influence of spring and autumn phenological transitions on forest ecosystem productivity, Philos. T. Roy. Soc. B, 365, 3227–3246,, 2010. 

Richardson, A. D., Hufkens, K., Milliman, T., and Frolking, S.: Intercomparison of phenological transition dates derived from the PhenoCam Dataset V1.0 and MODIS satellite remote sensing, Sci. Rep.-UK, 8, 5679,, 2018a. 

Richardson, A. D., Hufkens, K., Milliman, T., Aubrecht, D. M., Chen, M., Gray, J. M., Johnston, M. R., Keenan, T. F., Klosterman, S. T., Kosmala, M., Melaas, E. K., Friedl, M. A., and Frolking, S.: Tracking vegetation phenology across diverse North American biomes using PhenoCam imagery, Scientific Data, 5, 180028,, 2018b. 

Roetzer, T., Wittenzeller, M., Haeckel, H., and Nekovar, J.: Phenology in central Europe: difference and trends of spring phenophases in urban and rural areas, Int. J. Biometeorol., 44, 60–66,, 2000. 

Ryu, Y., Baldocchi, D. D., Verfaillie, J., Ma, S., Falk, M., Ruiz-Mercado, I., Hehn, T., and Sonnentag, O.: Testing the performance of a novel spectral reflectance sensor, built with light emitting diodes (LEDs), to monitor ecosystem metabolism, structure and function, Agr. Forest Meteorol., 150, 1597–1606,, 2010. 

Sakamoto, T., Gitelson, A. A., Nguy-Robertson, A. L., Arkebauer, T. J., Wardlow, B. D., Suyker, A. E., Verma, S. B., and Shibayama, M.: An alternative method using digital cameras for continuous monitoring of crop status, Agr. Forest Meteorol., 154–155, 113–126,, 2012. 

Schaber, J.: Phenology in Germany in the 20th century: methods, analyses and models, University of Potsdam, Germany, 2002. 

Schaber, J. and Badeck, F. W.: Evaluation of methods for the combination of phenological time series and outlier detection, Tree Physiol., 22, 973–982,, 2002. 

Sims, D. A. and Gamon, J. A.: Relationship between leaf pigment content and spectral reflectance across a wide range species, leaf structures and development stages, Remote Sens. Environ., 81, 337–354,, 2002. 

Sonnentag, O., Hufkens, K., Teshera-Sterne, C., Young, A. M., Friedl, M., Braswell, B. H., Milliman, T., O'Keefe, J., and Richardson, A. D.: Digital repeat photography for phenological research in forest ecosystems, Agr. Forest Meteorol., 152, 159–177,, 2012. 

Soudani, K., le Maire, G., Dufrêne, E., François, C., Delpierre, N., Ulrich, E., and Cecchini, S.: Evaluation of the onset of green-up in temperate deciduous broadleaf forests derived from Moderate Resolution Imaging Spectroradiometer (MODIS) data, Remote Sens. Environ., 12, 2643–2655,, 2008. 

Soudani, K., Hmimina, G., Delpierre, N., Pontailler, J.-Y., Aubinet, M., Bonal, D., Caquet, B., de Grandcourt, A., Burban, B., Flechard, C., Guyon, D., Granier, A., Gross, P., Heinesh, B., Longdoz, B., Loustau, D., Moureaux, C., Ourcival, J.-M., Rambal, S., Saint André, L., and Dufrêne, E.: Ground-based Network of NDVI measurements for tracking temporal dynamics of canopy structure and vegetation phenology in different biomes, Remote Sens. Environ., 123, 234–245,, 2012. 

Sparks, T. H. and Carey, P. D.: The responses of species to climate over two centuries: an analysis of the Marsham phenological record, 1736–1947, J. Ecol., 83, 321–329,, 1995. 

Toda, M. and Richardson, A. D.: Estimation of plant area index and phenological transition dates from digital repeat photography and radiometric approaches in a hardwood forest in the Northeastern United States, Agr. Forest Meteorol., 249, 457–466,, 2018. 

Toomey, M., Friedl, M. A., Frolking, S., Hufkens, K., Klosterman, S., Sonnentag, O., Baldocchi, D. D., Bernacchi, C. J., Biraud, S. C., Bohrer, G., Brzostek, E., Burns, S. P., Coursolle, C., Hollinger, D. Y., Margolis, H. A., McCaughey, H., Monson, R. K., Munger, J. W., Pallardy, S., Phillips, R. P., Torn, M. S., Wharton, S., Zeri, M., and Richardson, A. D.: Greenness indices from digital cameras predict the timing and seasonal dynamics of canopy-scale photosynthesis, Ecol. Appl., 25, 99-115,, 2015. 

Wang, Q., Tenhunen, J., Dinh, N. Q., Reichstein, M., Vesala, T., and Keronen, P.: Similarities in ground- and satellite-based NDVI time-series and their relationship to physiological activity of a Scots pine forest in Finland, Remote Sens. Environ., 93, 225–237,, 2004. 

Wang, Q., Tenhunen, J., Schmidt, M., Kolcun, O., Droesler, M., and Reichstein, M.: Estimation of total, direct and diffuse PAR under clear skies in complex alpine terrain of the National Park Berchtesgaden, Germany, Ecol. Model., 196, 149–162,, 2006. 

Webster, C., Rutter, N., Zahner, F., and Jonas, T.: Measurement of incoming radiation below forest canopies: A comparison of different radiometer configurations, J. Hydrometeorol., 17, 853–864,, 2016. 

Wilson, T. B. and Meyers, T. B.: Determining vegetation indices from solar and photosynthetically active radiation fluxes, Agr. Forest Meteorol., 144, 160–179,, 2007. 

Wingate, L., Ogée, J., Cremonese, E., Filippa, G., Mizunuma, T., Migliavacca, M., Moisy, C., Wilkinson, M., Moureaux, C., Wohlfahrt, G., Hammerle, A., Hörtnagl, L., Gimeno, C., Porcar-Castell, A., Galvagno, M., Nakaji, T., Morison, J., Kolle, O., Knohl, A., Kutsch, W., Kolari, P., Nikinmaa, E., Ibrom, A., Gielen, B., Eugster, W., Balzarolo, M., Papale, D., Klumpp, K., Köstner, B., Grünwald, T., Joffre, R., Ourcival, J.-M., Hellstrom, M., Lindroth, A., George, C., Longdoz, B., Genty, B., Levula, J., Heinesch, B., Sprintsin, M., Yakir, D., Manise, T., Guyon, D., Ahrends, H., Plaza-Aguilar, A., Guan, J. H., and Grace, J.: Interpreting canopy development and physiology using a European phenology camera network at flux sites, Biogeosciences, 12, 5995–6015,, 2015.  

Wohlfahrt, G., Pilloni, S., Hörtnagl, L., and Hammerle, A.: Estimating carbon dioxide fluxes from temperate mountain grasslands using broad-band vegetation indices, Biogeosciences, 7, 683–694,, 2010. 

Wu, C., Peng, D., Soudani, K., Siebicke, L., Gough, C. M., Arain, M. A., Bohrer, G., Lafleur, P. M., Peichl, M., Gonsamo, A., Shiguang, X., Fang, B., and Quansheng, G.: Land surface phenology derived from normalized difference vegetation index (NDVI) at global FLUXNET sites, Agr. Forest Meteorol., 233, 171–182,, 2017. 

Yang, X., Tang, J., and Mustard, J. F.: Beyond leaf color: Comparing camera-based phenological metrics with leaf biochemical, biophysical, and spectral properties throughout the growing season of a temperate deciduous forest, J. Geophys. Res.-Biogeo., 119, 181–191,, 2014. 

Zhang, X., Friedl, M. A., Schaaf, C. B., Strahler, A. H., Hodges, J. C. F., Gao, F., Reed, B. C., and Huete, A.: Monitoring vegetation phenology using MODIS, Remote Sens. Environ., 84, 471–475,, 2003. 

Short summary
We present an exhaustive comparative survey of eight proximal methods to estimate forest phenology. We focused on methodological aspects and thoroughly assessed deviations between predicted and observed phenological dates and pointed out their main causes. We show that proximal methods provide robust phenological metrics. They can be used to retrieve long-term phenological series at flux measurement sites and help interpret the interannual variability and trends of mass and energy exchanges.
Final-revised paper