Evaluating MODIS vegetation products using digital images for quantifying local peatland CO2 gas fluxes

In peatlands plant growth and senescence affect annual ecosystem carbon dioxide (CO2) exchange, and CO2 fluxes show considerable inter‐annual variability. Phenology is a fundamental indicator of ecosystem carbon dynamics and can be measured from remote sensing systems, but the extent to which satellite products provide useful proxies of peatland vegetation phenology is not well known. Using MODIS vegetation products coupled with field observations of phenology from a basic camera system and measurements of Gross Primary Productivity (GPP) measured using a closed chamber method, we sought to establish the extent to which satellite observations capture phenological processes at a UK peatland site. Daily, true‐colour digital images were captured with a time‐lapse camera (Brinno TLC100) between 23‐Apr‐2013 and 29‐Oct‐2013 and converted into a Green‐Red Vegetation Index (GRVI). These were compared with a range of MODIS vegetation products at various spatial resolutions. We found that vegetation products with finer spatial resolution (<500 m) more accurately captured spring green‐up (e.g. Normalized Difference Vegetation Index 16‐day product), whereas those with 8‐day temporal resolution better captured whole‐season dynamics. The 8‐day Gross Primary Productivity (GPP8) and the fraction of absorbed photosynthetically active radiation (fPAR8) products had the strongest daily Pearson's correlations with camera‐derived GRVI (r > 0.90). The camera‐GRVI (P = 0.005, r = 0.98) and MODIS‐GRVI (P = 0.041, r = 0.89) products showed the strongest significant correlations with gross primary productivity measured using static chambers in the field. This work demonstrates that freely available MODIS data captured up to 92% of the daily variation in phenology over an upland peatland. This approach shows great potential for capturing ecosystem carbon dynamics which underpin carbon trading schemes, a budding funding source for peatland restoration projects.


Introduction
Peatlands are large carbon stores (Gorham 1991), many of which are damaged by poor land management (Littlewood et al. 2010), reducing their ability to provide a range of ecosystem services (Grand-Clement et al. 2013). With recent alterations to the Kyoto Protocol, carbon trading provides a means of funding peatland restoration projects (Birkin et al. 2011;Dunn and Freeman 2011).
Water level, vegetation composition and subsidence have been proposed as proxies for greenhouse gas emissions (Couwenberg et al. 2011). However, carbon dioxide (CO 2 ) fluxes are affected by water level in some peatlands (Laine et al. 2006;Riutta et al. 2007;Maanavilja et al. 2011) but not all (Urbanov a et al. 2013) and vary between years influenced by environmental conditions (Lafleur et al. 2003). The timing and amplitude of plant growth and senescence are known to influence gas exchange within a range of habitats (Richardson et al. 2013). Measurement of green leaf phenology over whole landscapes is therefore vital in understanding ecosystem carbon dynamics.
Regional canopy development models based on meteorological controls (Jolly et al. 2005) are insensitive to biotic (e.g. predation and disease) and some abiotic (e.g. frost and drought damage) controls. Repeat images from unattended in situ digital cameras or webcams (Richardson et al. 2007;Ahrends et al. 2008;Ide and Oguma 2010;Migliavacca et al. 2011;Sonnentag et al. 2012;Mizunuma et al. 2013;Inoue et al. 2015;Toomey et al. 2015;Wingate et al. 2015) have been used to track vegetation phenology in a range of ecosystems at a fine temporal resolution (daily to sub-daily) averaged over a greater spatial extent (10-100's of metres) than traditional, manual, methods. Vegetation products, based on the relative reflectance of green to red and/or blue in these digital images, are at a temporal and spatial scale where they can be related directly to measured CO 2 fluxes (Botta et al. 2000;Knorr et al. 2010;Migliavacca et al. 2011;Saitoh et al. 2012;Mizunuma et al. 2013;Westergaard-Nielsen et al. 2013;Zhou et al. 2013;Toomey et al. 2015). However, this method is site specific without standard protocols (Ide and Oguma 2010;Migliavacca et al. 2011) and is not straightforward in isolated upland areas such as peatlands.
Satellite data have repeated global coverage at a range of spatial and temporal scales making these a viable tool for monitoring vegetation phenology at regional scales in remote, unmonitored areas. In addition, the Moderate resolution Imaging Spectroradiometer (MODIS) database (2000-present) allows existing CO 2 measurements to be retrospectively related to vegetation products. Studies have shown the potential for using MODIS data combined with meteorological data in light use efficiency models (Schubert et al. 2010;Kross et al. 2013;Yan et al. 2015) as measures of CO 2 exchange. This study builds on such work by assessing a greater range of MODIS products within a maritime peatland of upland Britain. It is hypothesized that (1) MODIS vegetation products (Normalized Difference Vegetation Index (NDVI), Leaf Area Index (LAI8), fraction of absorbed Photosynthetically Active Radiation (fPAR8), Enhanced Vegetation Index (EVI16) and Gross Primary Productivity (GPP8)) can be used to model vegetation phenology as measured by an on-site time-lapse camera and (2) both time-lapse camera and MODIS-derived vegetation products correlate with temporal variations in CO 2 measured in situ (using gross primary productivity measured using a closed chamber method). Our aim was to determine which of these would provide the most useful measure of vegetation phenology for estimating peatland CO 2 budgets.

Digital camera
A battery powered Brinno TLC100 (Brinno Design, CA) time-lapse camera was installed at 0.95 m above the ground-facing magnetic north at an angle of 5°from the horizontal to maximize reflected back-scatter. The field of view was approximately 49.5°facing towards the centre of an experimental catchment with a uniform vegetation cover dominated by Molinea caerulea (Grand-Clement et al. 2014;Gatis et al. 2016;Luscombe et al. 2015aLuscombe et al. ,b, 2016. A daily true-colour image was taken at solar noon (13:00 British Summer Time) to minimize the angular effect of the canopy's bidirectional reflectance function (Ahrends et al. 2008). Between 23-Apr-2013 (DOY 113) and 29-Oct-2013 (DOY 302), 190 images were collected.
The camera records 256 possible brightness values (0-255) for three broad bands (red, green and blue) using an 8-bit charge-coupled device. A region of interest (ROI) was selected (Fig. 2) where the distribution of species was known to be spatially uniform with a dominant cover of Molinia caerulea tussocks. The ROI aimed to maximize spatial sampling of the foreground canopy without giving undue weighting to the closest vegetation. Background vegetation was excluded as it was more likely to have different light conditions, contain a mix of species and be increasingly affected by atmospheric moisture content. The exact areal coverage of the ROI is not known, and is difficult to measure precisely given the orientation of the camera lens coupled with landscape and topographic distortion. However, assuming that the area covered is a rhombus, with the horizontal axes in the foreground being 2 m in length, and at distance, 6 m in length, with an extent of 4 m, we have estimated that the size of this ROI is around 16 m 2 in area. The red and green digital numbers were extracted and averaged across the ROI using MATLAB (R2011b The Mathworks Inc, Natick, Massachusetts). From this the green-red vegetation index (camera-GRVI) (Tucker 1979) was calculated (eq 1).
where G DN and R DN are green and red brightness values, respectively. Chromatic coordinates and green excess index showed greater sensitivity to illumination conditions in initial trials. As blue light has greater day-to-day variation dependent on the quality of incident radiation (Richardson et al. 2007) an index insensitive to blue was selected. Rain, condensation on the lens and uneven illumination adversely affected image quality. A paired t-test was carried out to test the effect of condensation on GRVI. To remove noise, a third-order Fourier adjustment was fitted to the data, points outside the 95% confidence interval were excluded and replaced by interpolation. A third-order Fourier adjustment was chosen above Gaussian (Jonsson and Eklundh 2002) or logistic (Zhang et al. 2003) growth models as camera-GRVI showed a distinct secondary maxima possibly reflecting a decrease in the red band as purple flowers senesced. An 8-day moving average was calculated, selected to match the frequency of MODIS data (Supplementary information (SI) section 2 evidences these data processing steps). It was assumed that the effect of seasonal variation in solar altitude (Ahrends et al. 2008) was negligible.

Terra modis
Terra Earth Orbiting System's Moderate-resolution Imaging Spectroradiometer (MODIS) data were selected as they are freely available and a good compromise between spatial and temporal resolution. Tile h17v3 of MOD09A1 (land surface reflectance), MOD13Q1 (vegetation indices), MOD15A2 (LAI/fPAR) and MOD17A2 (Gross Primary Productivity (GPP8)) products were downloaded (USGS Earth Explorer, http://earthexplorer.usgs.gov) for the period 11-Mar-2013 to 24-Oct-2013. Spatial and temporal resolution details for the datasets are given in Table 1.
Vegetation products at this site with coarser spatial resolution (1000 m) contain mixed land cover content, encompassing both Molinia caerulea-dominated peatland (the target area imaged by the field camera) and areas of improved acidic grassland pasture ( Fig. 1 shows the coarse pixel containing both peatland and rough pasture). MODIS products selected were either estimated from Earth surface reflectance values (Vermote et al. 2002) or derived from these values using an observation operator such as radiative transfer model (Knyazikhin et al. 1999;Running et al. 1999). Radiometrically calibrated data were corrected by NASA for the effects of scattering and absorption by atmospheric gases and aerosols, contamination by thin cirrus clouds and bidirectional reflectance distribution function (Vermote and Vermeulen 1999).
NASA assesses data quality based on cloud cover, observation geometry (solar and satellite zenith angles) and the radiographic testing method was used. MOD9A1 is the best quality observation within the 8-day period. MOD13Q1 is the maximum Normalized Difference Vegetation Index (referred to as NDVI16) and the Enhanced Vegetation Index (EVI16) value within a 16-day period. The acquisition date (MOD9A1 and MOD13Q1), first day of the interval (MOD15A2 and MOD17A2) and data quality data were extracted. Band 1 (near infrared) and Band 2 (red) from MOD9A1 were used to calculate NDVI8, where NDVI8 = (Band 1-Band 2)/(Band 1 + Band 2). GRVI8 (Eq. 1) was also calculated from MOD9A1 data allowing a comparison with camera-GRVI, although one uses physically meaningless digital numbers (dependent on the camera sensor), the other (MODIS) uses surface reflectance; a physical property of the vegetation denoting the ratio of radiance out to irradiance NDVI has widely been used to track spring green-up (e.g. Zhang et al. 2003) as it is known to relate to Leaf Area Index (LAI) and total biomass. EVI was developed to overcome problems associated with NDVI such as saturation at high LAI and sensitivity to variation in background canopy/soil and atmospheric conditions. Both vegetation products have been included in the analysis as NDVI was found by others to be better at predicting the timing of onset of greenness and maturity in deciduous forest (Ahl et al. 2006), but worse than EVI at predicting GPP8 in a northern peatland (Schubert et al. 2010).
Two further 8-day Earth observation products, Leaf Area Index (LAI8) and fraction of absorbed Photosynthetically Active Radiation (fPAR8) were also used in the analysis. LAI is the one-sided green leaf area per unit ground area in broadleaf canopies, whereas fPAR is the fraction of photosynthetically active radiation (400-700 nm) absorbed by green vegetation. The MODIS fPAR8 and LAI8 product is an 8-day composite dataset with 500-m pixel size. The algorithm chooses the "best" pixel available from all the acquisitions of the Terra sensor from within the 8-day period (Myneni et al. 2015). These are known to vary seasonally in Molinia-dominated ecosystems (Nieveen and Jacobs 2002) and are important variables in estimating primary productivity (e.g. Street et al. 2007). LAI8 and fPAR8 (MOD15A2) products are reverse modelled from reflectance values using a biome classification (Knyazikhin et al. 1999), good quality data are then averaged over an 8-day period. The biome assigned was annual grassland, which although not appropriate for most peatlands was suitable for this Molinia caerulea-dominated peatland.
MOD17A2 GPP8 is an 8-day cumulative composite of values calculated using a radiation efficiency use model based on fPAR, LAI, minimum air temperature and vapour pressure deficit (Running et al. 1999). The satellite-derived GPP8 measure has shown potential for monitoring CO 2 fluxes (both GPP and NEE) in northern peatlands (Schubert et al. 2010;Kross et al. 2013).
As blanket peatlands typically form in wet areas (>1000 mm per year) (Lindsay et al. 1988) the probability of cloud cover is high. Therefore, a processing method where questionable data are not immediately omitted, but rather excluded based on an objective assessment against good quality data is required to minimize data loss. Poor quality data (cloudy, high aerosol concentrations or poor geometry) were given a weighting of zero and all other data a weighting of one. To minimize variation due to atmospheric conditions, illumination and observation geometry, a smoothing filter was applied. Logistic functions (Zhang et al. 2003;Fisher et al. 2006;Richardson et al. 2007), high-order polynomials (Cook et al. 2009), asymmetrical Gaussian functions (Jonsson and Eklundh 2002) and second-and third-order Fourier (Fontana et al. 2008) functions were tested. Third-order Fourier series were found to be the best fit (similar to camera-GRVI) as they allowed for the multiple senescence periods observed but did not assume that all data were exact and temporally accurate. Points outside the 99% confidence interval were excluded similar to camera-GRVI but not infilled due to their lower temporal resolution. All remaining points (those inside the 99% confidence interval) were then weighted equally, regardless of quality flag, and a Fourier third-order series fitted. A continuous time series was derived for the MODIS vegetation products (8-and 16-day) using a Fourier third-order series rather than a moving average, as for the camera-GRVI (daily), due to the lower temporal resolution.

Ancillary measurements
The rate of CO 2 accumulation in a closed chamber (55 9 55 cm) over 2 min was measured monthly between May and September (EGM-4 infrared gas analyser; PP Systems, Hitchin, UK), at 100, 60, 40, 10 and 0 % light levels (following Shaver et al. 2007) concurrent with photosynthetic active radiation (PAR) inside the chamber (PAR Quantum, Skye Instruments, Llandrindod Wells, UK). Measurements were taken at 18 Molinia caeruleadominated locations within 500 m of the camera location ( Fig. 1). Net ecosystem exchange (NEE) measurements from all locations were used to derive parameters (P max , K and R Eco ) for a hyperbolic light response curve (Eq. 2) for each sample campaign (monthly). Using these parameters gross primary productivity (GPP) at a PAR of 600 lmol Photons m À2 sec À1 (equivalent to a clear, sunny day) was estimated (eq. 2) for each sample campaign, giving five mean GPP values from 18 replicates across the peatland over the growing season.
where NEE is the net ecosystem exchange (lgC m À2 sec À1 ), R eco the ecosystem respiration (lgC m -2 sec -1 ), GPP is the gross primary productivity at 600 lmol Photons m À2 sec À1 (lgC m À2 sec À1 ), P max the light saturated photosynthetic rate (lgC m À2 sec À1 ), k the half-saturation constant of photosynthesis (lmol Photons m -2 sec À1 ) and PAR the photosynthetic active radiation (lmol Photons m -2 sec À1 ). See Gatis (2015) and Gatis et al. (2016) for more details. This enabled the relationship between mean GPP values and MODIS vegetation products for the mean GPP date to be assessed. A more detailed explanation of the NEE measurement and GPP modelling approach is given in the Supplementary Information (SI) section 1.

Analysis
Four phenological metrics were identified to enable an objective and quantifiable comparison of the timing of seasonal vegetation change; start, peak and end of growth and season length. Several methods have been proposed to identify the start and end of growth; the slope method selects the steepest part of the curve (e.g. Ahrends et al. 2008), the half-maximum selects the value half-way between minimum and maximum values (e.g. Coops et al. 2011), other threshold methods use site-specific percentages on a cumulative frequency curve (e.g. Fontana et al. 2008), with noisy or complex data, visual inspection has been used (Ide and Oguma 2010). The slope method is considered more ecologically meaningful (Zhang et al. 2004) as it does not use arbitrary thresholds and is independent of the amplitude of maximum and minimum values. However, in a direct comparison the half-maximum method was found to more accurately capture traditionally observed phenological events than the slope method (Studer et al. 2007).
The start and end of growth were identified as the DOY when the value first reached half the maximum value, each limb calculated separately. Peak was the DOY when maximum value occurred. The season length is the number of days between the start and end of growth. The average number of days difference between the three phenological metrics (start, peak and end of growth) estimated by camera-GRVI and the MODIS vegetation products were calculated (Days Out).
When GRVI goes from negative winter values to positive summer values (i.e. GRVI = 0) offers another possible metric of phenology (from here on referred to as the zero-threshold method (Motohka et al. 2010)), independent of maximum and minimum values, enabling dates estimated by these two methods to be compared.
To test for correlations between camera-GRVI and MODIS-derived vegetation products, a Pearson's test was applied to data from DOY 116-302. Two-tailed Pearson's tests were also carried out to test for correlations between vegetation products and GPP. A two-tailed approach was followed for testing the significance of the Pearson correlation coefficient.

Field measurements
Digital images from the field camera clearly captured phenological variation in the vegetation community (Fig. 3); through the first appearance of green shoots around DOY 143 (23-May-2013) to mid-summer (Peak Growth) and autumn senescence around DOY 302 (29-Oct-2013). A dramatic increase in green vegetation cover occurred between DOY 156 and 164 coinciding with start of growth by the zero-threshold method. Less obvious changes in vegetation occurred around peak and end of growth. Condensation can be seen affecting image quality (Fig. 3), for example, around DOY 291 and 234. Condensation occurred on the lens on 17 occasions. A paired ttest found no significant difference in GRVI between when condensation was present and the day immediately before (P = 0.857) or after (P = 0.369). Green and red digital numbers ranged between 47 and 204 and 50-213, respectively, indicating neither underexposure nor saturation occurred.
Digital camera-derived Green-Red Vegetation Index (camera-GRVI) increased from À4.9 on DOY 142, coincident with the first visual appearance of green shoots (DOY 143), to reach a maximum of 7.0 on DOY 201 (20-Jul-2013) (Fig. 4). There was a minima around DOY 233 (21-Aug-2013) possibly reflecting an increase in red due to purple flowers then a secondary maxima around DOY 250 as these flowers senesced before decreasing to À9.8 on DOY 302.
Removing outliers then infilling data gaps (Table 2) had no effect on the estimated DOY of peak growth. There was a maximum of 4 and 2 days difference in start and end of growth between the 8-day averaged raw data, data with outliers excluded and camera-GRVI time series used to represent vegetation phenology (outliers infilled). No data processing stage resulted in consistently earlier or later dates than the camera-GRVI time series used. This limited range in estimated threshold dates is indicative of the robustness of the method and confidence in these dates. The 8-day averaged GRVI still exhibited short-term variation (Fig. 4) most likely due to changes in moisture and illumination conditions.

MODIS vegetation products
Visual inspection of the MODIS-derived vegetation products compared to camera-GRVI (Fig. 5A) showed Gross Primary Productivity (GPP8) greened-up earlier and peaked 15 days earlier than camera-GRVI. The magnitude of decrease following peak growth was consistent with   that observed in camera-GRVI, but the autumn maximum visible in camera-GRVI was expressed as a reduction in gradient in GPP8 (discussed in section 4.1). Despite this, comparing across the vegetation products (Fig. 5), GPP8 uniquely captured the asymmetrical shape of the phenology curve as indicated by ground-based photography and had the strongest day-to-day correlation (r = 0.92, Table 3) despite a 15-day offset. Fraction of absorbed photosynthetically active radiation (fPAR8) most accurately estimated peak growth, peaking 1 day later than camera-GRVI (Table 3), before decreasing but not as rapidly as camera-GRVI (Fig. 5C). The gradient increased markedly around DOY 255, reflecting rapid structural decay of leaves, declining to a seasonal minimum consistent with camera-GRVI. NDVI8, NDVI16 and EVI16 showed a similar smooth seasonal increase (Fig. 5D) and decrease, peaking later than fPAR that visually captured the seasonal dynamics of camera-GRVI poorly. Despite poor representation of mid-summer, fPAR8 captured spring green-up and autumn senescence well leading to a strong day-to-day correlation coefficient (r = 0.90, Table 3).
LAI8 was the only vegetation product to show a negative skew, peaking more than 9 days later than the other vegetation products (DOY 223). Despite this because of the similar timing of the rising limb LAI8 showed a strong day today correlation with camera-GRVI (r = 0.89, Table 3).
The 16-day Enhanced Vegetation Index (EVI16) and Normalized Difference Vegetation Index (NDVI8) and MODIS-GRVI8 started to green-up and peaked later than other vegetation indices, except LAI8. EVI16, NDVI16, NDVI8 and fPAR all showed a similar bimodal pattern with maxima occurring near peak growth and the autumn maximum. The relative amplitude of these two peaks is more even in the MODIS products compared to camera-GRVI. The 16-day, 250-m vegetation products (EVI16 and NDVI16) then showed a delayed but more rapid autumn decrease ( Fig. 5E and F) compared to the other vegetation products. Although MODIS-GRVI8 has a similar pattern to EVI16 and NDVI8 in the spring, instead of the secondary maximum observed in other vegetation products and the camera-GRVI there was a late autumn increase (Fig. 5G), resulting in the poorest day-to-day correlation (r = 0.74, Table 3).

Phenological growth metrics
Start and end of growth dates were estimated for camera-GRVI using the half-maximum method and the zero threshold method, these methods estimated start of growth within 2 days (Table 3). All MODIS vegetation products predicted start of growth within an 18-day period (DOY 157-175) with values (Table 3) estimated using 1000-m resolution MODIS vegetation products earlier than those estimated using the half-maximum method by 3-9 days. Vegetation products with finer spatial resolution (<500 m) were later by 2-9 days. The timing of the peak growth varied from DOY 186 (GPP8) to 223 (LAI8), a 37-day period with all vegetation products except GPP8 predicting peak growth later than camera-GRVI (Table 3). End of growth estimated using MODIS vegetation products ranged from DOY 244 (MODIS-GRVI8) to DOY 292 (NDVI16), a spread of 48 days, all later than end of growth estimated by the zero-threshold method (DOY 234). MODIS vegetation products were closer to end of growth estimated by the half-maximum ranging from 9 (LAI8) to 25 (NDVI416) days out. The average days out across the phenological three metrics indicated that fPAR and MODIS-GRVI8 were the best predictors of phenological timings across the whole season based on half maximum and zero threshold, respectively (Table 3).

Comparing vegetation products to gross primary productivity
Gross Primary Productivity (GPP) calculated using model parameters derived from closed chamber NEE measurements increased at the same time as all vegetation products increased indicating greater CO 2 uptake (Fig. 6). The highest modelled GPP value was measured only 4 days before peak growth as determined by camera-GRVI. (Fig. 6I). All vegetation products showed strong correlations (r > 0.75) but only MODIS-GRVI and camera-GRVI significantly correlated with modelled GPP (Fig. 5A and H). NDVI8, NDVI16 and EVI16 did not significantly correlate with modelled GPP (high P -values) as they remained high and relatively constant despite a large range in modelled GPP (Fig. 6). Vegetation indices which best captured the minima around DOY 234 (camera and MODIS-GRVI) showed the strongest correlation with modelled GPP (Figs 5 and 6). Modelled GPP was significantly underestimated by MODIS-GPP8; À41 to À84 lgC m -2 sec À1 compared to À11 to À304 lgC m À2 sec À1 measured.

MODIS vegetation products as measures of phenology
Vegetation products of 1000-m spatial extent (LAI8, GPP8 and fPAR8) predicted start of growth earlier than camera-GRVI and were the best predictors of peak (fPAR8) and end of growth (LAI8 and GPP8). This is most likely as the coarser-grained vegetation products contain a mix Molinia caerulea-dominated peatland and areas of improved acidic grassland pasture (Fig. 1) and will therefore reflect the phenology of both land cover types. Start of growth may therefore be measured earlier in the 1000-m vegetation products as pasture species sprout earlier than M. caerulea which dominates the peatland (Proffitt 1985). However, many species typical of upland pastures (e.g. Common Bent grass (Agrostis capillaris), Sheep's Fescue (Festuca ovina), Sweet Vernal Grass (Anthoxanthum odoratum) and Wavy Hair Grass (Deschampsia flexuosa)) are also present within the Molinia caerulea-dominated peatland, albeit the species dominance is different in the pasture sward compared to the peatland. The good day-to-day correlations and the ability of these 1000-m vegetation products to accurately depict peak and end of growth suggests that once growth had been initiated, the pasture and peatland vegetation responded similar to environmental conditions.
Comparing data across scales (field to landscape) introduces additional uncertainty ). As discussed above, the greater spatial extent of the MODIS data resulted in mixed pixels measuring the seasonality of a vegetation community rather than a single species. Given it is not feasible to ground validate at a scale equivalent to MODIS pixel size, the photographic method provides a larger footprint than traditional manual techniques and is seen by many as a "bridge" to assist scaling across data products (Coops et al. 2011;Hufkens et al. 2012;Klosterman et al. 2014). This study aimed to test whether useful correlations exist between these data sources despite their inherent uncertainties and an order of magnitude difference in the spatial extent being sampled.
This study tested both MODIS vegetation products with an acquisition date (NDVI8, GRVI8, NDVI16 and EVI16) and composite vegetation products averaged over an 8-day interval arbitrarily assigned to the first day of the interval (LAI8, fPAR8 and GPP8). Using the first date of the interval rather than acquisition date for 16-day NDVI has been shown to bias the resultant time series, reducing the accuracy of start of growth estimation (Testa et al. 2014). GPP8, LAI8 and fPAR8 estimated start of growth earlier than camera-GRVI (Table 3) indicating a potential bias in spring, possibly due to using the first date in the interval. However, these vegetation products both over-and underestimated peak and end of growth suggesting the parameters used to derive each vegetation product (e.g. leaf structure, chlorophyll content) are more influential than the choice of date assigned.
For ecosystems where green-up occurs rapidly, such as temperate woodlands (bud burst to canopy maturity in 10-12 days (Ahl et al. 2006)), a higher temporal resolution vegetation product (NDVI8) has been found to better predict start of growth as at least one data point lies within the greening-up period. If there are no data within the period of rapid change, for example, due to cloud cover, the accuracy of start of growth estimated dramatically decreased (Zhang et al. 2009). The deciduous grasses of Exmoor took around 75 days to go from small shoots to maximum green leaf extent (Fig. 5), comparable to the 40-60 days to maximum leaf height observed for alpine grasslands (Migliavacca et al. 2011). Therefore, neither 8day nor 16-day MODIS data are likely to "miss" green-up in these ecosystems where the green-up periods are longer. NDVI16 (250-m spatial resolution) estimated start of growth 2 days later than the half-maximum method suggesting spatial resolution is more important than temporal resolution for accurately predicting start of growth in these fragmented grassland ecosystems.
Vegetation products with 8-day temporal resolution (GPP8, fPAR8 and LAI8) had the strongest day-to-day correlations with camera-GRVI (Table 3). Indicating to capture full-season phenology rather than start of growth, finer temporal resolution becomes more important than spatial resolution.
Green vegetation was observed in the images from DOY 143 to 302. These dates coincide with the base of the rising and falling limbs (Fig. 4) but not the estimated start and end of growth dates. These thresholds (halfmaximum and zero-threshold) were selected to enable an objective and quantifiable comparison between vegetation products rather than for their ability to reflect more traditional phenological metrics. Both the zero-threshold (Motohka et al. 2010) and half-maximum (Fisher and Mustard 2007;Coops et al. 2011;Henneken et al. 2013) methods showed agreement for both start and end of growth (2 and 1 days), excluding the first crossing of the zero threshold (Fig. 5) (DOY 234), indicating their robustness. However, the zero-threshold method can have multiple values if the GRVI fluctuates around zero. In addition, the gradual green-up and senescence in this ecosystem make the timing of start and end of growth estimated using the half-maximum method sensitive to maximum and minimum values. It also explains why there was greater uncertainty in the estimates of end of growth, with a shallower gradient, compared to start of growth within this study (Table 3). Filippa et al. (2015) found more data were required to accurately predict autumn phenology compared to spring phenology.
The range in dates for start, peak and end of growth is progressively larger (Table 3) particularly for higher-order vegetation products such as GPP8, LAI8 and FPAR8 derived using more complex models containing parameters such as chlorophyll and carotenoid content, leaf structure and angle, vapour pressure deficit, etc., that are more variable as the growing season progresses. GPP8 showed the most asymmetry with a small autumnal increase, possibly as meteorological variables are included in the algorithm (Running et al. 1999). GPP8 also showed the earliest end of growth (Fig. 6B, Table 3) reflecting a decline in chlorophyll content, whereas NDVI16, EVI16, NDVI8, fPAR8 showed a delayed but more rapid autumnal decline (Fig. 6D, E, F and G, Table 3) as leaves, which may not have photosynthesized for a number of weeks, finally degraded. These physiological controls on vegetation properties appear to have had a greater impact on the timing of peak and end of growth than either spatial or temporal resolution (Table 3).
MODIS-GRVI8 might be expected to have the strongest correlation with camera-GRVI as it has an 8-day temporal resolution and uses the same wavelengths of light albeit one measures camera-specific digital numbers the other reflectance. However, MODIS-GRVI8 showed an autumnal increase resulting in the weakest day-to-day correlation possibly in part because GRVI8 uses only visible wavelengths where the noise-to-signal ratio is greater than in other bands such as near infrared used by NDVI8. Studies using digital images to track phenology have typically used relative brightness of green (Richardson et al. 2007;Ahrends et al. 2008;Parihar et al. 2013;Alberton et al. 2014), greenness excess index (Ide and Oguma 2010;Coops et al. 2011;Hufkens et al. 2012;Saitoh et al. 2012) or both (Migliavacca et al. 2011;Mizunuma et al. 2013). However, these products are strongly affected by illumination conditions , with the blue band showing most variability. GRVI (Tucker 1979) was selected for this study as it does not use the blue band and is therefore less sensitive to variation in illumination conditions. This is in keeping with Zhou et al. (2013) who found the ratio of green to red better suited to monitor winter wheat phenology than relative brightness of green.
A key consideration when using digital camera imagery to validate MODIS vegetation products is the uncertainty in the digital camera measurement (Keenan et al. 2014). For example, as the elevation of the camera was limited, the viewing angle was low therefore, early green shoots were not observed within the previous year's dead biomass. As a result camera-GRVI started to increase only once shoots protruded above the dead biomass, resulting in a late start of growth date, a similar effect was observed in European grasslands by Wingate et al. (2015). As another example, LAI8 showed rapid spring growth to DOY 180 (Fig. 5B) then continued to increase to a late peak (DOY 238) before rapid senescence, this is consistent with the expected pattern of leaf growth (Nieveen et al. 1998) for Molinia caerulea but in contrast to the phenology measured by camera-GRVI. MODIS-LAI8 uses canopy structural attributes typical of annual grasslands (e.g. no stems, no understory, etc.) in the radiative transfer model for the pixel of interest, appropriate in this deciduous grass-dominated peatland. LAI8 continued to increase beyond maximum camera-GRVI probably as leaf area continued to increase but the canopy had closed in the direction of the camera. Wingate et al. (2015) found greenness to be insensitive to LAI > 1. Without ground validation it is not possible to assess if MODIS-LAI8 accurately measured the dynamics of LAI. However, as both leaf structure and pigment influence photosynthetic activity, it is the combined affect that is of interest.

Vegetation products as proxies for gross primary productivity
It is known that seasonal variation in physiological pigments as well as leaf area, both measurable by digital camera vegetation products, results in variation in photosynthesis (Muraoka and Koizumi 2005;Migliavacca et al. 2011;Saitoh et al. 2012;Mizunuma et al. 2013;Westergaard-Nielsen et al. 2013;Zhou et al. 2013;Toomey et al. 2015). This study investigated if MODIS vegetation products could also identify this seasonal variation.
MODIS-GPP8 significantly underestimated GPP when compared to measured values (Fig. 6B) possibly due to the overly strong influence of temperature and assumptions about the effect of water availability in the MODIS-GPP8 algorithm that are incorrect in a peatland ecosystem. This agrees with direct comparisons between GPP and MODIS-GPP8 in northern peatlands (Moore et al. 2006;Kross et al. 2013). In tropical deciduous woodland MODIS-GPP did not track seasonal dynamics as water availability rather than phenology controlled leaf-on GPP (Christian et al. 2015). Other studies have found MODIS-GPP8 to correlate more strongly with ecosystem respiration (Schubert et al. 2010) as the radiation efficiency use model used to derive MODIS-GPP8 includes air temperature.
NDVI8 and NDVI16 were non-significantly correlated with GPP (r > 0.77, P < 0.128) ( Fig. 6E and F). This reflects both the limited number of samples (n = 5) and aggregated short-term variation in illumination, wetness and temperature in this study. NDVI16 explained 39-71% of the temporal variation in 8-day averaged GPP (Kross et al. 2013) across four northern peatland sites, whereas NDVI16 and EVI16 explained 71% and 44% of the variation in GPP (Rossini et al. 2012) in alpine grasslands. It is likely correlations were lower in these studies as they compared average GPP, including cloudy periods, to MODIS data from cloud-free days, whereas in this study best quality (e.g. fPAR8) or greatest value (e.g. NDVI8) data have been compared to date-matched GPP for the equivalent of a clear, sunny day.
In this study, vegetation indices which best captured the minima about DOY 234, most likely reflecting a period of flowering, had the strongest correlation with GPP (Figs. 5 and 6). Biscoe et al. (1975) found photosynthesis from barley decreased 3 weeks before flowering demonstrating how GPP may be limited by biotic and abiotic factors which do not equally affect vegetation indices. It is not, therefore, proposed that these vegetation indices be used as a direct proxy for GPP, but rather a proxy for vegetation phenology to be used in conjunction with meteorological data to assess GPP. Where meteorological data such as PAR and/or temperature have been used in conjunction with MODIS data the explanatory power has increased markedly (Rossini et al. 2012). At a global scale fPAR has been found to be a key variable for spatially upscaling FLUXNET measurements influencing both seasonal and interannual variability (Jung et al., 2009). The ability to model interannual CO 2 fluxes would enable NEE for a given vegetation community, and thus the potential carbon saved, to be better characterized. In addition, the spatial coverage and temporal resolution of MODIS enables measured fluxes to be spatially and temporally extrapolated making the cost of monitoring any effects due to wetland restoration more practical. The results indicate, despite a small sample size that fPAR8, GPP8 and LAI8 show good day-to-day agreement with camera-GRVI (Table 3) highlighting their potential to be used as phenological proxies in CO 2 models. Bearing in mind the small sample size used in this study we recommend two areas for further work. Firstly, research to determine the best MODIS proxy for vegetation phenology across a wider range of grassland ecosystems. Secondly, determining whether a CO 2 model derived for a Molinia caerulea-dominated site, based on meteorological data and a MODIS vegetation product, is transferable to other sites thereby providing a means of upscaling CO 2 fluxes to landscape and global scales.

Conclusion
Start, peak and end of growth and season length were estimated using the Green-Red Vegetation Index (camera-GRVI) derived from time-lapse camera images. The Moderate resolution Imaging Spectroradiometer (MODIS)-derived vegetation products with the closest start, peak and end of growth dates were NDVI16 (2 days), fPAR8 (-1 day) and LAI8 (-9 days), respectively. Over the three metrics fPAR8 was the most accurate being 9 days out on average. In general, vegetation products with finer spatial resolution (250 and 500 m) more accurately captured spring green-up. Mid-summer and autumnal differences were observed between MODIS vegetation products reflecting their different controls (e.g. chlorophyll content, leaf structure, meteorology, etc.) with GPP decreasing earlier and more gradually than NDVI8, NDVI16 and EVI16. GPP8, fPAR8 and LAI8 had stronger day-to-day correlations (r > 0.89) with camera-GRVI as vegetation products with finer temporal resolution (8 days) better captured whole-season dynamics in this ecosystem. These products were capable of characterizing the phenology of this Molinia caerulea-dominated peatland and therefore offer the most potential as proxies of vegetation phenology at a landscape scale.

Supporting Information
Additional supporting information may be found online in the supporting information tab for this article.
Data S1. 1, Gross Primary Productivity (GPP) Method and 2, Evidence of fit between MODIS and in situ camera data. Table S1. Effect of data processing of Green-Red Vegetation Index (GRVI) thresholds derived from the time-lapse camera. Figure S1. Net ecosystem exchange rates (µg C m sec -1 ) measured (across 18 locations) for a range of photosynthetically active radiation (PAR) (µmol Photos m -1 sec -1 ) and the modelled hyperbolic light-response curve (SI Equation 1). Figure S2. GRVI data processing. Figure S3. MODIS data processing.