Interpreting canopy development and physiology using a European phenology camera network at flux sites

. Plant phenological development is orchestrated through subtle changes in photoperiod, temperature, soil moisture and nutrient availability. Presently, the exact timing of plant development stages and their response to climate and management practices are crudely represented in land surface models. As visual observations of phenology are la-borious, there is a need to supplement long-term observations with automated techniques such as those provided by digital repeat photography at high temporal and spatial resolution. We present the ﬁrst synthesis from a growing observational network of digital cameras installed on towers across Europe above deciduous and evergreen forests, grasslands and crop-lands, where vegetation and atmosphere CO 2 ﬂuxes are measured continuously. Using colour indices from digital images and using piecewise regression analysis of time series, we explored whether key changes in canopy phenology could be detected automatically across different land use types in the network. The piecewise regression approach could capture the start and end of the growing season, in addition to identifying striking changes in colour signals caused by ﬂowering and management practices such as mowing. Exploring the dates of green-up and senescence of deciduous forests extracted by the piecewise regression approach against dates estimated from visual observations, we found that these phenological events could be detected adequately (RMSE < 8 and 11 days for leaf out and leaf fall, respectively). We also investigated whether the seasonal patterns of red, green and blue colour fractions derived from digital images could be modelled mechanistically using the PROSAIL model pa-rameterised with information of seasonal changes in canopy leaf area and leaf chlorophyll and carotenoid concentrations. From a model sensitivity analysis we found that variations in colour fractions, and in particular the late spring ‘green hump’ observed repeatedly in deciduous broadleaf canopies across the network, are essentially dominated by changes in the respective pigment concentrations. Using the model we were able to explain why this spring maximum in green signal is often observed out of phase with the maximum period of canopy photosynthesis in ecosystems across Europe.


Introduction
Within Europe, continuous flux measurements of CO 2 , water and energy exchange between ecosystems and the atmosphere started in the early 1990s at a handful of forest sites (Janssens et al., 2001;Valentini et al., 2000). Nowadays, through the realisation of large European programmes such as EUROFLUX and CARBOEUROPE-IP amongst others, the number of natural and managed terrestrial ecosystems where the dynamics of water and CO 2 fluxes are monitored continuously has increased tremendously (Baldocchi, 2014;Baldocchi et al., 2001), and that number is set to be maintained in Europe for at least the next 20 years as part of the European Integrated Carbon Observation System (ICOS, www.icos-infrastructure.eu/). This long-standing coordinated European network, placed across several important biomes, has already documented large interannual variability in the amount of CO 2 sequestered over the growing season (Delpierre et al., 2009b;Le Maire et al., 2010;Wu et al., 2012) and witnessed both the short-lived and long-term impacts of disturbance (Kowalski et al., 2004), heat waves  and management practices Magnani et al., 2007;Soussana et al., 2007) on the carbon and water balance of terrestrial ecosystems. As a direct result of such an observational network, it is now possible to estimate with greater confidence how evapotranspiration (ET) and net ecosystem CO 2 exchange (NEE) have responded to changes in climate over recent years (Beer et al., 2010;Jung et al., 2010) and better constrain our predictions of how ecosystems are likely to respond in the future to changes in climate using land surface and biogeochemical cycle models (Friend et al., 2007;Krinner et al., 2005).
The growth of new leaves every year is clearly signalled in atmospheric CO 2 concentration records and exerts a strong control on both spatial and temporal patterns of carbon (C) sequestration and water cycling (Keeling et al., 1996;Piao et al., 2008). Hence, for the purpose of understanding patterns and processes controlling C and water budgets across a broad range of scales, there are obvious advantages in creating explicit links between flux monitoring, phenological observation and biogeochemical studies (Ahrends et al., 2009;Baldocchi et al., 2005;Kljun, 2006;Lawrence and Slingo, 2004;Richardson et al., 2007;Wingate et al., 2008). Leaf phenology has fascinated human observers for centuries and is related to external signals such as temperature or photoperiod (Aono and Kazui, 2008;Demarée and Rutishauser, 2009;Linkosalo et al., 2009). In the modern era, phenology has gained a new impetus, as people have realised that such records must be sustained over many years to reveal subtle changes in plant phenology in response to climate change (Rosenzweig et al., 2007) and improve our understanding of the abiotic but also biotic (metabolic and genetic) triggers that determine seasonal changes in plant development. Currently, descriptions of phenology in dynamic vegetation models are poor and need to be improved and tested against long-term field observations if we are to predict the impact of climate change on ecosystem function and CO 2 sequestration (Keenan et al., 2014b;Kucharik et al., 2006;Richardson et al., 2012).
Variations in the concentrations of pigments and spectral properties of leaves also provide a valuable mechanistic link to changes in plant development and photosynthetic rates when interpreted with models such as PROSAIL that combine our knowledge of radiative transfer through forest canopies and leaf mesophyll cells with leaf biochemistry (Jacquemoud and Baret, 1990;Jacquemoud et al., 2009). Automated techniques to detect and assimilate changes in the optical signals of leaves either near the canopy or remotely from space in conjunction with the coupling of radiative transfer models with biogeochemical models will help to improve the representation of leaf phenology and physiology in dynamic vegetation models (Garrity et al., 2011;Hilker et al., 2011;Klosterman et al., 2014).
This paper aims to synthesise data from flux sites across Europe where researchers have embraced the opportunity to establish automatic observations of phenological events by mounting digital cameras and recording daily (or even hourly) images of the vegetation throughout the seasons. We examine whether coherent seasonal changes in digital image properties can be observed at the majority of sites even if camera types and configuration are not yet harmonised and calibration procedures are not fully developed. In this work our two main objectives were to (1) assess how well digital images can be processed automatically to reveal the key phenological events such as leaf out, flowering, leaf fall, or land management practices such as mowing and harvesting using a piecewise regression approach; and (2) provide a mechanistic link between digital images, leaf phenology and the physiological performance of leaves in the canopy. To address the second objective, we adapted the model PROSAIL (Jacquemoud and Baret, 1990;Jacquemoud et al., 2009) to simulate the seasonal changes in red, green and blue signals detected by digital cameras above canopies and performed sensitivity analysis of these signals to variations in canopy structure (leaf area, structure and angles) and biochemistry (leaf pigment and water content).

Study sites and camera set-up
The European network of digital cameras currently covers over 50 flux sites across Europe ( Fig. 1 and Table 1). Table 1 demonstrates that a diverse selection of commercially available cameras are being used across the network. The cameras are installed in waterproof housing that is firmly attached to the flux tower some height above the top of the canopy. The cameras are generally orientated north, looking slightly down to the horizon to ensure that the majority of the image contains the vegetation of interest. Canopy images of the same scene are taken repeatedly each day, at a frequency that varies across the different sites from between 1 and 12 images per day and stored as 8-bit JPEG files (i.e. digital numbers rang-  Table 1. ing from 0 to 255). The archived images used in the present analysis are all taken between 11:00 and 13:00 local time (LT). The camera setup is specific to each camera type but a common requirement to observe the seasonal colour fraction time series is to set the colour balance to "fixed" mode (on the Stardot cameras) or the white balance to "manual" mode (for the Nikon Coolpix cameras) (Mizunuma et al., 2013(Mizunuma et al., , 2014.

ROI selection and colour analysis
For each site, a squared region of interest (ROI) was selected from a visual inspection of the images. The ROI had to be as large as possible and common to all images of the same growing season, while including as many plants as possible but no soil or sky areas. Image ROIs were then analysed using the open-source image analysis software, Image-J (Image-J v1.36b; NIH, MS, USA). A customised macro was used to extract red-green-blue (RGB) digital number (DN) values between 0 and 255 (n colour ) for each pixel of the ROI, and a mean value for all pixel values in a given ROI of each image was calculated. Various colour indices can be obtained from digital image properties (Mizunuma et al., 2011), but here we used only the chromatic coordinates, hereafter called colour fraction (Richardson et al., 2007): colour fraction = n colour /(n red + n green + n blue ), where colour is either red, green or blue. In the following we will thus refer to the "green fraction" as the mean green colour fraction for all pixels in the ROI, as www.biogeosciences.net/12/5995/2015/ Biogeosciences, 12, 5995-6015, 2015 opposed to the amount of pixels covered by vegetation in the entire image (Comar et al., 2012).

Data filtering procedure
Image quality is often adversely affected by rain, snow, low clouds, aerosols, fog and uneven patterns of illumination caused by the presence of scattered clouds. These influences often create noise in the trajectories of colour indices, and are an important source of uncertainty that can hamper the description of canopy seasonal variations and the derivation of robust phenological metrics from colour index time series. To remove problematic images that were affected by raindrops, snow or fog from the digital photograph analysis, we used a filtering algorithm based on the statistical properties of the time series; we used two steps to filter the raw data. The first filter, based on the deviation from a smoothed spline fit as described in Migliavacca et al. (2011), was used to remove outliers. Thereafter, we applied the method implemented in Sonnentag et al. (2012), to reduce the variability of the colour fractions (Fig. 2).

Piecewise breakpoint change analysis
We used a piecewise breakpoint regression approach to extract the main phenological events automatically, such as leaf emergence and senescence from the colour fraction time series (Fig. 2). The procedure is implemented in the R package strucchange (Zeileis et al., 2002(Zeileis et al., , 2003 and is used to detect breaks in a time series by identifying points where the multiple linear correlation coefficients shift from one stable regression relationship to another (Bai and Perron, 2003). The 95 % confidence intervals of the identified breakpoints were then computed using the distribution function proposed by Bai (1997). To obtain credible breakpoints in complex green fraction time series such as in highly managed sites (e.g. grazed or cut grasslands, multiple rotation crops) the possibility of identifying numerous peaks or breakpoints (up to eight breaks) may be necessary. However, such a high number of breakpoints would be excessive in natural ecosystems, and so we decided to set a maximum of five breakpoints per growing season, for both managed and natural ecosystems. We opted for a breakpoint approach over other commonly used methods to extract phenological transitions from time series (i.e thresholds, derivative methods) because it can be considered as more robust and less affected by noise in the time series (Henneken et al., 2013). This is particularly relevant for our application that encompasses a large data set consisting of many different camera set-ups (camera type, target distance, image processing).

Determining phenophases by visual assessment
In order to relate the breakpoint detection method to phenological phases we also visually examined images from broadleaf forest ecosystems for leafing out, senescence and leaf fall. Six pre-trained observers looked through the same daily images and used a common protocol to identify dates when (1) the majority of vegetation started leafing out (i.e. when 50 % of the ROI contains green leaves), (2) the canopy first started to change colour, first to non-green colours such as yellow and orange, in autumn and (3) the last day when a few non-green leaves were still visible on the canopy before the day the branches became bare. These visually assessed dates were then averaged across observers and compared to the relevant breakpoints, identifying the same phenological stage. The leafing out phase was associated to the first automatically detected breakpoint, leaf senescence to the penultimate breakpoint and leaf fall to the last breakpoint. Based on this classification the key dates identified by the algorithm and visual inspection were consistently correlated with one another (Fig. 3). However, there was a tendency for the automatic algorithm to identify all the phenological transitions before the visually assessed dates, by about a week. Also, visual inspections had larger standard deviations, especially during canopy senescence. Because of this systematic difference between the two methods, the breakpoints indicating the start of the growing season were in agreement with visually inspected dates to within 9 days only, and an even lower accuracy was found for leaf senescence and leaf fall (RMSE of 15.9 and 13.3 days, respectively). These RMSE values fall, however, within the range recently found by Klosterman et al. (2014)

Radiative transfer modelling
To interpret the colour fraction time series produced by cameras in the network mechanistically, we combined the bidirectional radiative transfer model PROSAIL (Jacquemoud and Baret, 1990;Jacquemoud et al., 2009) with some basic spectral properties of the photodetectors and the transmittances of the optical elements and filters used in the cameras, all combined into the so-called spectral efficiency of the RGB colour channels (G RGB , in DN W 1 , defined as the digital number per watt). This version of PROSAIL has been made available to the research community on the following repository (https://bitbucket.org/jerome_ogee/webcam_ network_paper). This version of PROSAIL is required to make the necessary link between the DN values measured by the camera sensors to the reflectances simulated by PRO-SAIL. It is also necessary to take into account the camera set-up and specific characteristics. Examples of these properties are shown in Supplement (Figs. S1 and S2) for two types of camera from the network: (1) the Stardot NetCam CS5, in use at the majority of the European sites (Table 1) and the dominant camera used in the North American PhenoCam Network (http://phenocam.sr.unh.edu/webcam/) (Toomey et al., 2015) and (2) the Nikon Coolpix 4500, in use at two sites within the European camera network, and ca. 17 sites within the Japanese Phenological Eyes Network (PEN) (http: //pen.agbi.tsukuba.ac.jp/index_e.html) (Nasahara and Nagai, 2015). The PROSAIL model combines the leaf biochemical model PROSPECT (PROSPECT-5) that simulates the directional-hemispherical reflectance and transmittance of leaves over the solar spectrum from 400 to 2500 nm (Jacquemoud and Baret, 1990) and the radiative transfer model SAIL (4SAIL; Verhoef, 1984). The SAIL model assumes that diffusers are randomly distributed in space (turbid medium assumption) and thus may not be applicable for very clumped canopies such as sparse forests or crops (e.g. vineyards orchards . . . ). Furthermore, the radiative transfer model assumes only one type of foliage, and therefore cannot deal well with species mixtures. However, for mixed forests, PROSAIL can still be used to interpret RGB signals if the ROI selected on the images is dominated by one single species. The version of the model used here (version 5B for IDL http://teledetection.ipgp.jussieu.fr/prosail/) requires 11 parameters from PROSAIL (leaf area, leaf angle, leaf mass and chlorophyll, carotenoid, water or brown pigment contents, hotspot parameter, leaf structural parameter and dry soil fraction, percentage of diffuse light) as well as geometrical parameters (sun height, view zenith angle and sun-view azimuthal difference angle). The percentage of diffuse light (ϕ diffuse ) is necessary to calculate the amount of direct and diffuse incoming radiation spectra at the top of the canopy and is estimated here from observed incoming global radiation (R g , in W m −2 ) using a procedure developed by Reindl et al. (1990). The incoming radiation spectra at the top of the canopy are then estimated from R g , ϕ diffuse and mean, normalised spectra for direct (I direct ) and diffuse (I diffuse ) radiation derived by François et al. (2002) using the 6S atmospheric radiative transfer model (Fig. S3). The 6S simulations performed by François et al. (2002) considered a variety of aerosol optical thicknesses at 550 nm (corresponding to a visibility ranging from 8.5 to 47.7 km), water vapor content (from 0.5 to 3.5 g cm −2 , corresponding to most situations encountered in mid-latitudes), solar incident angle (from 0 to 78.5 • ) and standard values for ozone, CO 2 and other atmospheric constituents but did not consider the presence of clouds. Therefore the derived spectra are only valid for cloud-free conditions and, following François et al. (2002), values of ϕ diffuse below 0.5 were considered to represent such cloud-free sky conditions. PROSAIL was then used to estimate the amount of light reflected by the canopy in the direction of the camera for each wavelength E(λ) (W m −2 nm −1 ) from which we compute the RGB signals according to where λ UV and λ IR are the UV and IR cut-off wavelengths of the camera sensor and filter (see Fig. S1) and B RGB is a constant factor that accounts for camera settings (mostly colour balance) and was manually adjusted for each RGB signal and each camera/site using a few days of measurements outside the growing season. The modelled RGB colour fractions were then computed in a similar fashion as in Eq.
(1) but with I RGB instead of n colour . In practice, G RBG is often normalised to its maximum value rather than expressed in absolute units and for this reason, I RGB is not a true digital number, but this has no consequence once expressed in colour fractions. Also, we are aware that the image processing of real (nonuniform) scenes is far more complex than Eq.
(2) (Farrell et al., 2012) but from the preliminary results presented below, this simplistic formulation is robust enough to describe time series of average colour fractions over a large and fixed ROI, measured with different camera settings.

CO 2 flux analysis
A recent study by Toomey et al. (2015) demonstrated that the seasonal variability in daily green fraction measured at 18 different flux sites generally correlated well with GPP estimated from co-located flux towers over the season. However, they found that the correlation between daily GPP and daily green fraction was better for some plant functional groups (grasslands) than for others (evergreen needleleaf forests or deciduous broadleaf forests). Thus for some of the flux sites presented in Table 1, we also qualitatively related changes in vegetation indices to gross primary productivity (GPP) time series. Net CO 2 fluxes were continuously measured at each site using the eddy-covariance technique described in Aubinet et al. (2012). Level 4 data sets were both qualitychecked and gap-filled using online eddy-covariance gapfilling and flux-partitioning tools provided by the European flux database cluster Reichstein et al., 2005). Full descriptions of the flux tower set-ups used to compare digital images and GPP are provided in Wilkinson et al. (2012) for Alice Holt (deciduous forest), Wohlfahrt et al. (2008) and Galvagno et al. (2013) for Neustift and Torgnon (sub-alpine grasslands), Vesala et al. (2010) for Hyytiälä (evergreen forest) and Aubinet et al. (2009) for Lonzee (cropland).

Needleleaf forest ecosystems
Differences in the seasonal evolution of canopy green fractions are presented for individual years at a selection of needleleaf forest flux sites spanning a latitudinal gradient of approximately ca. 30 • (31 • 20 -61 • 50 N) ( Fig. 4 and Table 1). It was often difficult to automatically determine the start of the growing season for the coniferous sites, either because of snow cover (noticeable changes in the colour signals were often associated with the beginning and end of snow cover) or because of problems caused by the set-up of the camera (either too far away from the crown or the view contains too much sky). However seasonal changes in the amplitude of the canopy green fraction of needleleaf trees were generally conservative across sites and often displayed a gentle rise in green fraction values during the spring months (Fig. 4a) and a gentle decrease during the winter months. In contrast, the deciduous Larix site in Italy showed a steep and pronounced start and end to the growing season that lasted approximately 5 months (Fig. 4b).
The main evergreen conifer species Pinus sylvestris L. and Picea abies L. exhibited different green fraction variations during spring. In Picea species, new shoots contain very bright, light green needles that caused a noticeable increase in the green colour fraction during the months of May and June (Fig. 4, e.g. Tharandt and Wetzstein). In contrast, the new shoots of Pinus species primarily appeared light brown as the stem of the new shoot elongated (Fig. S4) and caused a small reduction in the green fraction, only detectable at a number of sites and years (Brasschaat in 2009(Brasschaat in , 2012Norunda in 2011;Hyytiälä in 2008Hyytiälä in , 2009Hyytiälä in , 2010Hyytiälä in , 2012, followed by an increase in the green fraction values once needle growth dominates at the shoot scale ( Fig. 4 and S4).
The breakpoint analysis could identify a change in the green fraction when the canopy became consistently snowfree or covered in snow, but also when it experienced sustained daily mean air temperatures of above 0 • C (Fig. 5). For instance, at the Finnish site Hyytiälä, the green fraction exhibited rapid changes during final snowmelt (day 18, bp1) and first snowfall (day 346, bp5 confirmed with visual inspection of the images), but also around days 100-110 (bp2) as daily mean temperatures increased from −10 • C and stabilised for several days just above 0 • C and GPP rates increased (Fig. 5). Towards the end of the growing season significant changes in the green fraction were again detected by the piecewise regression approach, one around day 267 (bp3), coinciding with the period when minimum air temperatures began to fall below 0 • C and another around day 323 (bp4) coinciding with a short period of warmer temperatures (Fig. 5). Interestingly from the first breakpoints at the beginning of the year, day 115 (bp2)  323 (bp4) coincided with a stable period of no photosynthesis (GPP = 0 gC m −2 d −1 ) and a transient increase of daily mean temperature above 0 • C. However, the start of this stable period with no photosynthetic uptake was not detected by a breakpoint. During this period the green fraction tracked fluctuations in the daily mean temperature. The recovery of biochemical reactions and the reorganisation of the photosynthetic apparatus in green needles is known to be triggered when air temperature rises above 0 • C (Ensminger et al., 2004(Ensminger et al., , 2008) and at about 3-4 • C at this particular site in Finland (Porcar-Castell, 2011;Tanja et al., 2003). Our results suggest that green fractions from digital images seem sensitive enough to detect these changes in the organisation of the photosynthetic apparatus of the coniferous evergreen needles as they acclimate to cold temperatures.

Grassland and cropland ecosystems
Many land surface models still lack crop-related plant functional types and often substitute cropland areas with the characteristics of grasslands (Osborne et al., 2007;Sus et al., 2010). Our initial results from the European phenology camera network demonstrate the difficulty of teasing apart the seasonal and interannual developmental patterns as they are often complicated by co-occurring agricultural practices (e.g. cutting, ploughing, harvesting and changes in animal stocking density) (Fig. 6).
Overall we found that there was surprisingly little difference between the onset dates of growth for most of the permanent grassland sites, despite being located in very different locations across Europe; however, the onset of the green fraction signal was considerably delayed at the sub-alpine grassland site in Torgnon compared with the other sites (Fig. 6). Torgnon also had the most compressed growing season of all the sites, encompassing a period of less than 100 days, some 100 days shorter than the other sites. These differences in growing season length between permanent grassland sites are caused by elevation (most grassland sites are situated at 1000 ± 50 m a.s.l., while the Italian site, Torgnon, is located at ca. 2160 m a.s.l.) that induced differences in temperature and snow cover and consequently the RGB signals.
Management practices such as the cutting of meadows (e.g. Neustift and Früebüel) and changes in animal stocking  Fig. S5). In contrast, the blue signal tended to increase in strength during flowering periods, making the impact of the mowing events dramatic when they occurred (Fig. 7). If the piecewise regression algorithm was set to allow the detection of up to eight breakpoints at those sites, flowering events were often identified as well as the mowing events but this made the detection of the start and end dates of the growing season even more challenging using automated algorithms, as manual inspection of each breakpoint would be necessary to identify the dates of first leaf growth and leaf death. For example when eight breakpoints were observed in grassland ecosystems the first breakpoint was frequently caused by the start of the snow-free period, as opposed to the start of growth. Subsequent breakpoints typically indicated the phenology and management of the vegetation, particularly mowing and flowering, and suggest that a visual inspection of images may still be necessary to clarify the nature and management causes behind breakpoints in some grassland sites. Thus having the option to detect up to eight breakpoints or more could be an advantage for ecological studies on the phenology of flowering, pollination and community dynamic responses to environmental changes.
In addition, at some sites it was not always easy to discern visually from images when a grassland started its first growth, as often fresh shoots are hidden by litter and dead material from the previous year. This particular problem may lead to a slight overestimation of the start date of growth and additionally lead to a potential temporal mismatch between GPP and green-up signals. For example at the subalpine site in Torgnon the green fraction and GPP peaked at the same time of year (bp2), but the onset of the green signal in spring lagged behind the onset of photosynthesis by about 7-10 days (Fig. 8). Given these results at Neustift, Torgnon and the other sites within the network, grassland green signals can provide additional information on the variations in GPP over the season and between years as suggested by Migliavacca et al. (2011) and more recently Toomey et al. (2015), but an underestimation of the growing season length may also occur using our automatic breakpoint detection approach.
As mentioned above, the use of grassland characteristics to describe the phenology of different crops is common in some land surface models. However, at the same site we can see that the start of the growth period for each different crop varies widely over commonly applied rotations. For example, at Lonzee in Belgium we found that for individual crop types such as potato and winter wheat the patterns of the green signal were fairly similar between different years (Fig. 9). However, we could still observe variations between years such as a slightly shorter green period for the potatoes in 2014 rela- tive to 2010 or a slightly later start to the season for winter wheat in 2013 relative to 2011.
As huge areas of Europe are dedicated to the production of crops (326 Mha) and grasslands (151 Mha) (Janssens et al., 2003), and are known to be one of the largest European sources of biospheric CO 2 to the atmosphere at a rate of about 33 TgC yr −1 (Schulze et al., 2009), it is increasingly important that more field observations of developmental or phenological transitions are obtained to constrain the timing and developmental rate of plants in situ and improve model simulations. Our initial results from the European phenology camera network show that digital repeat photography may provide a valuable assimilation data set in the future for providing useful indicators of developmental transitions and agricultural practices (e.g. cutting, ploughing, harvesting and changes in animal stocking density).

Broadleaf forest ecosystems
Within the European phenology camera network, the majority of broadleaf forest sites are deciduous and are identified easily by the seasonal RGB signals they produce. Typically the start of the temperate growing season coincides with a strong increase in the relative green fraction (Fig. 10). Across the network the timing of this spring green-up varies Figure 8. Time series for the alpine grassland, Torgnon in Italy demonstrating (a) the daily mean air temperature (b) the gross primary productivity and (c) the calculated breakpoints for the green fraction using the piecewise regression approach. Vertical solid lines show major breakpoint changes, identifying important transitions in the green fraction over the growing season slightly with latitude, usually occurring first in the southern sites and moving north, with the British sites starting later than continental sites at similar latitude (for the same years data not shown) (Fig. 10). The end of the growing season, identified clearly as a decrease in the green signal, varied considerably across the network. The more continental sites such as the Hainich and Lägeren deciduous forests exhibited the shortest growing season lengths, whilst the oceanic and Mediterranean sites had far longer growing seasons. A high degree of variability in the timing for colour changes is commonly found in autumn across temperate deciduous ecosystems (Archetti et al., 2008) despite the fact that, at least in the case of European tree species, changes in photoperiod and air temperature are usually considered as the main drivers of the colouration of senescent leaves (Delpierre et al., 2009a;Keskitalo et al., 2005;Menzel et al., 2006).
Interestingly, the evergreen broadleaf forests at the Mediterranean sites in Spain and France displayed similar RGB seasonal variations (Fig. 10b). However, the peak in the green fraction values were observed somewhat later in spring compared to those of the deciduous broadleaf sites. These maximum green fraction values are most probably linked to the production of new leaves that typically occurs at this pe- 2003). Interestingly, a strong decrease in the green fraction was observed prior to the peak at the Spanish site, Majadas del Tieter. Similar patterns in the NDVI time series have been observed around the same period, at the Puechabon site but for different years to the one studied here (2006-2008) (Soudani et al., 2012). This drop in NDVI was explained by the shedding of old leaves coinciding with the period of leaf sprouting in spring. However, on inspection of the Spanish site photos (Fig. S6) we found that the canopy during this period was covered in conspicuous male catkin-type flowers that appear yellow-brown in the images (Fig. 10). Thus, in the case of the Spanish site at least, the strong decrease in the green fraction seemed dominated by a male flowering event, in addition to the shedding of old leaves. In the case of the Puechabon site, visual inspection of the photos did not detect a strong flowering event, however, the camera is located slightly further away from the canopy, making it difficult to detect flowers easily by eye. However, phenological records maintained at the site indicate that the period between bp2 and bp3 when the green signal slightly decreases, coincides with the start and end of the male flowering period, as well as leaf fall. In contrast, the signal between bp3 and bp4 indicates the period of leaf flushing for this year. Further studies comparing NDVI signals with digital images should allow us to understand the observed variations in both signals better and their link to phenological events such as flowering and litterfall in evergreen broadleaves and how these vary between years in response to climate. For broadleaf deciduous (and to some extent evergreen) species, our breakpoint approach also detected a significant decline in the green fraction a few weeks after leaf emergence and well before leaf senescence (Figs. 2 and 10). This pattern in the green fraction has also been observed for a range of deciduous tree species in Asia and the USA ( Keenan et al., 2014a;Nagai et al., 2011;Saitoh et al., 2012;Sonnentag et al., 2012;Toomey et al., 2015;Yang et al., 2014). In addition, this feature has also been observed at the leaf scale using scanned images of leaves (Keenan et al., 2014a;Yang et al., 2014) and at the regional scale for a number of deciduous forest sites using MODIS surface reflectance products Keenan et al., 2014a). The causes that underlie the shape and duration of this large peak at the beginning of the growing season remain unclear at present. Recent camera studies that have also measured either leaf or plant area index at the same time have found no dramatic reductions in leaf area during this rapid decline in the green fraction following budburst (Keenan et al., 2014a;Nagai et al., 2011). At two of the deciduous sites within our network, Alice Holt and Søro (Fig. 10), daily photosynthetically active radiation (PAR) transmittance was also measured, providing a suitable proxy for changes in canopy leaf area. In both cases no decrease in the leaf area index (LAI) proxy was detected during the decrease of the green signal shortly after budburst. If this decrease in the green signal after leaf growth is not caused by a reduction in the amount of foliage in most cases, it is likely associated with either changes in the concentration and phasing of the different leaf pigments or changes in the leaf angle distribution. These different hypotheses are tested in the next section.

Sensitivity analysis of model parameters
Using the PROSAIL model as described above, with the camera sensor specifications of the Alice Holt oak site (see Fig. S1) we performed three different sensitivity analyses of the simulated RGB fractions to the 13 model parameters (Table 2). All sensitivity analyses consisted of a Monte Carlo simulation of between 2000 and 10000 runs each. For the first analysis the model was allowed to freely explore different combinations of the parameter space over the range of values commonly found in the literature and with no constraints on how the parameters were related to each other (all parameters being randomly and uniformly distributed). The results from this initial sensitivity analysis indicated that the RGB signals were sensitive to four parameters: the leaf chlorophyll ([Chl]), carotenoid ([Car]) and brown contents (C brown ) and the leaf structural parameter (N) (see Supplement, Fig. S7). In contrast, the simulated RGB signals were relatively insensitive to leaf mass (LMA), leaf water content (EWT) and, to some extent, to LAI (above a value of ca. 1). This sensitivity analysis also nicely demonstrates how measurements of NDVI made above canopies are most strongly influenced by LAI and to a slightly lesser extent by leaf pigment contents. Recent studies have also discussed whether changes in leaf inclination angle have a strong affect on the green fraction over the growing season (Keenan et al., 2014a). Our simulations for two different camera view angles (one looking across the canopy and another looking down (results not shown)) indicate that this particular parameter did not have a strong impact on the RGB fractions (Fig. S7).
The model also demonstrated that the impact of diffuse light had a negligible impact on the green fraction. In contrast, an increase in the blue fraction alongside a decrease in the red fraction was observed as the percentage of diffuse light increased (Fig. S7). This result was anticipated given the different spectra prescribed for incoming direct and diffuse light (Fig. S3). However, it is also important to bear in mind that colour signals throughout this manuscript are expressed as relative colour fractions (Eq. 1), and thus strong increases in one or more colour fraction, influenced by strong changes in a particular parameter (e.g. [Chl]. . . ), must also be compensated for by decreases in one or more of the other colour signals.
In a second sensitivity analysis conducted over the same time frame, we refined our assumptions on how certain parameters were likely to vary with one another in spring during the green-up. For this, we fixed all parameters to values typical for English oak during spring conditions (Demarez et al., 1999;Kull et al., 1999), except for LAI and the concentrations of chlorophyll and carotenoid. We then imposed two further conditions stating that (1) leaf chlorophyll contents increased in proportion to LAI (i.e. the ratio of [Chl] / LAI was normally distributed) and (2) carotenoid and chlorophyll contents also increased proportionally (i.e. [Car] / [Chl] was normally distributed around 30 ± 15 %). This ratio between pigment contents is commonly found in temperate tree species. For example, Feret et al. (2008, see Fig. 3 therein) showed a strong correlation between the concentrations of chlorophyll a/b and carotenoids for a range of plant types. This second sensitivity analysis revealed clearly how the RGB fractions would likely respond to LAI, chlorophyll and carotenoid contents during the spring green-up (Fig. 11). Firstly, we observed that the RGB signals were more sensitive to LAI, [Chl] and [Car] when the analysis was constrained by our a priori conditions (compare panels from Fig. 11 with those of Fig. S7). In addition we found that the RGB fractions varied in sensitivity across the range of LAI variations typically found in deciduous forests. Most of the sensitivity in the green signal was found at very low values of LAI (< 2), whereafter the signal became insensitive. In contrast, the NDVI signal remained sensitive throughout the full range of prescribed LAI values. For the range of likely [Chl] taken from the literature for oak species (Demarez et al., 1999;Gond et al., 1999;Percival et al., 2008;Sanger, 1971;Yang et al., 2014), our simulations indicated that the sharp increase in the green signals observed by the camera sensors during leaf out are mostly caused by an increase in [Chl]. More interestingly, and contrary to our previous analysis where changes in [Chl] and [Car] were not correlated (compare the same panels in Fig. 11 with those presented in Fig. S7), this new analysis clearly shows that when [Chl] Jacquemoud et al. (2009) reached ca. 30 µg cm −2 , the green signal begins to respond negatively to a further increase in [Chl]. This is because in this simulation, an increase in [Chl] is accompanied by an increase in carotenoids and the green fraction responds negatively to an increase in [Car] (Figs. 11 and S7). Another interesting feature of this sensitivity analysis is how the dependence of NDVI on pigment content was greater when we imposed the two constraints described above. This investigation with the PROSAIL model therefore suggests that the sharp reduction in the green fraction commonly observed after budburst in deciduous forests is mostly driven by an increase in leaf pigment concentrations.

Modelling seasonal RGB patterns mechanistically
Using the PROSAIL model and a few assumptions on how model parameters varied over the season (Table 2, Fig. 12) we investigated the ability of the PROSAIL model to simulate seasonal variations of the RGB signals measured by two different cameras installed at the Alice Holt site. For this, a proxy for seasonal variations in LAI was obtained by fitting a relationship to values of LAI estimated from measurements of transmittance as described in Mizunuma et al. (2013), while the expected range and seasonal variations of [Chl], [Car], C brown and N were estimated from several published studies on oak leaves (Demarez et al., 1999;Gond et al., 1999;Percival et al., 2008;Sanger, 1971;Yang et al., 2014), as summarised in Table 2. We also manually adjusted the parameter B RGB in Eq. (4) to match the RGB values measured by the camera during the winter period prior to budburst (highlighted by the grey bar in Fig. 12). Using this simple model parameterisation, the model was able to capture the seasonal pattern and absolute magnitude of the three colour signals relatively well, including the spring "green hump" (Fig. 12). As anticipated from our sensitivity analysis (Fig. 11) we found that at the beginning of the growing season the initial rise in the green signal is dominated by small changes in LAI and [Chl], whilst the small decrease in the green signal that follows is most likely caused by the further increase of [Chl] and [Car]. This period is also characterised by a strong decrease in the red signal and a strong increase in the blue fraction that distinguishes it clearly from the effect of changes in LAI (as blue and green signals remain constant at LAI > 2). This relationship is expected as the absorption by chlorophyll and carotenoids are not the same in all wavelengths (Feret et al., 2008, see Fig. 8 therein) and can account for most of the observed colour signal trends. Interestingly the maximum green signal observed by the camera coincides with the time when LAI reached half of its maximum value (dotted line on graphs), rather than when [Chl] is maximum. This theoretical result is consistent with pigment concentration data from other recent studies (Keenan et al., 2014a;Yang et al., 2014). Our analysis also indicates that the time when maximum [Chl] and/or [Car] is attained coincides more with the end of the green hump (i.e. several weeks after the peak in the green fraction) and also corresponds to when ecosystem GPP reaches its maximum. Recent pigment concentration measurements performed in the same forest canopy in 2012 and 2013 seem to confirm this theoretical result (data not shown).
Our observations across the deciduous forest sites of the network also highlight a second strong hump in the red and blue signals towards the end of the growing season linked to senescence of the canopy. This autumnal change in colour fractions starts with a decrease in both the green and blue signals and a strong increase in the red signal. To understand this combined signal we conducted a third sensitivity analysis with the model. In this Monte Carlo simulation we maintained all parameters of the model constant except those likely to be important during autumn, mainly LAI,  Figure 11. Sensitivity of modelled RGB fractions and NDVI for the NetCam camera at the Alice Holt deciduous broadleaf forest, as predicted by the PROSAIL model and constrained by Chl : Car and Chl : LAI ratios (see text and compare to those in Fig. S7). All other parameters are set to standard budburst values and the solar elevation is fixed at 60 • . The NDVI is computed using the camera view angle and the same wavebands as for MODIS NDVI (545-565 nm for red and 841-871 nm for near infrared). The bottom panel represents the frequency of estimates for a given value of a particular parameter. One can see that compared with Fig. S7 there are more likely to be values skewed to the lower pigment concentrations with our constraints compared with the no-constraints analysis of Fig. S7.
[Chl], [Car], C brown and N ( other hand were kept as for the sensitivity analysis shown in Fig. 11. The results from this third sensitivity analysis suggest that the start of the autumnal maximum is most likely caused by a reduction in [Chl] and [Car] of the foliage and a simultaneous increase in the structural parameter N and brown Biogeosciences, 12, 5995-6015, 2015 www.biogeosciences.net/12/5995/2015/ pigment concentration (Fig. 12). The degradation of green chlorophyll pigments in autumn has been observed in oak and other temperate forest species before. Interestingly, the rate of degradation in [Chl] and [Car] during the autumn breakdown is not always proportional and can often lead to a transient increase in the [Car] / [Chl] that allows the yellow carotenoid pigments to remain active as nutrient remobilisation takes place. Leaves of temperate deciduous broadleaf forests in Europe have an overall tendency to yellow in autumn and eventually turn brown before dropping (Archetti, 2009;Lev-Yadun and Holopainen, 2009). This subsequent browning of the foliage is the likely cause for the sharp decrease in the red signal and increase in the blue signal around day 325 (Fig. 12). Visual inspection of the digital images confirmed this was the case and was reflected in the modelling with the increase in C brown , the only parameter that can reconcile the strong reduction in the red signal and the increase in blue strength observed (Fig. 13). Thereafter, reductions in LAI may also contribute further to this trend as leaves fall from the canopy. Results shown in Figs. 11-13 have been obtained using spectral properties of the Stardot camera, that we obtained directly from the manufacturer (D. Lawton, personal communication, 2012). However, at Alice Holt, both a Stardot camera and a Nikon camera have been operating since 2009 (Mizunuma et al., 2013). Using a similar approach but with the Nikon camera, we then tested the idea that colour fractions seen by this other camera would also suggest similar variations in leaf pigment and structural parameters over the season. Spectral characteristics of the Nikon Coolpix 4500 were characterised at Hokkaido University in Japan and are shown in the Supplement (Fig. S2). We then used these camera specifications and the same parameterisation of canopy structural properties (LAI and N ) and leaf pigment concentrations as the those used in Fig. 12. We also manually adjusted B RGB using the same procedure as for the Stardot camera.
From this camera comparison we found that, in order to match the independent RGB camera signals with the same radiative transfer model, only the value of C brown during the growing season had to be adjusted (Fig. S8). On inspection of the camera images this may be justified as the scenes and ROI captured by each camera are very different despite looking at the same canopy, as the Stardot looks across the canopy, whilst the Nikon has a hemispherical view looking down on the canopy. Based on the need for different levels of C brown , it is suggested that more brown (woody material) occupies the ROI of the Stardot camera during the vegetation period compared with that in the Nikon image, as confirmed in the images. Besides this difference in C brown , the same model parameterisation reasonably captured the seasonal features of all three colour signals.  Figure 13. Results from the PROSAIL sensitivity analysis as performed in Fig. 12 except that the leaf brown pigment (C brown ) and structural paramer (N) are now allowed to vary independently of the other parameters, as expected during senescence.

Technical considerations for the camera network
As demonstrated in this paper and elsewhere (Keenan et al., 2014a;Migliavacca et al., 2011), researchers are now exploring the link between canopy colour signals and plant physiology in order to maximise the utility of this relatively inexpensive instrument. However, for this step to proceed further, the network should address several technical considerations. Firstly, the digital cameras in our European network are currently uncalibrated instruments unlike other commonly used radiometric instruments. In addition, these cameras are deployed in the field and are often exposed to harsh environmental conditions. Thus their characteristics may drift over time. For example although the CMOS sensors (commonly found in the cameras of our network) do not age quickly over periods of several years, the colour separation filters on the sensor plate may age after time because of UV exposure. Most commercial cameras already contain UV filter protection and in addition they are protectively housed from the elements and sit behind a glass shield that protects the camera www.biogeosciences.net/12/5995/2015/ Biogeosciences, 12, 5995-6015, 2015 from the damage by UV or other environmental conditions, such as water. Also, the studies of Sonnentag et al. (2012) and Mizunuma et al. (2013Mizunuma et al. ( , 2014 demonstrated that different camera brands and even different cameras of the same brand that produce different colour fraction time series reflecting differences in their spectral response still produce coherent phenological metrics. Nonetheless, it would be useful to develop a calibration scheme by digital photography of radiometrically characterised colour sheets such as those used to detect the health and nutritional status of plants (Mizunuma et al., 2014) whilst the camera is in the field and exposed to a variety of light conditions. However, the deployment of such a colour checker for long-term continuous monitoring is problematic as the spectral quality of the colour checker will alter over time as particles, such as dust and insects, accumulate on its surface. An alternative solution could involve routinely checking for sensor drift using images taken during the winter months of different years. For example, at least in the case of deciduous broadleaf species, 1 or 2 months before the beginning of the growing season (see grey bar in Fig. 12), the RGB signals tend to display relatively constant values. Provided the camera does not move and the colour balance settings are not changed, this period could be used to detect problems with the camera when compared over several years. At Alice Holt where two cameras have been operating for at least 5 years we have found no evidence for such drift using winter time values (data not shown). Ideally, it would be desirable to use only digital cameras from manufacturers that provide information on the spectral characteristics of the sensors, filters and algorithms used, or to measure routinely the spectral response of individual cameras within the network as the cameras age. Secondly, digital camera studies have typically focussed on using the green fraction to deduce the dates of canopy green-up and senescence. The green fraction has been preferred because its signal-to-noise ratio is higher in vegetated ecosystems. For this very reason, and because of inter-pixel dependencies caused by the cameras built-in image processing, we recommend setting the camera up so that the images do not contain too much of the sky area, ideally less than 20 %. In addition, as demonstrated in this study, combining the information of all three colour signals may provide more useful information on canopy physiology, phenology and management impacts. However, to look at these particular signals we must understand the effects of other environmental factors that can impact the day-to-day variability of the RGB colour fractions. In particular, colour signals are sensitive to the spectral properties of the incoming light and thus to the percentage of diffuse radiation. Using the PRO-SAIL model we explored the impact of diffuse radiation on the RGB signals at the Alice Holt site and found that the red and blue fractions were much more affected by rapid changes in sky conditions than the green fraction. When simulating the impact of diffuse light on the RGB response we make use of direct and diffuse light spectra presented in Fig. S3. From the spectra shown in Fig. S3 it is clear that increases in diffuse light will result in more energy received in the blue wavelengths and less energy in the red, thereby influencing the red and blue fractions directly. This is why there is very little change in the green fraction over the range of diffuse light (%). Furthermore, by incorporating the day-to-day variations in diffuse radiation at the site the model did not reproduce better the red and blue fractions, even for cloud-free conditions (Fig. 14), demonstrating that the influence of diffuse light was not easy to account for in the model. This can be problematic as the percentage of diffuse light, unlike other variables such as leaf area index or pigment concentrations, can change dramatically from one day to the next. For this reason, the green fraction is probably the best-suited signal for detecting rapid changes in leaf area and pigment concentrations, at least within their lower range of values (0-2 for LAI and 0-30 µg cm −2 for [Chl]); however, these results will need to be checked thoroughly for the different cameras across the entire network. Lastly, cameras within the network are so far not thermally regulated. Our preliminary results using the two most commonly used cameras in the networks (Stardot and Nikon Coolpix) seem to indicate that images are not sensitive to contrasting temperatures. This point will certainly need to be addressed more thoroughly in future studies and could be addressed by taking black images periodically to assess the level of instrument noise and its relationship with temperature. If this could be achieved then it would be possible to at least reduce noise in the signals caused by temperature. , 12, 5995-6015, 2015 www.biogeosciences.net/12/5995/2015/

Conclusions
This synthesis analysis of the European camera network demonstrates that using digital repeat photography at a daily resolution can aid the automatic identification of interannual variations in vegetation status, such as the emergence of new foliage (i.e. budburst, regrowth), flowering, fruit development, leaf senescence and leaf abscission, driven by climate. Furthermore, agricultural practices are captured well by the camera, providing a useful archive of images and colour signal changes that can be interrogated with complementary flux data sets. In the long term such data sets collected across the different networks of flux sites will become invaluable for investigating in detail the connections between climate and growing season length, and will contribute to a better understanding of the underlying controls on plant development and how these vary between plant functional types, species and location. Currently, large-scale estimates of photosynthesis are principally derived from models with highly uncertain parameterisations. In particular, large-scale vegetation models have difficulties predicting when the growing season starts and ends , how this relates to the carbon uptake period and the magnitude of photosynthetic rates attained during peak summertime. Traditionally, the remote sensing NDVI products are generally used for linking largescale variations of vegetation phenology and terrestrial primary productivity to climate change (Running et al., 2004;Zhou et al., 2003). However the spatial resolution of such measurements is typically low (> 100 m) and the temporal resolution often inhibitive for tracking the dynamic transitions of canopy cover over a mosaic of heterogeneous land surfaces, leading to problems in reconstructing the seasonal changes in canopy cover robustly (Hagolle et al., 2010, Quaife andLewis, 2010).
Our results indicate that by combining variations in the red, green and blue (RGB) colour fractions and mechanistic radiative transfer models, these digital archives can be used to quantify changes in the plants' physiological status and constrain estimates of NDVI (see Fig. S7). Near-surface measurements could thus provide RGB ground-truth estimates for similar remotely sensed multi-spectral satellite products  and provide the link between remotely sensed products and ecosystem physiology measured at flux towers (Migliavacca et al., 2011;Richardson et al., 2007;Wingate et al., 2008;Mizunuma et al., 2013). This technological breakthrough will provide a means of increasing our understanding of how canopy pigment contents vary between ecosystems and climates, and improving predictions of the CO 2 sequestration period and potential of terrestrial ecosystems at large scales.
The Supplement related to this article is available online at doi:10.5194/bg-12-5995-2015-supplement.