Vegetation phenology from Sentinel-2 and field cameras for a Dutch barrier island

(1)

Contents lists available atScienceDirect

Remote Sensing of Environment

journal homepage:www.elsevier.com/locate/rse

Vegetation phenology from Sentinel-2 and

ﬁeld cameras for a Dutch barrier

island

Anton Vrieling

a,⁎

, Michele Meroni

b

, Roshanak Darvishzadeh

a

, Andrew K. Skidmore

a,c

,

Tiejun Wang

a

, Raul Zurita-Milla

a

, Kees Oosterbeek

d

, Brian O'Connor

e

, Marc Paganini

f

a_{University of Twente, Faculty of Geo-information Science and Earth Observation, P.O. Box 217, 7500 AE Enschede, The Netherlands} b_{European Commission, Joint Research Centre, Directorate D - Sustainable Resources, Via E. Fermi 2749, I-21027 Ispra, (VA), Italy} c_{School of Environmental Sciences, Macquarie University, NSW 2019, Australia}

d_{Sovon Dutch Centre for Field Ornithology, Sovon-Texel, Den Burg, The Netherlands}

e_{UN Environment}_{– World Conservation Monitoring Centre, 219 Huntingdon Road, Cambridge CB3 0DL, UK} f_{European Space Agency - ESRIN, Via Galileo Galilei, Casella Postale 64, 00044 Frascati, (RM), Italy}

A R T I C L E I N F O

Keywords: Phenology

Multi-temporal analysis NDVI time series Sentinel-2 Spatial resolution Radiative transfer modelling Landscape variability Salt marsh

Digital repeat photography

A B S T R A C T

Remote sensing studies of vegetation phenology increasingly benefit from freely available satellite imagery acquired with high temporal frequency atfine spatial resolution. Particularly for heterogeneous landscapes this is good news, given the drawback of medium-resolution sensors commonly used for phenology retrieval (e.g., MODIS) to properly represent thefine-scale spatial variability of vegetation types. The Sentinel-2 mission ac-quires spectral data globally at 10 to 60 m resolution everyfive days. To illustrate the mission's potential for studying vegetation phenology, we retrieved phenological parameters for the Dutch barrier island Schiermonnikoog for a full season of Sentinel-2A observations in 2016. Overlapping orbits resulted in two ac-quisitions per 10 days, similar to what is achieved globally since the launch of Sentinel-2B. For eight locations on the island's salt marsh we compared greenness chromatic coordinate (GCC) series derived from digital repeat RGB-cameras with vegetation index series derived from Sentinel-2 (NDVI and GCC). For each series, a double hyperbolic tangent model wasfitted and thresholds were applied to the modelled data to estimate start-, peak-, and end-of-season (SOS/PS/EOS). Variability in Sentinel-2 derived SOS, when taken as the midpoint between minimum and peak NDVI, was well-explained by camera GCC-based SOS (R2_{= 0.74, MSD = 8.0 days,}

RMSD = 13.0 days). However, EOS estimates from camera GCC series were on average almost two months before NDVI-based estimates. This could partially be explained by the observed exponential relationship be-tween GCC and NDVI, as well as by the combined effect of viewing angle differences and the presence of non-photosynthetic elements in the vegetation canopy. A two-layer canopy radiative transfer model incorporating reduced chlorophyll levels in the upper layer provided a physically-based explanation of the viewing angle effect. Finally, we applied the phenology retrieval approach to NDVI series for all pixels of the island in order to map spatial patterns of phenology atfine resolution. Our results demonstrate the potential of the Sentinel-2 mission for providing spatially-detailed retrievals of phenology.

1. Introduction

The timing of periodic events in plants, like budburst orflowering, is generally referred to as plant phenology. Food availability for ani-mals is affected by this timing, and as such plant phenology impacts on animal migration and breeding patterns (Durant et al., 2007; Miles et al., 2017;Shariati Najafabadi et al., 2015). Plant periodic events have a strong relationship with seasonal changes in moisture availability, temperature, and radiation (Stöckli et al., 2008), but are also affected

by extreme weather events (Jentsch et al., 2009). To better understand how climate change may impact on plant and animal populations, ac-curate monitoring of plant phenology in space and time is important (Cleland et al., 2007).

The traditional approach of monitoring plant phenology uses fre-quent visual observations to detect the timing of speciﬁc events like ﬂowering, for example by trained personnel in phenological gardens (Schnelle and Volkert, 1964) or by volunteers (Newman et al., 2012; Schwartz et al., 2012;van Vliet et al., 2014). An alternative approach,

https://doi.org/10.1016/j.rse.2018.03.014

Received 18 July 2017; Received in revised form 16 November 2017; Accepted 11 March 2018

⁎_{Corresponding author at: University of Twente}_{– Faculty ITC, P. O. Box 217, 7500 AE Enschede, The Netherlands.}

E-mail address:a.vrieling@utwente.nl(A. Vrieling).

Remote Sensing of Environment xxx (xxxx) xxx–xxx

(2)

requiring less frequentfield presence, is the use of fixed-position digital cameras that take several photographs of the vegetation each day. This so-called digital repeat photography has been implemented in various networks in Australia (Moore et al., 2016), Europe (Wingate et al., 2015), Japan (Nasahara and Nagai, 2015), and North America (Klosterman et al., 2014). Besides identifying discrete events, analysis of temporal changes in the relative brightness of red, green, and blue (RGB) channels allows for a more continuous tracking of canopy greenness (Richardson et al., 2009). Studies showed that this can be achieved effectively with relatively cheap consumer-grade RGB cam-eras (Nijland et al., 2014;Sonnentag et al., 2012), even if an additional near-infrared spectral band may provide complementary information similar to satellite-derived vegetation indices (Petach et al., 2014). In the past ten years, greenness series from digital repeat photography have been used extensively for assessing phenology and its relationship with primary productivity and climate, predominantly in forests (e.g., Klosterman et al., 2014; Menzel et al., 2015), but also in other eco-systems like alpine grasslands (Migliavacca et al., 2011), tropical sa-vannahs (Alberton et al., 2014), or Arctic vegetation (Anderson et al., 2016).

Even if tower-mounted or mountain slope-viewing cameras may observe larger areas, landscape-scale observation of phenology can only be achieved from satellites. Because this equally requires frequent ob-servation, satellite-based phenology studies mostly relied on medium-to coarse-resolution (250 m medium-to 8 km) optical imagery from sensors like the Moderate Resolution Imaging Spectroradiometer (MODIS) that is acquired at daily intervals (Reed et al., 1994; Vrieling et al., 2011; Zhang et al., 2003). The high frequency is needed to ascertain a su ﬃ-cient number of cloud-free observations throughout the growing season. Since at the medium spatial resolution individual plants cannot be discerned, we generally refer to“land surface phenology” to describe the aggregate temporal behaviour of the multiple plant species and vegetation communities present within a grid cell (de Beurs and Henebry, 2005). Although regional and global patterns of land surface phenology can be eﬀectively captured by these medium-resolution sensors, such sensors may not properly represent the actual phenolo-gical variability, particularly in areas with heterogeneous land cover (Melaas et al., 2013;Vrieling et al., 2017).

Finer-resolution (< 30 m) satellite acquisitions usually come at the expense of longer revisit times, generally resulting in too few cloud-free observations to effectively describe seasonal changes in greenness. One option to obtainfiner-resolution estimates of land surface phenology is to fuse sparse fine- with frequent medium-resolution acquisitions (generally Landsat and MODIS, e.g.,Frantz et al., 2016;Walker et al., 2014;Zhang et al., 2017). Although image fusion has received much attention in recent years, it remains sub-optimal for representing rapid and subtle changes, particularly in heterogeneous landscapes (Emelyanova et al., 2013; Gao et al., 2015; Zhu et al., 2010). An al-ternative is to combine vegetation indices from fine-resolution data, typically Landsat, from multiple years into a single synthetic year (Fisher et al., 2006). From this, multi-annual average estimates of phenological transition dates can be made that can be adjusted annually based on sparser individual-year observations (Melaas et al., 2013; Melaas et al., 2016). This approach has proven to be successful for deciduous forests at locations where two Landsat tiles overlap, but cloud cover around transition dates limited the possibility for annual adjustment of the estimates, particularly for other land covers that have little seasonal variability in greenness and/or strong year-to-year var-iations its temporal behaviour (Nijland et al., 2016; Vrieling et al., 2017).

Taking advantage of new fine-resolution optical sensors with re-duced revisit times, recent studies retrieved phenology directly from fine-resolution imagery without image fusion or combining data from different years. For example,Pan et al. (2015)assessed crop phenology of winter wheat and maize in central China using three years of 30-m resolution data from the two optical HJ-1 (Huan Jing) satellites. For

forest sites in Vermont (United States),White et al. (2014)estimated start of season from combined Landsat TM and ETM+ imagery of two overlapping orbits and compared results withfield estimates of bud burst. For the same site as the present study (Schiermonnikoog, The Netherlands),Vrieling et al. (2017) combined RapidEye and SPOT5 imagery of 2015 to retrieve spring phenology estimates. In that study, comparison with coarser resolution series (e.g., MODIS) demonstrated the additional information content of thefine-resolution retrievals, but ground data to assess the accuracy of the retrievals was not available. This study builds on those efforts and offers a first attempt to re-trieve vegetation phenology from the Multi Spectral Instrument (MSI) of the Sentinel-2 mission. This is thefirst mission that provides free publicly-accessible data at 10-60 m resolution (10 m in four spectral bands) everyfive days with global coverage (Drusch et al., 2012). Al-though thefive-day repeat frequency is only achieved since the launch of Sentinel-2B in March 2017, for Schiermonnikoog we already ob-tained two acquisitions per 10 days in 2016 with Sentinel-2A using two overlapping orbits, thereby more closely simulating the Sentinel-2A and -2B configuration. The goals of this study are (1) to compare retrievals of phenology from digital repeat photography and Sentinel-2A time series for the dynamic salt marshes of Schiermonnikoog (the Nether-lands); (2) to explain differences between both retrievals; and (3) to map phenology at 10 m resolution from Sentinel-2A time series. 2. Study area and data

2.1. Study area

Schiermonnikoog is a barrier island located in the north of the Netherlands between the North Sea and the intertidal Wadden Sea (Fig. 1). About 85% of the 40 km2-large island is a National Park, with the remainder containing grasslands, maize, and built-up area. The National Park contains dunes, marshes, and beaches, which together host a rich fauna and flora. The island's salt marshes are regularly flooded during high tides. Small altitudinal differences affect the fre-quency and duration offlooding and as such the vegetation composition (Olff et al., 1988). Tides, wind, and grazing together make the natural vegetation dynamic and spatially variable. Herbs and grasses of dif-ferent species and height are found across the island, while the dune area in addition contains forest and shrub vegetation.

2.2. Field camera time series

We installed eight Bushnell Trophy Cam Essential (model 119736) trail cameras across the salt marsh of Schiermonnikoog (Fig. 1) that collected 3-megapixel RGB photographs in JPEG format. White balance and exposure are determined automatically by these cameras, and cannot be adapted by the user. Despite that these cameras do not provide calibrated radiance measurements, Sonnentag et al. (2012) showed that phenology-relevant information is retained irrespective of camera choice and JPEG compression. We programmed the cameras to collect ten photos per day, i.e. every half hour between 10:30 and 15:00 CET. Six of the eight cameras (A-F) were placed along a transect from the low marshes to the higher dune area, purposively selecting rela-tively homogenous areas. Camera G was mounted to view a patch of reed vegetation that was not present near the transect. Cameras A–G were installed on 3 March 2016, while camera H was installed already in 2015 at a grassland site with better accessibility. Cameras A–G can only be reached by foot, and the salt marsh area is closed for public during the bird breeding season from 15 April to 15 July, hence limiting interference by humans.Fig. 2provides sample photos for each camera. All cameras faced north to avoid direct incoming sunlight and were mounted on poles at an approximate height of 2.2 m above the ground surface with a depression angle of ~8°, resulting in a detailed view of the area between 10 and 30 m north of the camera location. Metal pins wereﬁxed on top of the poles to avoid disturbance by birds. We noted

(3)

that the camerafield-of-view showed small deviations over time, but because the areas of interest were near the camera this had little effect on the observed area. Image series were collected by replacing memory flash cards on 26 July 2016 and 7 March 2017, although camera H was serviced more frequently. Four cameras did not function as expected for the full time period due to moisture entering the camera (B, E, F) or technical failure after card replacement (G), but nonetheless provided useful data for part of the season. We manually removed any photos that were blurred (due to humidity on the lens), over- or underexposed, or that contained surface water in the scene due to tidal inundation.

For all but one image frames we identified a single area of interest. For camera H two distinct areas were identified including a non-grazed fenced-off foreground area and a background area grazed by cattle. Per photograph and for all pixels within the area of interest we calculated the greenness chromatic coordinate (GCC90C) as follows:

= + + GCC G R G B c (1) where R, G, and B are the digital numbers corresponding to the red, green and blue channels respectively, contained in the JPEG image (Gillespie et al., 1987), and the subscript c refers to the fact that here GCC90Cwas calculated from camera data. Although other indices exist, GCCcis commonly applied to camera time series for phenology studies (Migliavacca et al., 2011) and can well suppress eﬀects of variable scene illumination (Sonnentag et al., 2012). We averaged the GCCc values for each area of interest, and to reduce noise we calculated the 90th percentile of all averaged GCCc values within non-overlapping three-day windows (GCC90c) followingSonnentag et al. (2012).

Fig. 1. The island of Schiermonnikoog with its main natural vegetation communities (Pranger and Tolman, 2012) and agriculturalﬁelds derived from the Dutch public BRP (Basisre-gistratie Gewaspercelen) data set for 2016. Contrary toVrieling et al. (2017), we here include the relatively recent pioneer-zone marsh in the southwest of the island. White areas within the island boundary refer to non-vegetated surfaces (built-up, bare, water) and vegetation communities of limited spatial extent. Letters A to H indicate the locations of the cameras used in this study. The inset shows the location of the island (red) within the Netherlands. (For interpretation of the references to colour in thisﬁgure legend, the reader is referred to the web version of this article.)

Fig. 2. Images acquired by the cameras A-H on 16 June 2016. The red polygon indicates the area of interest used for the extraction of GCC series. For camera H an additional region of interest (blue polygon) is used. (For interpretation of the references to colour in thisﬁgure legend, the reader is referred to the web version of this article.)

(4)

2.3. Sentinel-2 imagery

We downloaded all 33 available (partially) cloud-free acquisitions of Sentinel-2A over Schiermonnikoog between February 2016 and February 2017, i.e. starting and ending in winter with minimal green cover. This resulted on average in 27.6 ( ± 1.7) valid observations per pixel for the island. The acquisitions were made from two different orbits separated by three days, having similar viewing zenith angles (5.2 vs 7.8°) and opposite viewing azimuth angles (104.4 vs 288.6°). We atmospherically corrected the imagery using the Sen2Cor processor (version 2.3.0) to obtain Level-2A surface reflectance products (Louis et al., 2016). For each image we used Sen2Cor's Scene Classification output to discard any pixels classified as saturated or defective, cloud shadow, medium to high cloud probability, thin cirrus, or snow/ice. From the 10 m resolution spectral bands 2, 3, 4, and 8, we then cal-culated NDVIs and GCCs(subscript s for satellite), which were selected here because they represent commonly-used vegetation indices for re-spectively satellite- and camera-based phenology studies. It is noted that while satellite indices are computed with surface reflectance, camera indices are computed with raw digital numbers that scale with reflected radiances.

3. Methods

3.1. Phenology retrieval from index time series

A double hyperbolic tangent function wasﬁtted to all vegetation index (VI) time series, i.e., NDVIs, GCCsand GCC90c. The function is equivalent to the double logistic model (Beck et al., 2006) and can be written as (Meroni et al., 2014):

= + − ∗ + + − ∗ + − t a a t a a a t a a a VI( ) tanh[( ) ] 1 2 tanh[( ) ] 1 2 0 1 2 3 4 5 6 4 (2) where t is time (days). We ﬁtted the function with the Levenberg-Marquardt least squares method as implemented in the IDL-routine MPFIT (Markwardt, 2008), whereby we constrained parameters to re-main within a reasonable range. This routine requires an initial esti-mate of the parameter values. The following list describes the para-meters and indicates in italics how the initial parameter estimates were obtained:

- a0: the minimum VI value. The minimum VI in theﬁrst half of the time series;

- a1[a4]: the VI amplitude of the green-up [senescence] phase. The diﬀerence between the maximum VI and the minimum VI in the ﬁrst [second] half of the time series;

- a2 [a5]: the inﬂection point (days) for the green-up [senescence phase]. The midpoint between the start [end] of the time series and the time of maximum VI;

- a3, a6: control the slope at the inﬂection point for both phases (day−1). Initialized at 0.02 [−0.02].

Because VIs tend to be biased towards lower values due to re-maining atmospheric effects, undetected clouds, and inundation, an iterative procedure was applied tofit the function to the upper envelope of the data following a similar weighting procedure as described in Chen et al. (2004). Wefitted the function only for the support temporal range of 11 February 2016 to 15 February 2017. This corresponded to respectively the first and last Sentinel-2A image used in the satellite series, which was taken as a full year ranging between two minimum (winter) VI values (see alsoFig. 3). Based on thefitted model and the temporal range, four phenological parameters were determined using a threshold approach (White et al., 1997) as follows:

- SOS20: start of season, defined as the first moment when the fitted

function reaches 20% of the amplitude between maximum VI and theﬁtted value for 11 February 2016;

- SOS50: start of season, similar as SOS20but using a 50% threshold; - PS90: peak season, theﬁrst moment when the function reaches 90%

of the amplitude;

- EOS50: end of season, defined as the first moment of the senescence phase when thefitted function reaches 50% of the amplitude be-tween maximum VI and thefitted value for 15 February 2017. Two SOS measures (SOS20 and SOS50) were tested following the various SOS thresholds used in literature, with SOS20representing an early green-up signal and SOS50corresponding to the fastest increase in green-up (White et al., 2014). Rather than using the moment of max-imum VI of the season, we selected PS90as a more stable measure given that thefitted function for satellite VI series frequently resulted in a stable plateau of relatively large VI values during summer that persisted between one andfive months.

This retrieval procedure was applied to NDVIsand GCCsseries of each Sentinel-2 pixel within the island boundaries (Fig. 1), including those Sentinel-2 pixels corresponding to the areas of interest observed by theﬁeld cameras, as well as to the ﬁeld camera GCC90cseries. For two cameras (B and G) no data were acquired after August 2016. Be-cause only the green-up phase could be modelled for those cases, we applied a single hyperbolic tangent model as inVrieling et al. (2017) taking as the support temporal window 11 February to 1 August 2016. In this way we could still use the data from those cameras to estimate SOS20, SOS50, and PS90.

3.2. Radiative transfer modelling

We hypothesized that satellite and camera (GCCsand GCC90c) time series may have diﬀerent temporal trajectories due to:

1) the diﬀerence in view zenith angle only (VZA, being near-nadir for the satellite and about 70° for the camera);

2) the different VZA in combination with the observed presence of dry material during the green-up phase (i.e. the dead phytomass of the previous year) and of non-photosynthetic organs (i.e.flowers and dry reproductive organs) after maximum leaf development; 3) the difference in spectral sensitivity between the Sentinel-2A MSI

and the Bushnell cameras.

With regard to the third hypothesis, the actual response of the Bushnell cameras could not be reproduced because the exposure and white balancing for this type of cameras are fixed to automatic. In addition the camera spectral response function for the RGB channels was not available. Therefore, we restricted the analysis to evaluate what Sentinel-2A's MSI would observe from the camera viewpoint. We did not specifically address possible differences in sun elevation angles, because a) the satellite acquisition time (i.e., between 11:40 and 11:56 CET) was within the daily observation window of thefield camera, b) we found no evidence of daily GCCctrends, but rather found that cloud-induced illumination differences had the strongest effect on daily GCCc variation, and c) the GCC90cmeasure results from various observations and is not related to a single timing.

To test thefirst two hypotheses we used the two-layer canopy bi-directional radiative transfer (RT) model 4SAIL2 (Verhoef and Bach, 2003) coupled with the leaf optical properties model PROSPECT5B (Féret et al., 2008; Jacquemoud and Baret, 1990) that simulates the reflectance and the transmittance of a leaf as a function of its con-stituents. The 4SAIL2 model follows a four-stream concept and it is the two-layer version of the model SAILH (Verhoef, 1984), a one-dimen-sional bidirectional turbid medium radiative transfer model that also describes the hot spot effect in plant canopy reflectance (Kuusk, 1991). We simulated the temporal evolution of GCC90Cfor the nine locations corresponding to areas of interest defined on the camera frames using

(5)

diﬀerent assumptions on RT set-up and model parameters in an attempt to provide a possible physical explanation for the mismatch.

For each camera location we ran three simulation sets with different parameterization intended to represent the two working hypotheses: 1) a canopy described by a single layer with homogeneous leaf and canopy properties (RT-1layer); 2) an addition of a second top layer with re-duced chlorophyll content that reflects the dominance of dry vegetation carried over from the previous year and the appearance offlowers and grass stems later in the season (RT-2layer); and 3) a further reduction of the chlorophyll content of the upper layer to represent the continuing presence of non-photosynthetic elements during the season (RT-2layer*). Simulations 2 and 3 differ by the value of the maximum chlorophyll content of the upper layer. With RT-1layer we tested if the different view angle alone could explain the observed differences when considering the development of a simple and homogeneous green ca-nopy. With RT-2layer and RT-2layer* we still tested the viewing angle effect, but attempted to provide a more realistic description of the ca-nopy development as depicted by the camera image time series. Camera pictures show that relatively tall and dry phytomass material dominant at season start is progressively less visible as vegetation grows. Close to the timing of maximum biomass expansion,flowering, stem elongation, and seed formation take place above the canopy plane. This in-florescence progressively dries out during autumn. To represent such dynamics in our simplified RT modelling scheme, we modelled a lower layer with spherical leaf angle distribution and an upper layer occu-pying 20% of the total height with erectophile leaf angle distribution. While chlorophyll content was kept constant for the lower layer, we initiated the chlorophyll content of the upper layer at 0 for February

2016 and increased it in proportion to LAI to reach the same level as the lower layer at the time of maximum LAI to mimic the gradual re-placement of dead phytomass with new green vegetation. After max-imum LAI we reduced the upper level chlorophyll using an exponential decay function with a half-life of 30 days to represent the rapid increase of non-photosynthetic material related to flowering and seed/stem formation. For RT-2layer* we further reduced the maximum chlorophyll level of the top layer by 25%. We acknowledge that this modelling set-up can only be a rough representation of the real vegetation dynamics, which in reality also vary between camera locations. However, our main purpose was not to construct a model that accurately mimics the temporal dynamics of the observed GCC time series, but rather to assess if the observed differences in retrieved phenological parameters (i.e., due to viewing angle differences between satellite and field camera) can be reproduced with a physically-based explanation.

Canopy reflectance is computed for each available Sentinel-2A overpass using the illumination geometry (i.e. sun zenith and azimuth angles) at time of overpass. Satellite VZA isfixed to 0° (i.e. nadir-view) and the camera VZA to 70° (i.e., the average view angle of the camera ROI). The camera view azimuth angle was set to 0° as all cameras face North. In all simulations, the background reflectance is considered Lambertian and set to the average of three measurements acquired in situ with an ASD-Fieldspec 3 and a nadir view. To represent a plausible temporal vegetation development at the various sites we used the LAI retrieved from radiative transfer model inversion on Sentinel-2A data over the period of interest (O'Connor et al., 2017). Other model para-meters were set to the camera site mean retrieved value. A summary of the used parametrization is reported inTable 1.

Fig. 3. Time series of NDVIsand GCCsderived from Sentinel-2A, and GCCcand GCC90cderived from theﬁeld cameras. Each panel relates to the corresponding camera identiﬁer. For

camera H two areas of interest are defined, corresponding to panels h1 and h2 in this figure. The fitted double hyperbolic tangent functions are plotted for each series (single for camera series E and G), and the corresponding retrievals for SOS20, SOS50, PS90, and EOS50are indicated as vertical lines at the bottom of each panel.

(6)

Following the simulations, the spectral simulation outputs (at 1 nm interval) were resampled to Sentinel-2A's bands 2 (blue), 3 (green), and 4 (red) using the spectral response functions of the instrument (European Space Agency, 2017). Based on the simulated reﬂectances for these bands, we then calculated GCC for both VZA's. Finally, we applied the same procedure as in Section 3.1 to retrieve phenological parameters from the simulated GCC90Ctime series.

3.3. Statistical analyses

To assess if thefitted function (Eq.2) closely matched observations, we assessed the modelfit using the Pearson correlation coefficient (r) between the observed andfitted vegetation index values for each time series, as well as the corresponding root mean square deviation (RMSD) and mean signed deviation (MSD).

To compare phenological retrievals fromfield camera and Sentinel-2 series, we calculated the R2, RMSD, and MSD for each phenological parameter. The same was done for the comparison of the phenology estimates from simulated GCC data at different viewing angles. To evaluate if the difference between camera- and satellite-derived phe-nological retrievals can be reproduced by RT-model simulations, we focused on the MSD between the phenological timing derived from the two viewing geometries. In other words, we were interested in testing if the observed shift between camera- and satellite-derived timings could be explained by the different modelling set-ups. For this, differences in MSD were tested by ANOVA Randomized Block Design and Tukey's honest significant difference (HSD) test using the conventional 0.05 significance level.

To assess to what extent the selected VI may impact the phenology retrievals, we empirically modelled the relationship between Sentinel-2 derived GCCsand NDVIs. For this, an exponential function was selected after visual inspection of the scatterplot of the two variables. We also calculated the root mean square error (RMSE) and mean signed error (MSE) for this relationship; because the dependent variables are ob-served rather than estimated values, we replace‘deviation’ with ‘error’ here.

Finally, we assessed how reduced Sentinel-2 image availability may impact the retrieval of phenological parameters. Weﬁrst split the data into acquisitions by relative orbit 8 and 51 to examine whether ac-quisitions from a single orbit would have resulted in similar retrievals. Subsequently, we simulated reduced image availability due to periodic persistency of cloud cover by iteratively excluding all observations from

a single month. We focussed this analysis on Sentinel-2 derived NDVI only and limited it to the pixels corresponding to the nine areas-of-interest observed by theﬁeld cameras, as deﬁned inFigs. 1 and 2. We calculated the MSD and RMSD to summarize deviations between re-trievals from the reduced data set versus the full data set.

4. Results

4.1. Phenological retrievals: Camera versus satellite

Temporal profiles of camera-derived GCCc and GCC90c, and sa-tellite-derived NDVIsand GCCsare displayed inFig. 3, together with the fitted functions and the four extracted phenology parameters for each series. GCCcis plotted on the secondary y-axis as it is smaller than GCCs. The cause of this difference can be that GCCsis computed from re-flectance values while GCCcis derived from raw digital number pro-portional to reflected radiance. In addition, cameras do apply an au-tomatic scaling (i.e. the white balance) to the RGB channels.

Fig. 3illustrates that Sentinel-2 observations are reasonably spread throughout the year. Only three observation gaps of more than a month existed; 40 days between 7 June and 17 July, 40 days between 25 Oc-tober and 4 December (for camera H, slightly shorter for others), and 43 days between 4 December and 16 January. The temporal behaviour could be accurately modelled with Eq. (2) resulting in an average Pearson correlation coefficient (r) of 0.88 for NDVIs and 0.89 for GCCs betweenfitted and original values. For the GCC90cseries, r was 0.95. ForFigure 3h1 all indices show a decrease in greenness in July/August followed by a subsequent increase; this is caused by intrusion of cattle in the normally fenced-off area resulting in a trampling and reduction of green grass cover followed by regrowth. Despite this“intervention” we fitted the same model to the series.

Fig. 3shows that the green-up phase is more similarly timed be-tween camera and satellite series, as compared to the reduction of greenness during or after summer. The GCC90cseries generally display an earlier and faster decrease in greenness as compared to the satellite series. This may be partially attributed to the development of non-photosynthetic elements like flowers, grass stems, and seed heads, which are dominant in the oblique view of thefield camera but less important in nadir view because of their predominantly vertical elon-gation and higher position in the canopy.Fig. 4illustrates this with selected photographs for cameras B, C, and D. The numbers inFig. 4 show that NDVIs remains above 80% (of the green-up amplitude) throughout September for the selected camera locations, while GCC90c reduces much quicker. For camera B this is mostly due to the presence of common sea-lavender (Limonium vulgare), for which the purple flowering peaks in August but its flowers remain dominant on the photograph also after fading. For camera C, sea couch (Elytrigia atherica) grass forms stems with vertically-oriented seed heads that elongate and gradually lose their green colour from August onwards. Camera D contains a mixture of the effects of grass stems and heads (Festuca rubra), andflowering (Juncus gerardii and Limonium vulgare).

Phenological parameters retrieved from GCC90cseries are plotted against their equivalents retrieved from the NDVIsand GCCstime series in Fig. 5. On average, satellite retrievals are earlier than camera re-trievals for SOS20(MSD = 19 days for NDVIsand 1.7 days for GCCs), and gradually become later than the camera retrievals going from SOS50 via PS90 to EOS50. For NDVIs-retrievals the RMSD is lowest (13.0 days) for SOS50, for which the satellite retrievals explain 74% of the variability of camera retrievals. For GCCs, SOS20is the retrieved parameter that corresponds closest to camera retrievals (R2= 0.72, RMSD = 6.7 days). NDVIs-based retrievals are generally earlier than GCCs-based retrievals by 9 to 18 days, except for EOS50for which the GCCs-based retrieval is about a month earlier on average (MSD =−26.6 days for GCCs, −55.0 days for NDVIs). The largest RMSD is also found for EOS50(61.6 days for NDVIs, 35.0 days for GCCs), in line with our earlier observations (Figs. 3 and 4).

Table 1

RT model parametrization.

Parameter Value

Leaf (PROSPECT5B)

Chlorophyll RT-1layer: 35μg cm−2(RT-1layer), RT-2layer: scaled between 0 and 35μg cm−2

for upper layer (RT-2layer).

RT-2layer*: reduced by 25% with respect to RT-2layer.

Leaf mass per unit area 0.01 g cm−2 Equivalent water thickness 0.01 mm PROSPECT structure coeﬃcient 2.09 Canopy (4SAIL2)

LAI Variable in time and per camera location Leaf angle distribution lower

layer

Spherical

Leaf angle distribution upper layer

Erectophile

Hot spot parameter 0.1

Fraction of lower layer 1.0 in 1layer, 0.8 in 2layer and RT-2layer*

Dissociation 1.0 Clumping parameters Cv and

zeta

(7)

We found that GCCshas an exponential relationship with NDVIs (Fig. 6) that does not vary over time. GCCsis more sensitive to diﬀer-ences in NDVIsfor large NDVIs values (i.e., larger green biomass) as compared to smaller values. To evaluate if the exponential relationship can explain part of the diﬀerences between NDVI- and GCC-derived phenological parameters, we transformed NDVIsinto GCCs* using the Equation inFig. 6a.Table 2indicates that the transformation indeed makes phenological retrievals more similar to those derived from GCCs. For each phenological parameter the R2 _{increases, while MSD and} RMSD decrease. This demonstrates that a non-linear relationship

between VIs is translated into a diﬀerent shape of the VI temporal evolution and hence of the retrieved phenological parameters. 4.2. Viewing angle eﬀects from radiative transfer modelling

Fig. 7illustrates for two camera locations (camera C and D) the simulated GCCc and GCCsseries and theirﬁtted proﬁles for the RT-1layer, RT-2layer and RT-2layer* models. Assuming that the LAI re-trievals are accurate, the RT-1layer model resulted in a poor description of the GCC temporal pattern (Fig. 7a, d), compared to the observed

Fig. 4. Sample photographs for cameras B, C, and D around day 4 ( ± 2) of each month. The numbers on each photo indicated theﬁtted value for camera-derived GCC90 (top) and satellite-derived NDVI (bottom), expressed as a percentage of the amplitude between the maximum VI and theﬁtted value for 15 February 2017.

Fig. 5. SOS20(a), SOS50(b), PS90(c), and EOS50(d) derived from the GCC camera series on the x-axis plotted against the retrievals from Sentinel-2 NDVI and GCC series on the y-axis. The

dashed black line shows the 1:1 relationship. A positive MSD means that the satellite-retrieval is earlier than the camera-retrieval.

(8)

values in Fig. 3c and d. The addition of a top layer with varying chlorophyll levels over the season (Fig. 7b, c, e, f) produced more realistic patterns. As expected, these still do not match perfectly with observations because of our rough assumptions about depth of the top layer and chlorophyll dynamics that are uniformly applied to all loca-tions and because of inherent difficulties in modelling an actual canopy with a turbid medium approach. Nonetheless, the model results illus-trate that the addition of this less-photosynthetic top layer offers a physically sound explanation for the earlier EOS50observed for GCCcas compared to GCCs. For the green-up phase, retrievals from observed of SOS20, SOS50, and PS90 are earlier for GCCcthan those for GCCsfor camera locations C and D (Fig. 3c, d). While the sign of this difference is reflected in the retrieved SOS20, SOS50, and PS90from the two-layer

Fig. 6. Relationship between NDVI and GCC derived from Sentinel-2 time series for (a) pixels corresponding to camera locations, and (b) all pixels on the island.

Table 2

Comparison of phenological retrievals from GCCsseries with NDVIsand GCCs* series for

the camera locations. GCCs* refers to the transformed NDVIsseries applying the empirical

relationship ofFig. 6a.

SOS20 SOS50 PS90 EOS50

NDVI GCC* NDVI GCC* NDVI GCC* NDVI GCC*

R2 _0.54 _0.69 _0.83 _0.87 _0.60 _0.86 _0.77 _0.82

MSD 17.33 −5.44 15.11 −5.22 9.00 −6.33 −29.89 3.56 RMSD 19.47 13.68 16.81 10.62 16.50 9.99 33.04 12.35

Fig. 7. Temporal profiles of GCC from a nadir satellite perspective and the oblique camera perspective derived from the radiative transfer modelling for the locations of camera C (a, b, c) and D (d, e, f). Panels a and d show results for RT-1layer, panels b and e for RT-2layer, and panels c and e for RT-2layer*. For GCC, the lines show thefitted double hyperbolic tangent model. In addition, the simulation input values for LAI and Chl are shown. For RT-2layer simulation panels, the Chl profile refers to the Chl in the upper layer, Chl in the lower layer being kept constant as in panels a and d.

(9)

models, the magnitude of the diﬀerences vary. The further reduction of chlorophyll levels in the RT-2layer* model resulted only in small changes for the retrieved phenological parameters.

Combining all camera locations,Table 3summarizes the differences between GCCc- and GCCs-based phenological retrievals for observed and simulated data. It confirms that the RT-1layer model could not reproduce the differences of the observed data; MSD between camera-and satellite-view retrievals deviate by at least 20 days compared to observations. For the two-layer models, the match between MSD values is much better, particularly for the 50%-threshold parameters (SOS50 and EOS50). For example, for these parameters the MSD and RMSD for observed and RT-2layer* match within a day. For all phenological parameters, the MSD of RT-1layer is significantly different from the MSD of observed data, RT-2layer and RT-2layer*, based on Tukey's HSD test. On the contrary, and again for all phenological parameters, the MSD of RT-2layer* is not significantly different from the MSD of ob-served data. The same is true for the MSD of RT-2layer with the ex-ception of SOS20. This indicates that even though we used a single and simple representation of the top layer's evolution, it can explain a large part of the variability between phenological retrievals derived from nadir versus oblique viewing angles.

4.3. Mapping of phenological parameters

Fig. 8 shows the result of applying the phenological retrieval ap-proach to all individual Sentinel-2 pixels within the island boundaries. For 77.5% of the pixels the retrieval was successful, i.e. it could be

initiated and converged. Particularly for classes with limited NDVI variability throughout the year (Kpp and Sv) retrieval of phenological parameters often failed, and also the Pearson correlation coeﬃcient for successful retrievals is smallest of all classes (Table 4). On average for the island, SOS20is reached on 20 April, SOS50on 12 May, PS90on 13 June, and EOS50on 16 November. However, there is a large spatial variability both within and between vegetation communities (Table 4). Average dates for SOS20, SOS50, and PS90 increase going from the middle-high marshes (Km) to the pre-pioneer zone marshes (Kpp). Maize behaves distinct from the other classes; it has the smallest stan-dard deviation for SOS50and PS90of all classes, a late green-up and the earliest EOS50 due to the ripening and harvest early October. As a consequence, all maizeﬁelds are clearly visible inFig. 8.

The difference between phenology retrievals from NDVIsand GCCs series is shown inFig. 9. The maps confirm our earlier findings (e.g., Table 2) that SOS20, SOS50, and PS90are generally earlier for NDVIs-based retrievals (on average by 7.0, 7.3, and 6.7 days respectively), while EOS50is earlier for retrievals from GCCs(5.4 days). Substantial spatial variability exists between these differences. NDVI amplitude could explain a small fraction of this variability with significant Pearson correlation coefficients between green-up amplitude and SOS20, SOS50, and PS90of−0.230, −0.296, and − 0.236 respectively.

4.4. Eﬀects of reduced image availability

Table 5summarizes how reduced image availability would have affected the retrievals from Sentinel-2 derived NDVI for the nine camera locations. The largest effects are found when only acquisitions from a single relative orbit are used, because this reduces data avail-ability to about half of the full data set. The relative large RMSD values suggest that single-orbit acquisitions for Schiermonnikoog miss out important elements of the changes in greenness, and are therefore in-sufficient to accurately retrieve phenology. Excluding all acquisitions of a single month to mimic persistent cloud cover had varying effects, depending on which month was excluded. RMSD values were largest when February imagery (both at start and end of the period) was omitted, and second largest for SOS20and SOS50when removing August imagery. This underlines the importance for phenology studies of image availability not only during vegetation growth, but also during winter (minimum NDVI) and peak season. Instead for EOS50, removing a single

Table 3

MSD and RMSD between phenological parameters retrieved from observed series of sa-tellite GCCsandﬁeld camera GCCc; simulated GCCsand GCCcderived from a one-layer

radiative transfer model, a two-layer model, and a two-layer model with reduced chlor-ophyll concentration in the upper layer.

SOS20 SOS50 PS90 EOS50

MSD RMSD MSD RMSD MSD RMSD MSD RMSD

observed 1.7 6.7 −7.1 10.5 −16.4 21.6 −26.6 35.0 RT-1layer −19.7 23.6 −28.3 32.7 −41.3 51.5 15.6 19.5 RT-2layer −9.4 10.9 −10.0 11.7 −10.3 14.6 −30.8 37.6 RT-2layer* −8.2 9.8 −7.3 9.5 −5.8 10.7 −27.0 34.3

Fig. 8. Retrievals of SOS20(a), SOS50(b), PS90(c), and EOS50(d) from Sentinel-2 NDVI series for the island of Schiermonnikoog. The outer circle of the legend represents theﬁrst 10-day

period of each month.

(10)

data point in November or December resulted in a relatively large RMSD. Excluding images in spring months (March to June) resulted in RMSD values for SOS-measures of < 5.0, suggesting that if part of the spring season would have above-normal cloud cover, the eﬀect on SOS remains within limits. This can be attributed to the fact that the green-up phase (from minimum to maximum NDVI) usually takes place over a period of at least two months in the natural salt marsh areas (Fig. 3). However, when green-up is more rapid, like on agricultural plots or in

semi-arid rangeland systems, missing a month of data may have a stronger eﬀect.

5. Discussion

Our analysis showed that the retrieval of vegetation phenology at ﬁne spatial resolution is possible with frequent Sentinel-2 imagery. Our NDVIs-based retrievals were on average within 10 days of in-situ GCCc

Table 4

Summary statistics of the phenological retrievals from Sentinel-2 NDVI series for the entire island and per vegetation community (Fig. 1). Mean ± standard deviation are reported for the phenological parameters. Fit statistics (Pearson r, MSD, and RMSD) are calculated betweenfitted and original NDVI values; #images indicate the average number of Sentinel-2 observations used perfit, and %pixels specifies for what percentage of pixels within each class, a model could successfully befitted. Per column, values are colour-coded from light (small) to dark (large).

%pixels 94.5 92.1 93.8 95.8 93.7 98.5 96.4 72.2 28.5 99.8 99.9 1.1 77.5 #images 28.1 28.2 27.8 27.7 28.0 27.8 27.3 25.4 22.9 27.1 26.8 27.6 27.6 RMSD 0.073 0.060 0.072 0.062 0.062 0.068 0.071 0.082 0.068 0.083 0.109 0.053 0.072 MSD 0.017 0.013 0.017 0.015 0.014 0.015 0.017 0.020 0.017 0.021 0.024 0.014 0.017 r 0.84 0.91 0.88 0.92 0.91 0.92 0.88 0.82 0.80 0.84 0.92 0.80 0.87 EOS50 316 ± 22 327 ± 12 308 ± 22 312 ± 23 329 ± 14 325 ± 15 318 ± 17 338 ± 26 352 ± 31 341 ± 18 284 ± 15 357 ± 30 321 ± 24 PS90 157 ± 27 169 ± 29 157 ± 22 170 ± 20 169 ± 29 158 ± 21 174 ± 22 198 ± 38 212 ± 38 159 ± 20 200 ± 6 238 ± 71 165 ± 30 SOS50 127 ± 19 130 ± 25 132 ± 17 136 ± 18 131 ± 24 124 ± 17 137 ± 21 159 ± 41 178 ± 38 137 ± 22 185 ± 13 198 ± 66 133 ± 26 SOS20 107 ± 19 104 ± 24 114 ± 18 112 ± 19 106 ± 23 102 ± 16 110 ± 22 132 ± 47 152 ± 47 123 ± 27 174 ± 21 166 ± 67 111 ± 28 Vegetation community Dd: dry dunes Ddk: dune mosaick Dv: humid dunes Kb: brackish marsh Kh: high marsh Km: middle-high marsh Kl: low marsh Kp: pioneer-zone marsh Kpp: pre-pioneer-zone marsh Gr: grassland Ma: maize

Sv: beach and embryo dunes All pixels island

Fig. 9. Diﬀerence (in days) between phenology parameters derived from Sentinel-2 NDVI and GCC series for SOS20(a), SOS50(b), PS90(c), and EOS50(d). Negatives mean that the NDVI

(11)

retrievals of SOS50, and PS90, while early green-up (SOS20) and senes-cence (EOS50) had an MSD of 19.0 and− 55.0 days, respectively (Fig. 5). The better correspondence between camera andﬁne-resolution satellite-retrievals for green-up than for senescence is in line with phenology studies in forest systems. For example, using multi-year Landsat analysisNijland et al. (2016)reported an RMSE of 7 days for SOS50and 14 days for EOS50for mixed and coniferous forests in the Rocky Mountains. Similar ﬁndings were obtained by Melaas et al. (2016)for deciduous forest sites in North America.

Part of the deviation between phenology retrievals from NDVIsand GCCc can be explained by the use of different spectral indices. The exponential relationship (Fig. 6) makes that GCCsincreases later in the green-up phase and stabilizes earlier in the decay phase, as compared to NDVIs. GCCs-based phenology retrievals confirm these expected timing differences, but result only for SOS20and EOS50in a slightly stronger agreement with GCCc-based retrievals (Fig. 5). The exponential re-lationship also suggests that GCC is more sensitive to higher biomass levels. This empiricalfinding is opposite toWingate et al. (2015)who found that GCC was insensitive to LAI values above 2 based on a sen-sitivity analysis with PROSAIL, whereas NDVI remained sensitive for a wide range of LAI values.

An important difference between field camera and satellite ob-servations is the viewing angle. Particularly the presence of a non-photosynthetic top layer (flowers, stems, seed heads) may reduce the visibility of green vegetation from an oblique view. Contrary to the proverb, we found that grass is not always greener on the other side of the fence. This effect applies to various periods in the phenological cycle (Moore et al., 1991), including during early vegetative growth due to presence of dead grasses and stems from the previous year (e.g., Liao et al., 2008;Ter Heerdt et al., 1991). Even if the top layer would be green, reflection in the visible bands is much affected at off-nadir viewing angles, as demonstrated for example for green wheat panicles (Zipoli and Grifont, 1994). As such, the sensor viewing angle affects spectral vegetation indices (e.g.,Deering et al., 1992) and their tem-poral trajectory. The results of our two-layer radiative transfer model confirmed this dependency, even if the multiple vegetation components could at best be roughly represented in these models (see alsoWenhan, 1993). By adding a top layer of vegetation with reduced chlorophyll content in our model, we could further explain observed differences

between phenological parameters when retrieved from the camera versus satellite viewing geometry.

For the relatively inexpensive Bushnell cameras the spectral re-sponse function is not available and we could not adjust the automatic white balance and exposure. The use of more advanced camera systems for which the spectral response functions are known, could help in better characterizing the camera bands as input to radiative transfer modelling (Wingate et al., 2015). However, we did attempt replacing the Sentinel-2 spectral sensitivity with that from a Nikon RGB-camera, but found only minor differences in simulated GCCcvalues, suggesting little effect on the phenological results (data not shown). Their low cost allowed for installing multiple cameras in publicly-accessible sites without the need for additional security measures. Various studies have shown that cameras with different spectral responses and white balance settings resulted in similar phenological metrics (e.g.,Mizunuma et al., 2013;Sonnentag et al., 2012;Zhao et al., 2012). Our study did high-light that viewing angle has an effect on the retrieval of phenological metrics. In this regard, camera systems observing at nadir could provide a stronger link with satellite recordings, although their observed area would be smaller compared to a tilted camera mounted at the same height above the canopy.

Despite differences between our satellite- and camera-based phe-nological retrievals, the density of Sentinel-2A observations allows ex-traction of phenological parameters with a consistent approach across the island atfine spatial resolution. Application of a consistent meth-odology can provide useful information on spatial patterns and inter-annual trends of phenology, even if literature suggests that phenolo-gical retrievals from satellite data remain typically within 10–30 days from visually-observed transition dates of individual plant species (White et al., 2014). This can partially be attributed to the occurrence of multiple species within a grid cell resulting in mixed spectra (Vrieling et al., 2017;Zhang et al., 2017). Thefiner spatial resolution of Sentinel-2 reduces this effect, as compared to medium-resolution series like MODIS. Nonetheless, most of the salt marsh is floristically diverse (Bockelmann et al., 2002;Pranger and Tolman, 2012), which implies that phenological metrics for 10 m grid cells still relate to the phenology of the vegetation community rather than of individual plants. Instead, for more uniform areas like cropland, Sentinel-2 retrievals can be re-lated to the phenology of specific crop species, such as the timing of

Table 5

Impact of reduced image availability on NDVI-based phenology retrievals for the camera locations as compared to the full data set; #images indicates the average number of Sentinel-2 observations used perﬁt. The rows with the relative orbits indicate that only observations from that orbit were considered. The rows with individual months indicate that all observations for the indicated month were excluded from the retrieval. RMSD is colour-coded based on the overall minimum and maximum values (same for #images); MSD is colour-coded with blues indicating positive values (retrieved date for reduced data set is later than for the full data set) and red negatives.

#images MSD RMSD

Reduced data set SOS20 SOS50 PS90 EOS50 SOS20 SOS50 PS90 EOS50

Relative orbit 8 14.3 5.7 –0.3 –6.9 13.9 13.1 9.4 10.4 21.7 Relative orbit 51 13.4 2.6 –1.5 –3.8 –18.9 12.6 11.5 13.4 20.9 January 26.9 –0.3 0.0 0.3 –2.2 0.6 0.5 0.7 4.0 February 24.9 –0.8 –0.8 –0.7 2.3 9.6 5.1 2.6 12.0 March 24.7 –1.2 –0.1 1.2 –0.4 4.7 3.0 2.7 1.1 April 23.1 0.3 0.6 0.4 –0.4 4.4 2.6 2.9 0.7 May 24.7 0.3 0.7 0.4 –0.4 2.9 3.4 6.5 0.7 June 26.7 –0.4 0.0 0.7 –1.0 3.1 1.6 2.8 2.1 July 24.7 0.9 –0.7 –2.4 –0.1 1.8 2.2 5.5 3.4 August 25.4 3.4 2.2 0.3 3.0 8.1 4.4 3.4 5.6 September 24.1 –0.2 –0.1 –0.1 1.1 1.1 0.9 2.0 3.3 October 25.9 –0.6 0.1 0.4 –1.0 0.9 0.7 1.2 5.1 November 26.7 –0.3 0.3 0.6 2.8 0.7 0.9 1.6 8.6 December 26.7 –0.2 –0.6 –0.6 –2.2 1.2 1.0 2.0 6.7

(12)

crop emergence (Claverie et al., 2012;Gao et al., 2017).

The phenological retrievals for the entire island revealed a large spatial variability (Fig. 8), confirming our earlier findings with 2015 imagery of RapidEye and SPOT5 (Vrieling et al., 2017). Similar to Vrieling et al. (2017),Fig. 8shows very distinct and uniform pheno-logical retrievals for the agricultural area (see alsoFig. 1,Table 4). The natural vegetation communities reveal a greater heterogeneity, whereby pioneer-zone marshes generally green-up later than the higher marshes (see alsoWatkinson and Davy, 1985). The heterogeneity of phenology within natural classes may be partially explained by the large number of vegetation species with varying dominance with each vegetation community (Vrieling et al., 2017). The spatial patterns of the 2015 retrievals have a significant correlation with the 2016 Sentinel-2 retrievals with Pearson correlation coefficients of 0.47 for SOS20, 0.48 for SOS50, and 0.40 for PS90, and MSD of 20.4, 19.3, and 18.4, re-spectively. The positive MSD indicates that retrieved dates were ~19 days later for 2015, which can be explained by the relatively cold spring that year. As such, the future build-up of an archive of short-interval Sentinel-2 imagery may assist in understanding varying phe-nological responses to weather patterns for heterogeneous vegetation communities.

Our study demonstrated that sufficient Sentinel-2 observations could be obtained for Schiermonnikoog to model the seasonal devel-opment of vegetation greenness throughout 2016; only a few gaps of more than one month occurred. Although we could benefit from overlapping orbits, the launch of Sentinel-2B in March 2017 resulted in the same repeat frequency globally, while increasing the frequency further towards the poles. Cloud cover may alter image availability depending on location and season. As a consequence, single-year phe-nology retrievals may not be successful always and everywhere, even if our study site experiences a medium-to-high cloud cover relative to other locations across the globe (King et al., 2013). We showed that a reduced image availability during a specific month can have varying effects depending on how that month is placed in the vegetation cycle, and that the effects may be small when sufficient other data points remain to describe the vegetation growth. Further enhancing image availability may require incorporating Landsat imagery (Wang et al., 2017), utilize dedicated missions for short-repeat observation of small targeted areas like Venμs (Dedieu et al., 2006), or using commercial constellations of many small satellites (Strauss, 2017). Nonetheless, we expect accurate single-year retrieval of phenology atfine resolution to become feasible for large parts of the globe using Sentinel-2 only. 6. Conclusions

We showed that phenological parameters can now be retrieved at 10 m resolution from Sentinel-2 image time series, providing new op-portunities for landscape-scale phenology research. For areas with fine-scale spatial heterogeneity of vegetation communities, like Schiermonnikoog, Sentinel-2 allows retrieval and analysis of vegetation phenology with higher accuracy as compared to the commonly-used medium-resolution sensors like MODIS. Comparison of phenological parameters from satellite-derived NDVI series with camera-derived GCC series revealed that, despite their correlation, large differences existed between the two. We demonstrated that part of these differences could be explained by the vegetation index used and the different viewing angles for satellite andfield camera. In general, spring green-up was closer spaced in time between the two sources as compared to timing of senescence. We argued that this may be partially affected by factors that are difficult to model like flowering and presence of non-photo-synthetic elements like grass seed heads, which are more dominant in the oblique camera view as compared to the closer to nadir satellite view. Two-layer radiative transfer modelling including a top layer with reduced chlorophyll concentrations could explain part of the observed differences due to viewing angle. Future investigation in other vege-tation systems monitored by phenological cameras will allow to further

appreciate Sentinel-2 and similarﬁne-resolution satellites as key tools for studying vegetation phenology.

Acknowledgements

This work was partially ﬁnanced through the European Space Agency's Innovators-III project “Remote Sensing for Essential Biodiversity Variables” (ESRIN Contract No 4000113386). We thank Willem Nieuwenhuis (University of Twente) for his assistance with the atmospheric correction of the imagery. We thank Henk Kloosterman (University of Twente) for sharing his knowledge on the ﬂora of Schiermonnikoog. We much appreciate the constructive criticism of the three anonymous reviewers and the guest editors that helped to im-prove this work.

References

Alberton, B., Almeida, J., Helm, R., Torres, R.S., Menzel, A., Morellato, L.P.C., 2014. Using phenological cameras to track the green up in a cerrado savanna and its on-the-ground validation. Eco. Inform. 19, 62–70.

Anderson, H., Nilsen, L., Tømmervik, H., Karlsen, S., Nagai, S., Cooper, E., 2016. Using ordinary digital cameras in place of near-infrared sensors to derive vegetation indices for phenology studies of High Arctic vegetation. Remote Sens. 8, 847.

Beck, P.S.A., Atzberger, C., Høgda, K.A., Johansen, B., Skidmore, A.K., 2006. Improved monitoring of vegetation dynamics at very high latitudes: a new method using MODIS NDVI. Remote Sens. Environ. 100, 321–334.

de Beurs, K.M., Henebry, G.M., 2005. Land surface phenology and temperature variation in the International Geosphere-Biosphere Program high-latitude transects. Glob. Chang. Biol. 11, 779–790.

Bockelmann, A.-C., Bakker, J.P., Neuhaus, R., Lage, J., 2002. The relation between ve-getation zonation, elevation and inundation frequency in a Wadden Sea salt marsh. Aquat. Bot. 73, 211–221.

Chen, J., Jönsson, P., Tamura, M., Gu, Z.H., Matsushita, B., Eklundh, L., 2004. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golayﬁlter. Remote Sens. Environ. 91, 332–344.

Claverie, M., Demarez, V., Duchemin, B., Hagolle, O., Ducrot, D., Marais-Sicre, C., Dejoux, J.F., Huc, M., Keravec, P., Béziat, P., Fieuzal, R., Ceschia, E., Dedieu, G., 2012. Maize and sunﬂower biomass estimation in southwest France using high spatial and tem-poral resolution remote sensing data. Remote Sens. Environ. 124, 844–857.

Cleland, E.E., Chuine, I., Menzel, A., Mooney, H.A., Schwartz, M.D., 2007. Shifting plant phenology in response to global change. Trends Ecol. Evol. 22, 357–365.

Dedieu, G., Karnieli, A., Hagolle, O., Jeanjean, H., Cabot, F., Ferrier, P., Yaniv, Y., 2006. Venμs: a joint Israel–French Earth Observation scientiﬁc mission with high spatial and temporal resolution capabilities. In: Sobrino, J.A. (Ed.), Recent Advances in Quantitative Remote Sensing. University of Valencia, Valencia, Spain, pp. 517–521.

Deering, D.W., Middleton, E.M., Irons, J.R., Blad, B.L., Walter-Shea, E.A., Hays, C.L., Walthall, C., Eck, T.F., Ahmad, S.P., Banerjee, B.P., 1992. Prairie grassland bidirec-tional reﬂectances measured by diﬀerent instruments at the FIFE site. J. Geophys. Res.-Atmos. 97, 18,887–818,903.

Drusch, M., Del Bello, U., Carlier, S., Colin, O., Fernandez, V., Gascon, F., Hoersch, B., Isola, C., Laberinti, P., Martimort, P., Meygret, A., Spoto, F., Sy, O., Marchese, F., Bargellini, P., 2012. Sentinel-2: ESA's optical high-resolution mission for GMES op-erational services. Remote Sens. Environ. 120, 25–36.

Durant, J.M., Hjermann, D.Ø., Ottersen, G., Stenseth, N.C., 2007. Climate and the match or mismatch between predator requirements and resource availability. Clim. Res. 33, 271–283.

Emelyanova, I.V., McVicar, T.R., Van Niel, T.G., Li, L.T., van Dijk, A.I.J.M., 2013. Assessing the accuracy of blending Landsat–MODIS surface reﬂectances in two landscapes with contrasting spatial and temporal dynamics: a framework for algo-rithm selection. Remote Sens. Environ. 133, 193–209.

European Space Agency, 2017. Sentinel-2 Spectral Response Functions (S2-SRF) issue 2.0 (COPE-GSEG-EOPG-TN-15-0007). In.

Féret, J.B., François, C., Asner, G.P., Gitelson, A.A., Martin, R.E., Bidel, L.P.R., Ustin, S.L., le Maire, G., Jacquemoud, S., 2008. PROSPECT-4 and 5: advances in the leaf optical properties model separating photosynthetic pigments. Remote Sens. Environ. 112, 3030–3043.

Fisher, J.I., Mustard, J.F., Vadeboncoeur, M.A., 2006. Green leaf phenology at Landsat resolution: scaling from theﬁeld to the satellite. Remote Sens. Environ. 100, 265–279.

Frantz, D., Stellmes, M., Röder, A., Udelhoven, T., Mader, S., Hill, J., 2016. Improving the spatial resolution of land surface phenology by fusing medium- and coarse-resolution inputs. IEEE Trans. Geosci. Remote Sens. 54, 4153–4164.

Gao, F., Hilker, T., Zhu, X., Anderson, M., Masek, J., Wang, P., Yang, Y., 2015. Fusing Landsat and MODIS data for vegetation monitoring. In: IEEE Geoscience and Remote Sensing Magazine. vol. 3. pp. 47–60.

Gao, F., Anderson, M.C., Zhang, X., Yang, Z., Alﬁeri, J.G., Kustas, W.P., Mueller, R., Johnson, D.M., Prueger, J.H., 2017. Toward mapping crop progress atﬁeld scales through fusion of Landsat and MODIS imagery. Remote Sens. Environ. 188, 9–25.

Gillespie, A.R., Kahle, A.B., Walker, R.E., 1987. Color enhancement of highly correlated images. II. Channel ratio and“chromaticity” transformation techniques. Remote Sens.

(13)

Environ. 22, 343–365.

Jacquemoud, S., Baret, F., 1990. PROSPECT: a model of leaf optical properties spectra. Remote Sens. Environ. 34, 75–91.

Jentsch, A., Kreyling, J., Boettcher-Treschkow, J., Beierkuhnlein, C., 2009. Beyond gra-dual warming: extreme weather events alterﬂower phenology of European grassland and heath species. Glob. Chang. Biol. 15, 837–849.

King, M.D., Platnick, S., Menzel, W.P., Ackerman, S.A., Hubanks, P.A., 2013. Spatial and temporal distribution of clouds observed by MODIS onboard the Terra and Aqua satellites. IEEE Trans. Geosci. Remote Sens. 51, 3826–3852.

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

Kuusk, A., 1991. The hot spot eﬀect in plant canopy reﬂectance. In: Myneni, R.B., Ross, J. (Eds.), Photon-Vegetation Interactions: Applications in Optical Remote Sensing and Plant Ecology. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 139–159.

Liao, C.Z., Luo, Y.Q., Fang, C.M., Chen, J.K., Li, B., 2008. Litter pool sizes, decomposition, and nitrogen dynamics in Spartina alterniﬂora-invaded and native coastal marsh-lands of the Yangtze Estuary. Oecologia 156, 589–600.

Louis, J., Debaecker, V., Pﬂug, B., Main-Knorn, M., Bieniarz, J., Mueller-Wilm, U., Cadau, E., Gascon, F., 2016. Sentinel-2 Sen2Cor: L2A Processor for Users. European Space Agency, (Special Publication) ESA SP.

Markwardt, C.B., 2008. Non-linear least squaresﬁtting in IDL with MPFIT. In: Lewis, J., Argyle, R., Bunclarck, P., Evans, D., Gonzales-Solares, E. (Eds.), ASP Conference Series. Vol. XXX. Quebec.

Melaas, E.K., Friedl, M.A., Zhu, Z., 2013. Detecting interannual variation in deciduous broadleaf forest phenology using Landsat TM/ETM+ data. Remote Sens. Environ. 132, 176–185.

Melaas, E.K., Sulla-Menashe, D., Gray, J.M., Black, T.A., Morin, T.H., Richardson, A.D., Friedl, M.A., 2016. Multisite analysis of land surface phenology in North American temperate and boreal deciduous forests from Landsat. Remote Sens. Environ. 186, 452–464.

Menzel, A., Helm, R., Zang, C., 2015. Patterns of late spring frost leaf damage and re-covery in a European beech (Fagus sylvatica L.) stand in south-eastern Germany based on repeated digital photographs. Front. Plant Sci. 6.

Meroni, M., Verstraete, M.M., Rembold, F., Urbano, F., Kayitakire, F., 2014. A phenology-based method to derive biomass production anomalies for food security monitoring in the Horn of Africa. Int. J. Remote Sens. 35, 2472–2492.

Migliavacca, M., Galvagno, M., Cremonese, E., Rossini, M., Meroni, M., Sonnentag, O., Cogliati, S., Manca, G., Diotri, F., Busetto, L., Cescatti, A., Colombo, R., Fava, F., Morra di Cella, U., Pari, E., Siniscalco, C., Richardson, A.D., 2011. Using digital re-peat photography and eddy covariance data to model grassland phenology and photosynthetic CO2 uptake. Agric. For. Meteorol. 151, 1325–1337.

Miles, W.T.S., Bolton, M., Davis, P., Dennis, R., Broad, R., Robertson, I., Riddiford, N.J., Harvey, P.V., Riddington, R., Shaw, D.N., Parnaby, D., Reid, J.M., 2017. Quantifying full phenological event distributions reveals simultaneous advances, temporal stabi-lity and delays in spring and autumn migration timing in long-distance migratory birds. Glob. Chang. Biol. 23, 1400–1414.

Mizunuma, T., Wilkinson, M.L., Eaton, E., Mencuccini, M.I.L., Morison, J., Grace, J., 2013. The relationship between carbon dioxide uptake and canopy colour from two camera systems in a deciduous forest in southern England. Funct. Ecol. 27, 196–207.

Moore, K.J., Moser, L.E., Vogel, K.P., Waller, S.S., Johnson, B.E., Pedersen, J.F., 1991. Describing and quantifying growth stages of perennial forage grasses. Agron. J. 83, 1073–1077.

Moore, C.E., Brown, T., Keenan, T.F., Duursma, R.A., Van Dijk, A.I.J.M., Beringer, J., Culvenor, D., Evans, B., Huete, A., Hutley, L.B., Maier, S., Restrepo-Coupe, N., Sonnentag, O., Specht, A., Taylor, J.R., Van Gorsel, E., Liddell, M.J., 2016. Reviews and syntheses: Australian vegetation phenology: new insights from satellite remote sensing and digital repeat photography. Biogeosciences 13, 5085–5102.

Nasahara, K.N., Nagai, S., 2015. Review: development of an in situ observation network for terrestrial ecological remote sensing: the Phenological Eyes Network (PEN). Ecol. Res. 30, 211–223.

Newman, G., Wiggins, A., Crall, A., Graham, E., Newman, S., Crowston, K., 2012. The future of citizen science: emerging technologies and shifting paradigms. Front. Ecol. Environ. 10, 298–304.

Nijland, W., de Jong, R., de Jong, S.M., Wulder, M.A., Bater, C.W., Coops, N.C., 2014. Monitoring plant condition and phenology using infrared sensitive consumer grade digital cameras. Agric. For. Meteorol. 184, 98–106.

Nijland, W., Bolton, D.K., Coops, N.C., Stenhouse, G., 2016. Imaging phenology; scaling from camera plots to landscapes. Remote Sens. Environ. 177, 13–20.

O'Connor, B., Skidmore, A., Darvishzadeh, R., Wang, T., McOwen, C., Vrieling, A., Harfoot, M., 2017. Remote sensing for essential biodiversity variables (RS4EBV):ﬁnal synthesis report. In: ESRIN Contract No 4000113386, (p. 58).

Olff, H., Bakker, J.P., Fresco, L.F.M., 1988. The effect of fluctuations in tidal inundation frequency on a salt-marsh vegetation. Vegetatio 78, 13–19.

Pan, Z., Huang, J., Zhou, Q., Wang, L., Cheng, Y., Zhang, H., Blackburn, G.A., Yan, J., Liu, J., 2015. Mapping crop phenology using NDVI time-series derived from HJ-1 A/B data. Int. J. Appl. Earth Obs. Geoinf. 34, 188–197.

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

Pranger, D.P., Tolman, M.E., 2012. Toelichting bij de Vegetatiekartering

Schiermonnikoog 2010 op basis van false colour-luchtfoto's 1: 10.000 [Explanation to the Vegetation Mapping Schiermonnikoog 2010 on the Basis of False Colour Aerial Photographs 1:10.000, in Dutch] (p. 128). Rijkswaterstaat (RWS-DID), Delft, The Netherlands.

Reed, B.C., Brown, J.F., Vanderzee, D., Loveland, T.R., Merchant, J.W., Ohlen, D.O., 1994. Measuring phenological variability from satellite imagery. J. Veg. Sci. 5, 703–714.

Richardson, A.D., Braswell, B.H., Hollinger, D.Y., Jenkins, J.P., Ollinger, S.V., 2009. Near-surface remote sensing of spatial and temporal variation in canopy phenology. Ecol. Appl. 19, 1417–1428.

Schnelle, F., Volkert, E., 1964. Internationale phänologische gärten Stationen eines grundnetzes für internationale phänologische beobachtungen. Agric. Meteorol. 1, 22–29.

Schwartz, M.D., Betancourt, J.L., Weltzin, J.F., 2012. From caprio's lilacs to the USA National Phenology Network. Front. Ecol. Environ. 10, 324–327.

Shariati Najafabadi, M., Darvishzadeh, R., Skidmore, A.K., Kölzsch, A., Vrieling, A., Nolet, B.A., Exo, K.-M., Meratnia, N., Havinga, P.J.M., Stahl, J., Toxopeus, A.G., 2015. Satellite- versus temperature-derived green wave indices for predicting the timing of spring migration of avian herbivores. Ecol. Indic. 58, 322–331.

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

Stöckli, R., Rutishauser, T., Dragoni, D., O'Keefe, J., Thornton, P.E., Jolly, M., Lu, L., Denning, A.S., 2008. Remote sensing data assimilation for a prognostic phenology model. J. Geophys. Res. Biogeosci. 113.

Strauss, M., 2017. Planet Earth to get a daily selﬁe. Science 355, 782–783.

Ter Heerdt, G.N.J., Bakker, J.P., De Leeuw, J., 1991. Seasonal and spatial variation in living and dead plant material in a grazed grassland as related to plant species di-versity. J. Appl. Ecol. 28, 120–127.

Verhoef, W., 1984. Light scattering by leaf layers with application to canopy reﬂectance modeling: the SAIL model. Remote Sens. Environ. 16, 125–141.

Verhoef, W., Bach, H., 2003. Simulation of hyperspectral and directional radiance images using coupled biophysical and atmospheric radiative transfer models. Remote Sens. Environ. 87, 23–41.

van Vliet, A.J.H., Bron, W.A., Mulder, S., van der Slikke, W., Odé, B., 2014. Observed climate-induced changes in plant phenology in the Netherlands. Reg. Environ. Chang. 14, 997–1008.

Vrieling, A., de Beurs, K.M., Brown, M.E., 2011. Variability of African farming systems from phenological analysis of NDVI time series. Clim. Chang. 109, 455–477.

Vrieling, A., Skidmore, A.K., Wang, T., Meroni, M., Ens, B.J., Oosterbeek, K., O'Connor, B., Darvishzadeh, R., Heurich, M., Shepherd, A., Paganini, M., 2017. Spatially detailed retrievals of spring phenology from single-season high-resolution image time series. Int. J. Appl. Earth Obs. Geoinf. 59, 19–30.

Walker, J.J., de Beurs, K.M., Wynne, R.H., 2014. Dryland vegetation phenology across an elevation gradient in Arizona, USA, investigated with fused MODIS and Landsat data. Remote Sens. Environ. 144, 85–97.

Wang, Q., Blackburn, G.A., Onojeghuo, A.O., Dash, J., Zhou, L., Zhang, Y., Atkinson, P.M., 2017. Fusion of landsat 8 OLI and sentinel-2 MSI data. In: IEEE Transactions on Geoscience and Remote Sensing.

Watkinson, A.R., Davy, A.J., 1985. Population biology of salt marsh and sand dune an-nuals. Vegetatio 62, 487–497.

Wenhan, Q., 1993. Modeling bidirectional reﬂectance of multicomponent vegetation canopies. Remote Sens. Environ. 46, 235–245.

White, M.A., Thornton, P.E., Running, S.W., 1997. A continental phenology model for monitoring vegetation responses to interannual climatic variability. Glob. Biogeochem. Cycles 11, 217–234.

White, K., Pontius, J., Schaberg, P., 2014. Remote sensing of spring phenology in northeastern forests: a comparison of methods,ﬁeld metrics and sources of un-certainty. Remote Sens. Environ. 148, 97–107.

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

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

Zhang, X., Wang, J., Gao, F., Liu, Y., Schaaf, C., Friedl, M., Yu, Y., Jayavelu, S., Gray, J., Liu, L., Yan, D., Henebry, G.M., 2017. Exploration of scaling eﬀects on coarse re-solution land surface phenology. Remote Sens. Environ. 190, 318–330.

Zhao, J., Zhang, Y., Tan, Z., Song, Q., Liang, N., Yu, L., Zhao, J., 2012. Using digital cameras for comparative phenological monitoring in an evergreen broad-leaved forest and a seasonal rain forest. Eco. Inform. 10, 65–72.

Zhu, X., Chen, J., Gao, F., Chen, X., Masek, J.G., 2010. An enhanced spatial and temporal adaptive reﬂectance fusion model for complex heterogeneous regions. Remote Sens. Environ. 114, 2610–2623.

Zipoli, G., Grifont, D., 1994. Panicle contribution to bidirectional reﬂectance factors of a wheat canopy. Int. J. Remote Sens. 15, 3309–3314.