Improved Estimates of Arctic Land Surface Phenology Using Sentinel-2 Time Series

: The high spatial resolution and revisit time of Sentinel-2A / B tandem satellites allow a potentially improved retrieval of land surface phenology (LSP). The biome and regional characteristics, however, greatly constrain the design of the LSP algorithms. In the Arctic, such biome-speciﬁc characteristics include prolonged periods of snow cover, persistent cloud cover, and shortness of the growing season. Here, we evaluate the feasibility of Sentinel-2 for deriving high-resolution LSP maps of the Arctic. We extracted the timing of the start and end of season (SoS and EoS, respectively) for the years 2019 and 2020 with a simple implementation of the threshold method in Google Earth Engine (GEE). We found a high level of similarity between Sentinel-2 and PhenoCam metrics; the best results were observed with Sentinel-2 enhanced vegetation index (EVI) (root mean squared error (RMSE) and mean error (ME) of 3.0 d and − 0.3 d for the SoS, and 6.5 d and − 3.8 d for the EoS, respectively), although other vegetation indices presented similar performances. The phenological maps of Sentinel-2 EVI compared well with the same maps extracted from the Moderate Resolution Imaging Spectroradiometer (MODIS) in homogeneous landscapes (RMSE and ME of 9.2 d and 2.9 d for the SoS, and 6.4 and − 0.9 d for the EoS, respectively). Unreliable LSP estimates were ﬁltered and a quality ﬂag indicator was activated when the Sentinel-2 time series presented a long period ( > 40 d) of missing data; discontinuities were lower in spring and early summer (9.2%) than in late summer and autumn (39.4%). The Sentinel-2 high-resolution LSP maps and the GEE phenological extraction method will support vegetation monitoring and contribute to improving the representation of Artic vegetation phenology in land surface models.


Introduction
Studies of the Arctic are becoming increasingly important, particularly in the context of the onset of an expected tipping point in the function of its ecosystems as a result of ongoing climate warming [1] that is lengthening the growing season and increasing vegetation productivity [2,3]. Rising levels of vegetation productivity in the Arctic, characterized as greening, result in reduced albedo that further drives the warming trend [4,5] and has consequences for the regional carbon cycle [6]. For example, permafrost is thawing at an accelerating pace in a process that releases greenhouse gases and, consequently, exacerbates the warming trend in the region [7].
Land surface phenology (LSP), the study of remotely-sensed seasonal patterns in vegetation growth [8], complements the sparse field observations at high latitudes and is essential in the assessment and monitoring of responses of arctic vegetation to climate warming. Reliable LSP maps of the Arctic could be used to support models that show positive feedbacks between trends in climate warming and a lengthening growing season [3], while studies of Arctic LSP may reveal a link between the advancement of spring onset, drought severity, and an increase in the incidence of artic fires over recent decades [9], and support reports of links between carbon uptake and phenology metrics related to vegetation greenness, such as the amplitude of the vegetation index [10].
Studies of Arctic LSP have tended to use moderate-resolution satellite data, such as the Advanced Very-High-Resolution Radiometer (AVHRR) (1.1 km) and the Moderate Resolution Imaging Spectroradiometer (MODIS) (500 m) [2,3]. At these high latitudes, there are two key challenges for LSP analysis [10]. Firstly, LSP methods tend to erroneously detect the end of a snow period as the start of season (SoS), particularly in methods that extract phenology metrics from greenness indices, such as the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI), where distinct changes in reflectance during the snowmelt period are incorrectly detected as the onset of vegetation growth. To account for this problem, vegetation indices that are insensitive to snow have been developed to improve phenology estimation [11,12]. Secondly, a lack of valid satellite observations due to persistent cloud cover hampers LSP estimation, so biome-specific algorithms based on the combination of multi-satellite data and spatio-temporal gap filling methods have been developed [3].
Recent advances in remote sensing technologies present opportunities for the estimation of LSP at greater spatial resolution. Sentinel-2 mission provides decametric images with frequent revisit times (<5 d), allowing the extraction of phenology metrics [13] and study of vegetation dynamics at the canopy scale, while the development of cloud-based platforms, such as Google Earth Engine (GEE) [14], allows processing of large volumes of satellite data for planetary-scale analysis [15,16] and increases the accessibility of high-resolution satellite archive data required for time series analyses.
The objectives of this study were (1) evaluating the feasibility of Sentinel-2 for LSP retrieval in the Arctic at a spatial resolution of 10 m, (2) proposing a novel and fast cloud computing implementation in GEE of the widely used threshold phenological extraction method, and (3) assessing the performances of the Arctic Sentinel-2 LSP maps for the SoS and end of season (EoS) for the years 2019 and 2020 based on the comparison with MODIS LSP and PhenoCam ground data.

Study Area
We generated SoS and EoS maps over regions classified as tundra in the RESOLVE Ecoregions dataset 2017 [17]. High latitudes of the Arctic, where winters are cold and summers are short, are mostly uninhabited by human populations; here, tundra soils contain a layer of permafrost that prevents the growth of trees, but supports the growth of grass and shrub vegetation during the short summer period between June and August, as revealed by vegetation index time series (e.g., NDVI and EVI).

Data
We generated LSP metrics for 2019 and 2020 using Sentinel-2 level-2A data, which provide daily top-of-canopy reflectance at 10, 20, and 60 m of spatial resolution; from these data, we used the 10-m resolution bands 2, 3, 4, and 8, and the 20-m band 12. Sentinel-2A and -2B multispectral satellites were launched in 2015 and 2017, respectively, and have a revisit time of 5 d at the equator that decreases with increasing latitude. The maximum revisit time between 1 May and 30 September 2019 for latitudes between 70 and 75 • was 1.7 d (average: 0.9 d) (Figure 1a). The LSP metrics estimated with Sentinel-2 were compared with the same metrics estimated from the 500 m MOD09GAv6 product [18].
Maximum discontinuity in the time series from spring and early summer (1 May to 15 July) and later summer and autumn (15 July to 30 September) was plotted after non-valid observations had been filtered using values of the quality band scene classification layer (SCL) provided in the Sentinel-2 level-2A, where 1 = saturated or defective; 2 = dark area pixels; 3 = cloud shadows; 6 = water; and, 7−10 = clouds and cirrus (Figure 1b,c).

Phenology Extraction
We adopted a widely used threshold method [13,19,20] that assigns the SoS and EoS as the first and last days of the season, respectively, on which a threshold u is exceeded; u may be a constant or defined dynamically for each pixel [21]. In this study, we estimated u as a dynamic value that depends on the annual amplitude of the time series (Equation (1)): where Vmin and Vmax are the minimum and maximum annual values in the time series, respectively, and p is a given proportion (%) of the amplitude. In this study, we used p = 0.5 as the mid-greenup and mid-greendown of the growing season. The threshold metrics estimated with 50% of the amplitude are less affected by biases due to discontinuities in the time series [13]. We applied two variants of the threshold method as follows: 1. Threshold method without smoothing. The threshold method was applied directly to the daily time series. For the SoS, we searched for the earliest date when the daily vegetation index exceeded u; then, we applied a linear interpolation between this first observation when the vegetation index was >u and the preceding observation from which to estimate SoS as the value >u ( Figure 2c). For the EoS, the linear interpolation was applied between the latest date when the vegetation index was >u, and the subsequent observation, where EoS corresponded to the linearly interpolated value >u. 2. Threshold method after smoothing. The time series data were smoothed prior to the extraction of LSP metrics (Figure 2d), as is common practice in LSP estimation to reduce noise and discontinuities of time series data [21]. The criteria for selection of the processing steps were based on the feasibility of their implementation in GEE, without comprising the recreation of the phenology curve (see GEE code in Supplementary Materials). Excessive smoothing of time series may lead to unrealistic recreations of the growing season. We first applied a moving average window, with an average radius of 10 d, every 20 d (Figure 2d); if a pixel in the 20 d composite window was empty due to a lack of valid observations, the window size was increased to 40 d. Next, we applied a cubic interpolation to convert the 20 d composites to a daily time series. The threshold was estimated from the amplitude of the interpolated time series, rather than with daily observations, and then the SoS and EoS were estimated as the first and last days, respectively, that exceeded the dynamic threshold in the interpolated time series.

Phenology Extraction
We adopted a widely used threshold method [13,19,20] that assigns the SoS and EoS as the first and last days of the season, respectively, on which a threshold u is exceeded; u may be a constant or defined dynamically for each pixel [21]. In this study, we estimated u as a dynamic value that depends on the annual amplitude of the time series (Equation (1)): where V min and V max are the minimum and maximum annual values in the time series, respectively, and p is a given proportion (%) of the amplitude. In this study, we used p = 0.5 as the mid-greenup and mid-greendown of the growing season. The threshold metrics estimated with 50% of the amplitude are less affected by biases due to discontinuities in the time series [13]. We applied two variants of the threshold method as follows: 1.
Threshold method without smoothing. The threshold method was applied directly to the daily time series. For the SoS, we searched for the earliest date when the daily vegetation index exceeded u; then, we applied a linear interpolation between this first observation when the vegetation index was >u and the preceding observation from which to estimate SoS as the value >u ( Figure 2c). For the EoS, the linear interpolation was applied between the latest date when the vegetation index was >u, and the subsequent observation, where EoS corresponded to the linearly interpolated value >u.

2.
Threshold method after smoothing. The time series data were smoothed prior to the extraction of LSP metrics (Figure 2d), as is common practice in LSP estimation to reduce noise and discontinuities of time series data [21]. The criteria for selection of the processing steps were based on the feasibility of their implementation in GEE, without comprising the recreation of the phenology curve (see GEE code in Supplementary Materials). Excessive smoothing of time series may lead to unrealistic recreations of the growing season. We first applied a moving average window, with an average radius of 10 d, every 20 d (Figure 2d); if a pixel in the 20 d composite window was empty due to a lack of valid observations, the window size was increased to 40 d. Next, we applied a cubic interpolation to convert the 20 d composites to a daily time series. The threshold was estimated from the amplitude of the interpolated time series, rather than with daily observations, and then the SoS and EoS were estimated as the first and last days, respectively, that exceeded the dynamic threshold in the interpolated time series.

Reclassification of Snow Observations in Green Chromatic Coordinate (GCC), Normalized Difference Vegetation Index (NDVI), and Enhanced Vegetation Index (EVI)
NDVI and EVI show a sharp increase after snowmelt in tundra vegetation (Figure 2a,b), because the previously snow-covered vegetation canopy is predominantly evergreen. The break in the time series attributed to the transition snow-to-vegetation and vegetation-to-snow is problematic in the estimation of LSP metrics. The SoS and EoS can be erroneously assigned to the snow transition dates instead to the actual vegetation dynamics. The reclassification of snow and post-thaw values is a common practice in Arctic LSP studies [13,24] to ensure a consistent threshold value (Equation (1)). We reclassified the snow observations to a fixed minimum value, which is specific for each vegetation index: GCC min = 0.31, NDVI min = 0.39, EVI min = 0.2 (See histograms in Supplementary Figure S2), and NDPI min = 0.24 (Supplementary Figure S1). This parameter was estimated as the mean value of the first Sentinel-2 snow-free observation of the year, extracted from 400 points randomly distributed. For the NDPI, the NDPI min corresponded to the mean NDPI value of snow observations for the optimal alpha (Supplementary Figure S1). The pixels that presented a maximum vegetation index in the time series below the minimum snow value were masked as non-vegetated pixels and the phenology metrics were not computed.

Implementation of Land Surface Phenology Algorithms in Google Earth Engine
The two variants of the threshold method were vectorized for their correct implementation in GEE. The vectorization is the process of transforming a code so that all components of an array are processed simultaneously [25]. This concept is contrary to the commonly used practice in LSP estimation, in which time series are processed separately, pixel by pixel, in a for loop. The use of for loops are, however, highly discouraged in GEE in preference for the recommended map functions. For instance, the moving average for the 20-day composition was implemented as a function that maps a list of dates (see code in Supplementary Materials). The function takes the dates as the input arguments to filter the Sentinel-2 collection with a window size of 20 days and, then, makes the average out of the selected images. The purpose of the map function is that each element of the array, in this case the dates, is processed separately and, consequently, each 20-day composition is generated at the same time.

Validation with PhenoCam
The SoS and EoS metrics extracted from Sentinel-2 were compared with the same metrics estimated from the near-surface reflectances of the PhenoCam network [26]. PhenoCam provides half-hourly images captured from digital cameras for 393 sites across North America and Europe, from which 15 cameras cover the tundra biome. At the moment of the analysis, GEE provided Sentinel-2 data from 2019 onwards, and only three sites presented available data in the Arctic for the

Implementation of Land Surface Phenology Algorithms in Google Earth Engine
The two variants of the threshold method were vectorized for their correct implementation in GEE. The vectorization is the process of transforming a code so that all components of an array are processed simultaneously [25]. This concept is contrary to the commonly used practice in LSP estimation, in which time series are processed separately, pixel by pixel, in a for loop. The use of for loops are, however, highly discouraged in GEE in preference for the recommended map functions. For instance, the moving average for the 20-day composition was implemented as a function that maps a list of dates (see code in Supplementary Materials). The function takes the dates as the input arguments to filter the Sentinel-2 collection with a window size of 20 days and, then, makes the average out of the selected images. The purpose of the map function is that each element of the array, in this case the dates, is processed separately and, consequently, each 20-day composition is generated at the same time.

Validation with PhenoCam
The SoS and EoS metrics extracted from Sentinel-2 were compared with the same metrics estimated from the near-surface reflectances of the PhenoCam network [26]. PhenoCam provides half-hourly images captured from digital cameras for 393 sites across North America and Europe, from which 15 cameras cover the tundra biome. At the moment of the analysis, GEE provided Sentinel-2 data from 2019 onwards, and only three sites presented available data in the Arctic for the years 2019 and 2020. The PhenoCam data used in the study was presented as provisional near-real time and subject to changes. The coordinates of these three PhenoCam cameras, used in the study, are shown in Supplementary Table S1. PhenoCam also provides daily time series of GCC derived from different regions of interest observed by the digital camera [26]. Such regions of interest cover the vegetation types in the site. For the selected sites, however, the camera only observed one vegetation type consisting of tundra grasses.
The GCC has been proved a good index for LSP estimation [22] since the time series does not present the sharp increase during the snowmelt. However, we also reclassified the snow values in PhenoCam GCC to the minimum value of the snow-free time series in order to be consistent with the forcing applied in the other Sentinel-2 vegetation indices. The snow observations in GCC were identified and selected by their similar values before and after the growing season (see Supplementary  Figures S3-S5).
We extracted the LSP metrics with a 50% threshold method from the daily GCC time series for the years 2019 and 2020 and compared them with the same LSP metrics estimated from the four Sentinel-2 vegetation indices (GCC, NDVI, EVI, and NDPI). The statistics that we reported were the mean error (ME) as the bias metric and the root mean squared error (RMSE) as the accuracy metric. Although the selected PhenoCam sites cover a homogeneous area, the location of the sites were manually relocated to ensure consistency with the area observed by the digital camera. The adjustment of the site coordinates was based on the orientation of the camera; the coordinates were relocated to 50 m from the original coordinates in the direction of the observation of the camera.

Comparison with the Moderate Resolution Imaging Spectroradiometer (MODIS) Land Surface Phenology
Our estimates of SoS and EoS using Sentinel-2 time series were compared with estimates based on MODIS time series data, using the same metrics and LSP extraction method; we also applied the reclassification of snow values to the MODIS time series. To compare both satellite datasets, we resized the 10 m Sentinel-2 LSP estimates to the 500 m MODIS projection. The Sentinel-2 was aggregated by averaging the pixels that lie within a 500 m MODIS pixel. The resampled SoS and EoS estimated with Sentinel-2 and the original MODIS LSP metrics, both at 500 m, were compared pixel-wise and the ME and RMSE were reported. To ensure reliable estimates in MODIS, we only considered the MODIS LSP estimates obtained from time series that presented a gap lower than 10 days. This condition did not result in an excessive filtering of pixels since the MODIS time series had a maximum gap that generally did not exceed 20 days (Supplementary Figure S6).
To test the reliability of our MODIS LSP estimates, we compared the SoS and EoS with the 'Mid_greenup' and 'Mid_greendown' layers, respectively, of the MODIS Land Cover Dynamics product (MCD12Q2v6) [27]. The MCD12Q2v6 product was also estimated using the threshold method, with a dynamic threshold of 50% in the 'Mid_greenup' and 'Mid_greendown' layers [19]. MCD12Q2v6 data for 2019 were unavailable at the time of this analysis, so we compared our MODIS LSP estimates for 2018 with the MCD12Q2v6 layers for 2018.

Results
The Sentinel-2 vegetation index that showed the best results in the comparison with PhenoCam was the EVI (SoS and EoS ME: −0.3 and −3.8 d, and SoS and EoS RMSE: 3.0 and 6.5 d) for the threshold method without time series smoothing (Figure 3). Results with the threshold method with smoothing were less conclusive (Supplementary Figure S7); RMSE and ME results were uneven between SoS and EoS for the same LSP method and none of the vegetation indices excelled in both ME and RMSE. Overall, the four Sentinel-2 vegetation indices compared well with PhenoCam, and the RMSE and ME were generally below 10 d for the two threshold methods (see time series of PhenoCam and Sentinel-2 in Supplementary Figures S3-S5). for Sentinel-2. The phenology metrics were extracted with a 50% threshold method without time series smoothing. The bias between PhenoCam and Sentinel-2 is reported with the mean error (ME) and the accuracy with the root mean squared error (RMSE).
LSP maps generated using Sentinel-2 EVI (Figure 4a,b) showed a high level of similarity at the continental scale compared with the same phenology metrics estimated with MODIS time series (Supplementary Figure S8). However, at the local scale, differences were apparent particularly in the EoS; similarly, the map of length of season (LoS) (Figure 5a), which represented the difference between EoS and SoS, showed spatial patterns and latitudinal gradients that are expected across the region, with shorter LoS in the surrounding areas of the Kara Sea and Arctic Archipelago (Figure 5b), and in elevated areas, such as the Central Siberian Plateau (Figure 5c). At the local level, the SoS and EoS maps derived from Sentinel-2 showed a high degree of detail that moderate scales of resolution were unable to capture. For instance, Figure 5d shows the positive relationship between LoS and elevation in the Scandinavian mountains. The altitudinal gradient is observed at the canopy scale, which allows the analysis of LSP metrics by vegetation type. Figure 5e shows another example of the effect of topography on the LoS in the Ural mountains. In this second example, the LoS is different in both sites of the mountain range depending on the presence of glaciers. Supplementary Figure S9 shows the different spring growth onset depending on the land cover type, tundra shrublands and herbaceous cover in the delta of the Lena River.  LSP maps generated using Sentinel-2 EVI (Figure 4a,b) showed a high level of similarity at the continental scale compared with the same phenology metrics estimated with MODIS time series (Supplementary Figure S8). However, at the local scale, differences were apparent particularly in the EoS; similarly, the map of length of season (LoS) (Figure 5a), which represented the difference between EoS and SoS, showed spatial patterns and latitudinal gradients that are expected across the region, with shorter LoS in the surrounding areas of the Kara Sea and Arctic Archipelago (Figure 5b), and in elevated areas, such as the Central Siberian Plateau (Figure 5c). At the local level, the SoS and EoS maps derived from Sentinel-2 showed a high degree of detail that moderate scales of resolution were unable to capture. For instance, Figure 5d shows the positive relationship between LoS and elevation in the Scandinavian mountains. The altitudinal gradient is observed at the canopy scale, which allows the analysis of LSP metrics by vegetation type. Figure 5e shows another example of the effect of topography on the LoS in the Ural mountains. In this second example, the LoS is different in both sites of the mountain range depending on the presence of glaciers. Supplementary Figure S9 shows the different spring growth onset depending on the land cover type, tundra shrublands and herbaceous cover in the delta of the Lena River.  for Sentinel-2. The phenology metrics were extracted with a 50% threshold method without time series smoothing. The bias between PhenoCam and Sentinel-2 is reported with the mean error (ME) and the accuracy with the root mean squared error (RMSE).
LSP maps generated using Sentinel-2 EVI (Figure 4a,b) showed a high level of similarity at the continental scale compared with the same phenology metrics estimated with MODIS time series (Supplementary Figure S8). However, at the local scale, differences were apparent particularly in the EoS; similarly, the map of length of season (LoS) (Figure 5a), which represented the difference between EoS and SoS, showed spatial patterns and latitudinal gradients that are expected across the region, with shorter LoS in the surrounding areas of the Kara Sea and Arctic Archipelago (Figure 5b), and in elevated areas, such as the Central Siberian Plateau (Figure 5c). At the local level, the SoS and EoS maps derived from Sentinel-2 showed a high degree of detail that moderate scales of resolution were unable to capture. For instance, Figure 5d shows the positive relationship between LoS and elevation in the Scandinavian mountains. The altitudinal gradient is observed at the canopy scale, which allows the analysis of LSP metrics by vegetation type. Figure 5e shows another example of the effect of topography on the LoS in the Ural mountains. In this second example, the LoS is different in both sites of the mountain range depending on the presence of glaciers. Supplementary Figure S9 shows the different spring growth onset depending on the land cover type, tundra shrublands and herbaceous cover in the delta of the Lena River.   The comparison between the LSP metrics estimated using MODIS and Sentinel-2 resized to 500 m, in which MODIS pixels were filtered when variance of the Sentinel-2 phenology estimates within the 500 m pixel were >5 d 2 , showed a high level of similarity in homogeneous landscapes (Figure 6), regardless of data smoothing (Figure 6a,b), and greater similarity for SoS than EoS (e.g., threshold method with smoothing SoS and EoS RMSE: 8.1 and 10.4 d, respectively). The comparison of our MODIS estimates and the MCD12Q2v6 product (Supplementary Figure S10) showed a slight bias towards EoS (threshold method with smoothing ME: 4.9 and 6.4 d, respectively) and similarity with Sentinel-2 for SoS and EoS (threshold method without smoothing RMSE = 10.2 and 9.8 d for SoS and EoS, respectively). The comparison between the LSP metrics estimated using MODIS and Sentinel-2 resized to 500 m, in which MODIS pixels were filtered when variance of the Sentinel-2 phenology estimates within the 500 m pixel were >5 d 2 , showed a high level of similarity in homogeneous landscapes (Figure 6), regardless of data smoothing (Figure 6a,b), and greater similarity for SoS than EoS (e.g., threshold method with smoothing SoS and EoS RMSE: 8.1 and 10.4 d, respectively). The comparison of our MODIS estimates and the MCD12Q2v6 product (Supplementary Figure S10) showed a slight bias towards EoS (threshold method with smoothing ME: 4.9 and 6.4 d, respectively) and similarity with Sentinel-2 for SoS and EoS (threshold method without smoothing RMSE = 10.2 and 9.8 d for SoS and EoS, respectively).
The phenology maps of the Arctic presented some unreliable LSP estimates, particularly when there were continuous gaps in the time series, due to the presence of clouds. The proportion of vegetated land pixels in the study area with a gap length that exceeded 40 days was 9.2% for the spring and early summer period and 39.4% for the late summer and autumn period. These percentages refer only to the pixels that presented an EVI value higher than 0.2 and, thus, were considered as vegetated land in the study. The percentage of pixels flagged as non-vegetated areas (EVI values lower than 0.2 during the entire time series) was 15.4%. Figure 7 illustrates the flags associated with the SoS estimates. The regions with low revisit times (Figure 1a) were prone to discontinuities in the time series. The phenology maps of the Arctic presented some unreliable LSP estimates, particularly when there were continuous gaps in the time series, due to the presence of clouds. The proportion of vegetated land pixels in the study area with a gap length that exceeded 40 days was 9.2% for the spring and early summer period and 39.4% for the late summer and autumn period. These percentages refer only to the pixels that presented an EVI value higher than 0.2 and, thus, were considered as vegetated land in the study. The percentage of pixels flagged as non-vegetated areas (EVI values lower than 0.2 during the entire time series) was 15.4%. Figure 7 illustrates the flags associated with the SoS estimates. The regions with low revisit times (Figure 1a) were prone to discontinuities in the time series.   The phenology maps of the Arctic presented some unreliable LSP estimates, particularly when there were continuous gaps in the time series, due to the presence of clouds. The proportion of vegetated land pixels in the study area with a gap length that exceeded 40 days was 9.2% for the spring and early summer period and 39.4% for the late summer and autumn period. These percentages refer only to the pixels that presented an EVI value higher than 0.2 and, thus, were considered as vegetated land in the study. The percentage of pixels flagged as non-vegetated areas (EVI values lower than 0.2 during the entire time series) was 15.4%. Figure 7 illustrates the flags associated with the SoS estimates. The regions with low revisit times ( Figure 1a) were prone to discontinuities in the time series.

Discussion
The high agreement between the phenology metrics estimated with Sentinel-2 and MODIS corroborates the feasibility of LSP estimation owing to the low revisit time of the Sentinel-2 time series [13] in a region prone to cloud coverage. We found that only 9.2% of pixels showed a discontinuity of >40 d in the time series during the green-up period. In contrast, the cloud coverage was more persistent in late summer and autumn (39.4% of pixels with >40-d gap), and this high level of cloud occurrence may explain the lower similarity between EoS estimated using Sentinel-2 and MODIS. Such discontinuities in time series make the use of gap-filling techniques and robust smoothing techniques necessary for the accurate estimation of LSP metrics [13]. Moreover, when time series present continuous gaps, the combined use of Sentinel-2 and Landsat-8 is recommended as, overall, the combination of sensors provides improved LSP estimates [28].
The comparison with PhenoCam further proves the feasibility of LSP estimation with Sentinel-2 with a variety of vegetation indices. Results showed that the EVI performed the best. It was shown by [28] that EVI was more suited than NDVI in temperate deciduous forests and [13] used EVI2 (two-band EVI) for continental-scale analysis, both using a combination of Landsat-8 and Sentinel-2. Despite this, an extensive analysis including more in situ observations is required to further justify the use of EVI in Arctic phenology studies. Moreover, the correct reclassification of snow values seems more determinant than the variable selection in high-latitudes.
The procedure for flagging pixels with discontinuities in the time series may be used to identify good-quality pixels for further analysis in phenology studies; however, this flag only works when clouds are successfully filtered. Non-valid observations that were not filtered, mainly cloud-contaminated pixels and cloud shadows that were not reflected in the quality band, are included in the time series and may underestimate the long gaps in the Sentinel-2 time series. Furthermore, such contaminated values may change the growing season curve in the vegetation index time series and lead to erroneous estimates of the phenology.
We found high levels of similarity between the phenology metrics estimated using Sentinel-2 and MODIS and PhenoCam, even though the comparison reflects changes in vegetation greenness that do not necessarily correspond to vegetation phenophases or vegetation productivity [8]. In situ databases, such as the National Phenology Network (NPN) [29], the Pan European Phenology database (PEP725) [30], and the FLUXNET network [31] are essential to ground-truth remotely sensed phenological changes in vegetation. Near-surface LSP, such as the PhenoCam network, supports the validation of satellite LSP estimates, as observed in the current study. Phenophase observations provide information on the timing of relevant stages in the annual life cycle of vegetation, such as the leaf-out, which may differ from the LSP metrics estimated from a satellite [20]. Similarly, the onset and the end of the photosynthetic activity, estimated from time series of carbon fluxes in FLUXNET, do not match with the LSP metrics depending on the forest ecosystem [32]. Phenophase, carbon fluxes, and near-surface LSP measurements are scarce in the Arctic, with the distribution of digital cameras generally restricted to northern areas of Alaska, so data are lacking for adequate validation of high-latitude LSP maps.
The use of cloud platforms, such as GEE, has made possible the estimation of LSP using Sentinel-2 time series for extensive regions, such as the Arctic. In this study, the generation of the phenology maps for the Arctic used approximately 133,000 Sentinel-2 images, solely for the year 2019; this type of platform is a useful tool that allows the scientific community to inspect such dense, high-spatial resolution time series. Furthermore, GEE allows sharing the code and data, and researchers can easily reproduce the algorithm and inspect the time series locally (see GEE code in Supplementary Materials).
The resulting phenology maps should be taken with caution when data availability is scarce, which occurs particularly during the EoS (39.4% of pixels with >40 d gap). The combination of Sentinel-2 and Landsat-8 is recommended when continuous gaps are present in the time series, and more elaborated phenological retrieval methods including gap-filling and robust smoothing techniques should be used instead. Despite this, sensor harmonization and temporal smoothing may improve the LSP retrievals but at the expenses of adding complexity to the preprocessing, which impedes the implementation and fast computation in GEE. Our SoS and EoS maps may benefit future Arctic phenology studies that aim to analyze spatial variability in vegetation dynamics and require a high degree of detail at the canopy level; test effects of spatial scaling on temporal changes in phenology; and identify homogeneous landscapes in which phenology dynamics are similar.

Conclusions
Here, we present the first high-spatial resolution maps (10 m) of SoS and EoS that fully cover the Arctic for the most recent years (2019-2020). We prove the feasibility of phenology metric extraction in the region solely with Sentinel-2 and using basic implementations of the threshold method in GEE; the high revisit time at high latitudes allows dense time series to be obtained under cloud-free conditions in 1-3 days and only 9.2% of pixels showed a discontinuity of >40 d in the time series during the green-up period. We propose a set of forcing values for the reclassification of snow observations in Sentinel-2 for three common vegetation indices used in LSP estimation: GCC, NDVI, and EVI; and also provide a re-adjusted parameterization of the NDPI specifically for Sentinel-2. Future work may use the adjusted parameters and forcing values along with the implementation of the threshold method in GEE for large-scale estimation of phenology metrics. The retrieved Sentinel-2 high-resolution LSP maps and the proposed GEE phenological extraction method will support monitoring vegetation changes at high-spatial resolution and are expected to contribute to the representation of Artic vegetation phenology in land surface models.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/22/3738/s1, Supplementary Materials contain one table, ten figures, the section Estimation of the optimal alpha in NDPI, and the GitHub link to the Google Earth Engine code. Table S1: PhenoCam sites used in the study, Figure S1: Estimation of the optimal alpha for the normalized difference phenology index (NDPI), Figure S2: Histogram of Sentinel-2 vegetation indices after snowmelt, Figure S3: Time series in the PhenoCam site imcrkrtussock, Figure  S4