Arctic-Boreal Lake Phenology Shows a Relationship between Earlier Lake Ice-Out and Later Green-Up

Satellite remote sensing has transformed our understanding of Earth processes. One component of the Earth system where large uncertainties remain are Arctic and boreal freshwater lakes. With only short periods of open water due to annual ice cover, lake productivity in these regions is extremely sensitive to warming induced changes in ice cover. At the same time, productivity dynamics in these lakes vary enormously, even over short distances, making it difficult to understand these potential changes. A major impediment to an improved understanding of lake dynamics has been sparsely distributed field measurements, in large part due to the complexity and expense of conducting scientific research in remote northern latitudes. This project overcomes that hurdle by using a new set of ‘eyes in the sky’, the Planet Labs CubeSat fleet, to observe 35 lakes across 3 different arctic-boreal ecoregions in western North America. We extract time series of lake reflectance to identify ice-out and green-up across three years (2017–2019). We find that lakes with later ice-out have significantly faster green-ups. Our results also show ice-out varies latitudinally by 38 days from south to north, but only varies across years by ~9 days. In contrast, green-up varied between years by 22 days in addition to showing significant spatial variability. We compare PlanetScope to Sentinel-2 data and independently validate our ice-out estimates, finding an ice-out mean absolute difference (MAD) ~9 days. This study demonstrates the potential of using CubeSat imagery to monitor the timing and magnitude of ice-off and green-up at high spatiotemporal resolution.


Introduction
High northern latitudes are rich in lakes [1,2] and experiencing dramatic warming [3] with global consequences. Arctic and boreal lakes are experiencing changes in the timing of ice-out [4][5][6] and primary productivity [7]. A recent study has shown that ice-out globally has advanced 8 days over the past 100 years [8] and 13 days since 2000 in parts of the Arctic [9]. Air temperature has been shown to be the primary control on ice-out timing [10,11] with ice-out therefore varying along latitudinal climate gradients.
Ice-out timing and duration in arctic and boreal lakes fundamentally structures lake ecology, including the onset and timing of lake productivity [12]. Because loss of ice cover can accelerate temperature-driven changes in water transparency, heat budgets and mixing regimes [13], climate change may have strong effects on lake productivity. Paleolimnological records from the Siberian Arctic verified that warmer, ice-free summer conditions were linked to the intensified aquatic primary productivity [14] and earlier ice-out with recent warming has been linked to earlier phytoplankton blooms and an earlier season start [15]. Arctic and boreal lakes are key players in global carbon cycles, and changes in the onset and duration of the ice-free growing period could alter their net contributions to carbon fluxes [16]. However, other studies show that impacts of climate on lake productivity are community-specific [17], and the relative influence of catchment versus climate processes on lake productivity remain an active area of research [18].
Despite the major recent advances in our earth observation capabilities, changes in the timing of ice-out and green-up, and links between the two, remain largely unknown, especially for smaller lakes (<1 km 2 ) [15]. A greater understanding of the drivers of earlier ice-out and increased lake productivity would clearly allow us to better understand and therefore predict the impacts of climate change on these important ecosystems. One barrier towards achieving this goal is the lack of spatially distributed observations of lake phenology, which studies have attempted to solve by using satellites to track lake ice-out and productivity. However, most past attempts at mapping trends in lake dynamics have relied on coarse temporal or spatial resolution satellites [19,20]. This is unfortunate, because the majority of lakes are small (<1 km) [1,2] and phenological shifts can happen over short timeframes.
Detecting fine-scale phenological shifts has only recently been made possible with the launch of European Space Agency's Sentinel-2 mission [21] and Planet Lab's PlanetScope constellation [22], which boast 5 and 2 day return periods respectively. Planet Lab has successfully launched over 1200 CubeSats (10 × 10 × 30 cm 3 ) since 2013. These spaceborne sensors are collecting optical data, in four or five bands, over the entire earth's surface each day at a spatial resolution of 2-5 m. Sentinel-2 (S2) consists of two satellites launched in 2015 and 2017 with a combined return time of 5 days. S2 has successfully shown high radiometric quality over inland waters at 10-30 m spatial resolution [23,24]. PlanetScope imagery suffers from inconsistencies in radiometric calibration due to the large number of sensors [25] in comparison to the publicly-funded Landsat and Sentinel series. Fusing Sentinel-2 and Planet is a strategic method that has been used in rangeland mapping to capture the benefits of each sensor while compensating for their respective drawbacks [26].
The objective of this study is to characterize patterns in variability and trends in lake color, as observed by PlanetScope CubeSats, across diverse arctic-boreal lakes. Lake color is designated an Essential Climate Variable by the Global Climate Observing System (GCOS) because it provides key insights into lake ecological and physical processes [27]. Lake color has been used to map gross primary productivity in arctic-boreal lakes [24], to detect changes lake area [28], and to map bathymetry [25,29,30] and water quality [29,31,32].
Here, we use high resolution satellite imagery time series to observe the timing of iceout and peak green-up in three different lake districts in western North America. We test the hypothesis that patterns in lake-ice are spatially coherent across latitudinal and geomorphic gradients. We also test for synchrony between ice-out and green-up across a latitudinal gradient. We compare PlanetScope surface reflectance to Sentinel-2 to assess the quality of the PlanetScope time series and we analyze climate reanalysis data to evaluate links between air temperature and phenological metrics. This study builds on the burgeoning interest in the impact of longer growing seasons lake ecological dynamics [33,34].

Study Lakes
A total of 35 lakes were sampled for lake color as part of the NASA Arctic and Boreal Vulnerability Experiment (ABoVE) project (Table 1). Study lakes were selected to represent a broad gradient of ecological, physical and climatic characteristics. Study sites were selected from three lake districts, including the Yukon Flats (Alaska, USA), the Peace-Athabasca Delta (Alberta, Canada) and the Yellowknife/Daring Lake Area (Northwest Territories, Canada) ( Figure 1). Of these lakes, 22 had a surface area > 0.1 ha and were included in the HydroLAKES dataset. The median surface area was 0.73 km 2 and the median lake depth was 1.7 m. Sites are abbreviated as follows: Yukon Flats (YF), Peace-Athabasca Delta (PAD), Yellowknife (YK), and Daring Lake (DL).

Satellite Observations
Planet Labs, Inc. operates a fleet of >100 PlanetScope CubeSats that capture imagery over the entire Earth's surface each day (Planet Labs 2020). This constellation of small

Satellite Observations
Planet Labs, Inc. operates a fleet of >100 PlanetScope CubeSats that capture imagery over the entire Earth's surface each day (Planet Labs 2020). This constellation of small satellites is in sun synchronous orbit and possess four spectral bands blue (455-515 nm), green (500-590 nm), red (590-670 nm), and NIR (780-860 nm) with a spatial resolution of 3-5 m and 12-bit radiometric resolution. The first PlanetScope CubeSat was launched in 2013. For near daily coverage, we used the surface reflectance product, which is produced using a radiometric calibration based in the 6SV 2.1 code [35].
Lake color time series were extracted ( Figure 2) from PlanetScope Level-2 imagery using a cloud-native serverless workflow tool called SWEEP [36] and from Sentinel-2 using Google Earth Engine [37,38]. SWEEP is a workflow management system that uses scatter and gather workflow constructs to achieve concurrency across tasks. The lake phenology workflow required executing the same tasks across many lakes. Using SWEEP thus reduced overall computational time and made the workflow more efficient. Time series of lake color (in red, green, blue and near-infrared bands) was extracted for each of the study lakes from PlanetScope near-daily overpasses. First, the PlanetScope archives were filtered to site coordinates (Table 1). Scenes were clipped to a 3 × 3 pixel area surrounding each sampling coordinate and raster values in each band were extracted at that sampling site. This process was repeated for every PlanetScope scene acquired over each lake from January 2017 to December 2019.

Time Series Modeling
The resulting time series (1 January 2017-31 December 2019) were extracted i format for daily time series modeling in R. Local weighted regression (loess) was u smoothen the resulting time series. The loess method defines a window to perfor smoothing assuming that the trend is linear within it. The stats package in R that con loess function was used with a span parameter value of 30%. Span of 30% translates proportion of data points that were used for the linear weighted fit [40]. Final visualiz and statistical analysis were performed in Python. A complete list of packages is giv Table S1.

Phenology Metrics
We identified peak green up and ice-out from the smoothed time series usin visible and near-infrared bands. In the visible bands, ice has high reflectance, appe as white or bright in contrast to open water, which can appear black, blue, brown the case of highly productive lakes, even green ( Figure 3). In some cases when ice is without snow, it can appear dark in the visible bands but will still have high reflec in the near-infrared. The European Space Agency's Sentinel-2 observation system consists of two satellites (Sentinel-2A and -2B) in a sun-synchronous orbit with a return period of 5 days at the equator [21]. With 12 bit radiometric resolution, the 13 spectral bands have a spatial resolution ranging from 10-60 m. The four spectral bands and their central wavelengths that correlate to PlanetScope include the blue (458-523 nm), green (543-578 nm), red (650-680 nm) and near-infrared (785-899 nm) bands that each possess 10 m spatial resolution. Level 2A products are orthorectified, radiometrically corrected surface reflectance produced using the sen2cor atmospheric correction [39]. Clouds and cloud shadows were masked using Scene Classification, Cloud Probability and Snow Probability Map bands in the Level 2A product. While Sentinel-2A and -2B were launched in 2015 and 2017, Level2A surface reflectance data has only been made available over our study area starting in 2019. Therefore, we can only compare PlanetScope and Sentinel-2 for the 2019 year. For this analysis, the blue (482 nm), green (560), red (665) and near-infrared (834 nm) bands were selected and will be referred to as B1 (blue), B2 (green), B3 (red) and B4 (NIR) for consistency with PlanetScope hereafter.

Time Series Modeling
The resulting time series (1 January 2017-31 December 2019) were extracted in .csv format for daily time series modeling in R. Local weighted regression (loess) was used to smoothen the resulting time series. The loess method defines a window to perform the smoothing assuming that the trend is linear within it. The stats package in R that contains loess function was used with a span parameter value of 30%. Span of 30% translates to the proportion of data points that were used for the linear weighted fit [40]. Final visualization and statistical analysis were performed in Python. A complete list of packages is given in Table S1.

Phenology Metrics
We identified peak green up and ice-out from the smoothed time series using the visible and near-infrared bands. In the visible bands, ice has high reflectance, appearing as white or bright in contrast to open water, which can appear black, blue, brown or, in the case of highly productive lakes, even green ( Figure 3). In some cases when ice is clear without snow, it can appear dark in the visible bands but will still have high reflectance in the near-infrared.  We first calculated the normalized difference vegetation index (NDVI) using the near-infrared (Rnir) and red (Rred) reflectance (NDVI = (Rnir − Rred)/(Rnir + Rred)) from the smoothed time series. NDVI can range from −1 to 1. Photosynthetic pigments reflect Rnir and absorb Rred [41], so positive values are typically associated with vegetation [42][43][44]. Pure water absorbs strongly in the near-infrared [45] and therefore clear deep water bodies typically have a negative NDVI. In complex coastal and inland waters, however, NDVI can be positive due to influences from water column constituents or, in the case of shallow water bodies, from bottom substrates. Factors contributing NIR reflectance in water bodies include turbidity [32,38,46], phytoplankton [47][48][49], cyanobacteria blooms [50,51], and submerged or emergent vegetation [52][53][54]. Field measurements of lake optical properties at these study sites have shown that these lakes are shallow and optically complex [24]  We first calculated the normalized difference vegetation index (NDVI) using the nearinfrared (R nir ) and red (R red ) reflectance (NDVI = (R nir − R red )/(R nir + R red )) from the smoothed time series. NDVI can range from −1 to 1. Photosynthetic pigments reflect R nir and absorb R red [41], so positive values are typically associated with vegetation [42][43][44]. Pure water absorbs strongly in the near-infrared [45] and therefore clear deep water bodies typically have a negative NDVI. In complex coastal and inland waters, however, NDVI can be positive due to influences from water column constituents or, in the case of shallow water bodies, from bottom substrates. Factors contributing NIR reflectance in water bodies include turbidity [32,38,46], phytoplankton [47][48][49], cyanobacteria blooms [50,51], and submerged or emergent vegetation [52][53][54]. Field measurements of lake optical properties at these study sites have shown that these lakes are shallow and optically complex [24] and that high benthic reflectance can lead to a positive NDVI over shallow, clear lakes. NDVI over shallow lakes is negative during ice-cover and positive after ice-out [55]. Therefore, ice-out was detected using as the first day of positive NDVI following the negative NDVI of the ice-cover period (Figure 3).
We selected the green band to identify peak greenness as it has both historically [45] and more recently [24] been shown to correlate with lake ecological dynamics. Green reflectance has been used for decades in models of phytoplankton and benthic productivity [56][57][58]. Peak green-up was defined as the maximum peak in green reflectance during the ice-free period ( Figure 3).
For each lake and each year, ice-out and green-up windows were identified by visually inspecting each reflectance time series. Window start and stop dates were selected for each region based on manual inspection. The local minimum (ice-out) and maximum (green-up) were selected from these windows for each lake. For each year, windows of selection were adjusted to account for year-to-year variability in timing of ice-out. Then, ice-out and green-up dates were plotted on the time series (Figures S5-S7) and visually inspected to make ensure windows had been selected appropriately (i.e. didn't miss the trough (ice-out) or peak (green-up).

Climate Data
Air temperature and precipitation records were extracted from the ERA5 dataset ( Figure S1) and compared to lake ice-out and green-up dates using European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis data. ERA5 data provides seamless hourly estimates of climate variables from 1979 to the present by combining historical observations with data assimilation and modeling. ERA5 data was released in March 2019, is available from the Copernicus Climate Change Service Data Store and was accessed through Google Earth Engine. Studies evaluating the ERA5 dataset have shown it's superiority to its predecessor ERA-Interim and other similar climate datasets [59,60].

Validation Data
For validation, ice-out was also identified manually by visually inspecting PlanetScope satellite time series. Ice-out was manually identified for each lake and year as the period when ice cover had receded from the lake (e.g., Figure 3b-e). Manual inspection dates were then compared to the NDVI-based ice-out to better quantify uncertainties related to these methods. For each lake, time series were manually inspected for the first ice-free date. The first ice-free date was determined by visually inspecting each scene and flagging when the lake sampling point was free of visible ice, including broken up or deteriorated ice cover (Figure 3d). Peak greenness could not be independently validated at this point because high resolution time series of in situ primary productivity data does not exist for these sites

Satellite Imagery Data Quantity and Quality
Over 1037 total PlanetScope scenes were acquired over our study area from 2017-2019, with an average of about 18-20 scenes per week per lake. The most observations occurred in August (~40 scenes) and the fewest observations were in January (~11 scenes) ( Figure S2). The lack of observations in snow-covered months is likely due to early quality control issues separating clouds from snow that led to fewer available images. Average PlanetScope summer reflectance in the green band (June-August) was 0.33 ± 0.22. Strong correlations were observed between all bands ( Figure S3) based on r 2 values ranging from 0.94-0.99. High correlations between visible bands has also been documented for Landsat [61] and presents both challenges and opportunities for modeling by introducing additional covariance but also presenting an opportunity to interpolate new bands [62].
From 2017-2019, 2930 Sentinel-2 Level 1C scenes that were acquired over the study area and 1246 Level 2A scenes were acquired during 2019 before quality control. Average Sentinel-2 summer reflectance in the green band (July-August) was 0.06 ± 0.15. Strong correlations, as indicated by r 2 ranging 0.84-0.99, between the visible and near-infrared bands were also observed in Sentinel-2 ( Figure S4). We note that the magnitude of surface reflectance values differs between PlanetScope and Sentinel-2, where PlanetScope has significantly higher reflectance (Figure 4). It is expected that reflectance values over open water should be very low (~10% or less) [63]. The lower reflectance of Sentinel-2 is due to higher radiometric fidelity, narrower bandpasses and a more restrictive cloud filter. The Sentinel-2 reflectance is closer to the surface reflectance values (0.01-0.08) measured in situ in the Yukon Flats using an ASD Handheld spectrometer [24].  The 2019 average mean absolute difference between PlanetScope and Sentinel-2 was 0.16 (unitless, surface reflectance) or 57% across all bands. The largest mean differences were observed in the NIR band (60%) and the smallest difference was observed in the green band (56%). We also calculated the difference on a weekly basis. We found differences were greater during the ice-free period than during the ice-covered period (Table S2, Figure S8). Seasonally, the differences were larger in summer and fall months than in spring and winter ( Figure 4 and Table S2).

Phenology Metrics
Ice-out dates ranged from as early as April 13th (Rat Lake, 2017) to as late as June 29th (Daring Lake, 2019) ( Figure S5-S7). Averaging sites, the mean annual ice-out date varied interannually only by 8 days (Figure 5a), falling on average on May 17th. As expected, ice-out occurs earlier at lower latitude sites (Figure 5c). Statistically significant differences between regions were found for the ice-out day-of-year (Kruskal Wallace H- The difference between PlanetScope and Sentinel were more pronounced in the summertime, where the MAPD between the sensors became as high as 55% (Figure 4a). Despite their differences in magnitude, a consistent relationship between 2019 weekly paired observations was observed for the two sensors (Figure 4b).
The 2019 average mean absolute difference between PlanetScope and Sentinel-2 was 0.16 (unitless, surface reflectance) or 57% across all bands. The largest mean differences were observed in the NIR band (60%) and the smallest difference was observed in the green band (56%). We also calculated the difference on a weekly basis. We found differences were greater during the ice-free period than during the ice-covered period (Table S2, Figure S8). Seasonally, the differences were larger in summer and fall months than in spring and winter ( Figure 4 and Table S2).

Phenology Metrics
Ice-out dates ranged from as early as April 13th (Rat Lake, 2017) to as late as 29 June (Daring Lake, 2019) ( Figures S5-S7). Averaging sites, the mean annual ice-out date varied interannually only by 8 days (Figure 5a), falling on average on 17 May. As expected, ice-out occurs earlier at lower latitude sites (Figure 5c). Statistically significant differences between regions were found for the ice-out day-of-year (Kruskal Wallace H-test, H = 47.92, p < 0.05, n = 105). Average ice-out at the lower latitude Peace-Athabasca sites fell on May 10th, more than a month earlier than the higher latitude Daring Lakes sites, which fell on 17 June on average (Figure 5c). The NDVI ice-out detection method showed a mean absolute difference of 9 days from the manual validation (n = 105). The NDVI method had a negative bias, meaning it detected an earlier ice-out than the manual method, and varied by region ( Figure S8). While this difference is large, it is still an improvement over the uncertainties introduced by the 16-day return period of Landsat. hydrologic connectivity, such as those in the PAD [64], may experience intense pulses of sediment in the spring during ice-out. Sediment has high reflectance in the near-infrared and may confound ice-detection [65]. Sediment could explain the negative model bias because it also reflects in the visible bands, thus making water appear "bright" like ice to the human eye even though ice may be gone. Another interesting deviation to the latitudinal trend is the earlier average ice-out at the Yukon Flats relative to Daring Lake (Figure 5c). Despite being 1.4 degrees south of the Flats, the Daring Lake sites are situated in the tundra with a lower mean annual temperature [66] than the boreal Yukon Flats sites [67]. The warmer annual temperature and consequent earlier ice-off at the Yukon Flats is consistent because ice-off is highly predicted by air temperature [68]. Green-up was more variable, with the average peak green-up occurring on July 29th (Figure 5b). Peak green-up between years varied by 22 days (Figure 5b). The earliest One possible source of this difference is hydrologic connectivity, or the degree to which a lake is connected to nearby surface waters (i.e. rivers, streams). Lakes with high hydrologic connectivity, such as those in the PAD [64], may experience intense pulses of sediment in the spring during ice-out. Sediment has high reflectance in the near-infrared and may confound ice-detection [65]. Sediment could explain the negative model bias because it also reflects in the visible bands, thus making water appear "bright" like ice to the human eye even though ice may be gone. Another interesting deviation to the latitudinal trend is the earlier average ice-out at the Yukon Flats relative to Daring Lake (Figure 5c). Despite being 1.4 degrees south of the Flats, the Daring Lake sites are situated in the tundra with a lower mean annual temperature [66] than the boreal Yukon Flats sites [67]. The warmer annual temperature and consequent earlier ice-off at the Yukon Flats is consistent because ice-off is highly predicted by air temperature [68].
Green-up was more variable, with the average peak green-up occurring on 29 July (Figure 5b). Peak green-up between years varied by 22 days (Figure 5b). The earliest green-up date was observed in Lake M2 near Yellowknife on 9 June 2017 and the latest green-up occurred on 20 August 2019 at Lake M10 near Yellowknife. Green-up showed a weaker relationship to latitude than ice-out as indicated by a lower r 2 value. (Figure 5c,d).
Peak green-up occurred on average 72 days (>2 months) after ice-out. This greening up was twice as fast (30 days faster than the mean) at the higher latitude Daring Lake sites (Figure 5d). Peak green onset the fastest in the sites with the latest ice-out (Figure 6a). This phenomenon also varied with latitude ( Figure 6b). However, we caution against over-interpreting these trends as independent validation data is not available to verify them. green-up date was observed in Lake M2 near Yellowknife on 9 June 2017 and the latest green-up occurred on 20 August 2019 at Lake M10 near Yellowknife. Green-up showed a weaker relationship to latitude than ice-out as indicated by a lower r 2 value. (Figure 5c,d).
Peak green-up occurred on average 72 days (>2 months) after ice-out. This greening up was twice as fast (30 days faster than the mean) at the higher latitude Daring Lake sites (Figure 5d). Peak green onset the fastest in the sites with the latest ice-out (Figure 6a). This phenomenon also varied with latitude ( Figure 6b). However, we caution against overinterpreting these trends as independent validation data is not available to verify them. Figure 6. (a) relationship between the ramp up, defined as the number of days between ice-out and peak green, and the ice-out day of year (type II regression, y = −1.65x + 307.06, r 2 = 0.35) (b) relationship between the number of days since iceout until peak green with latitude (type II regression, y = −6.95x + 486.59, r 2 = 0.29).

Climate Patterns and Processes
Mean annual temperature across all sites and years was −1.7 °C (standard deviation 16 °C) and mean annual precipitation, calculated as the average of the annual sum at each site averaged across sites and years, was 104 cm. Mean annual precipitation was 106, 87, 105, and 117 cm respectively for the PAD, Yukon Flats, Yellowknife and Daring Lakes. To investigate seasonal dynamics, we investigated relationships between ice-out and green-up and spring (March, April, May) and summer (June, July and August) dynamics. Ice-out showed stronger coherence with spring temperatures (Figure 7a). In contrast, green-up showed very little coherence with summer temperatures (Figure 7b). Interest- Figure 6. (a) relationship between the ramp up, defined as the number of days between ice-out and peak green, and the ice-out day of year (type II regression, y = −1.65x + 307.06, r 2 = 0.35) (b) relationship between the number of days since ice-out until peak green with latitude (type II regression, y = −6.95x + 486.59, r 2 = 0.29).

Climate Patterns and Processes
Mean annual temperature across all sites and years was −1.7 • C (standard deviation 16 • C) and mean annual precipitation, calculated as the average of the annual sum at each site averaged across sites and years, was 104 cm. Mean annual precipitation was 106, 87, 105, and 117 cm respectively for the PAD, Yukon Flats, Yellowknife and Daring Lakes. The To investigate seasonal dynamics, we investigated relationships between ice-out and green-up and spring (March, April, May) and summer (June, July and August) dynamics. Ice-out showed stronger coherence with spring temperatures (Figure 7a). In contrast, greenup showed very little coherence with summer temperatures (Figure 7b). Interestingly, the year with the warmest spring temperature (2017, mean spring T −2.3 deg C) also had the shortest distance from ice-out to green-up. In this warmer year, the ramp up period was 18 days shorter than in the cooler years of 2018 and 2019. This results support the hypothesis that increasingly warmer springs will lead to more rapid onsets of biological production [69].
Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 16 Figure 7. Ice-out day-of-year relative to spring temperature (a) and peak green day-of-year relative to summer temperatures (b).

CubeSat Satellites Provide Critical Insight Into Small Lake Ice-Out
Our study shows that it is possible to construct high resolution time series (daily, <30m) of lake ice-out and green up from PlanetScope imagery (Figure 3). Other studies have pointed out the error associated with PlanetScope reflectance relative to other sensors [70,71]. To explore the quality of these time series, we compared PlanetScope reflectance to Sentinel-2 reflectance (Figure 4). We discovered that while there was an average absolute difference of 36%, the shapes of the time series curves were similar. Sentinel-2 values were consistently lower than PlanetScope. This is likely because the cloud mask from PlanetScope imagery is not as restrictive as the Sentinel-2 masks. The difference in acquisition times between PlanetScope and Sentinel might be playing a role as lake surfaces or the atmosphere might have undergone changes. Ultimately, however, both products reconstruct ice-out in these relatively small lakes, even if there are differences in the absolute reflectance.

Ice-Out Showed Stronger Coherence with Temperature Than Green-Up
In comparing climate reanalysis data ( Figure S1) to ice-out patterns, we find a coherence between ice-out and air temperature gradients (Figures 5 and 7) consistent with other studies showing that variation in ice-out timing is controlled primarily by air temperature, which varies along latitudinal gradients [72][73][74]. Ice-out is primarily a physical process driven by temperature. In contrast, changes in the magnitude of green reflectance are primarily driven by the interaction of both physical and biological processes and indicate

CubeSat Satellites Provide Critical Insight Into Small Lake Ice-Out
Our study shows that it is possible to construct high resolution time series (daily, <30 m) of lake ice-out and green up from PlanetScope imagery (Figure 3). Other studies have pointed out the error associated with PlanetScope reflectance relative to other sensors [70,71].
To explore the quality of these time series, we compared PlanetScope reflectance to Sentinel-2 reflectance (Figure 4). We discovered that while there was an average absolute difference of 36%, the shapes of the time series curves were similar. Sentinel-2 values were consistently lower than PlanetScope. This is likely because the cloud mask from PlanetScope imagery is not as restrictive as the Sentinel-2 masks. The difference in acquisition times between PlanetScope and Sentinel might be playing a role as lake surfaces or the atmosphere might have undergone changes. Ultimately, however, both products reconstruct ice-out in these relatively small lakes, even if there are differences in the absolute reflectance.

Ice-Out Showed Stronger Coherence with Temperature Than Green-Up
In comparing climate reanalysis data ( Figure S1) to ice-out patterns, we find a coherence between ice-out and air temperature gradients (Figures 5 and 7) consistent with other studies showing that variation in ice-out timing is controlled primarily by air temperature, which varies along latitudinal gradients [72][73][74]. Ice-out is primarily a physical process driven by temperature. In contrast, changes in the magnitude of green reflectance are primarily driven by the interaction of both physical and biological processes and indicate ecosystem level change. The weaker relationship between green-up and latitude suggests that ecological productivity in these lakes is controlled in part by local and catchmentmediated processes such as microbial community composition, hydrologic connectivity, nutrient and water chemistry composition [18]. However, the higher latitude sites had the shortest time between when ice-out occurred and when greenness peaked (Figure 6b). This suggests that biological activity ramps up more quickly at higher latitudes, when ice-free season is shorter. In other words, the later the ice-out and the shorter the growing season, the faster biological activity ramps up. This has implications for total lake carbon cycling as longer ice-out periods have been shown to result in higher total chlorophyll-a in alpine lakes [75]. However, because these lakes are also at higher latitudes, we cannot conclusively rule out the effect could also be from latitudinal differences in photoperiod. While this paper can identify shifts in color, work remains to tease out the biological mechanisms from on the ground. In addition, while these trends may be coherent across large spatial scales, the relationship may vary at the individual lake scale. Individual lake characteristics such as depth, littoral zone extent, perimeter, surface area, and hydrologic connectivity exert control at the local level [76].
Here we present phenological 'fingerprints' using CubeSat data that can provide information on how the ice and productivity dynamics of small lakes change across heterogeneous landscapes. However, more research is needed to determine the mechanisms driving each individual lake phenology. We suggest that the summertime increases in greenness are associated with both plant and algal productivity within the lake once conditions are most conducive to primary production [24]. It is important to emphasize the range of peak green dates ( Figure 5). For example, peak greenness in one of our study regions, the Peace Athabasca Delta (PAD), occurred anywhere from June 30th to August 28th, covering nearly the entire summer growing season despite not varying in latitude by more than 4 degrees (Figures 1 and 5). It has been shown that proximity to riverine networks and floodplains can have strong coherence with lake hydrologic patterns within the PAD [77]. The stochastic nature of spring flooding events driven by the presence of ice jams along the Peace River flowing north may impose the largest influence on lowland lake phenology [64]. The frequency of ice jam flooding events within the PAD is decreasing, attributed to long-term decreases in precipitation and upper basin water management [78]. Changes in flow will also influence floodplain lakes in the Yukon Flats adjacent to the mainstem of the Yukon River as well, while the more isolated lakes in the Yellowknife/Daring Lake region may remain more stable. Developing a longer-term time series moving forward is critical to identify how hydrology will influence lake patterns and processes identified here.
These high-resolution time series revealed another interesting phenomenon: a double maximum in greenness in many sites ( Figure S5, Egg Lake for example). This is consistent with research conducted in Northwest Territory tundra ponds that showed two peaks in phytoplankton biomass and productivity driven by an early peak from diatoms followed by a later bloom driven by Chlorophyta [79]. A similar double-maxima was also observed over a 14 year period in a Finnish boreal lake driven first by a peak in Cryptomonas species followed by a peak in Gonyostomum semen [76]. While these reflectance spectra show a similar two-peak pattern, it must be emphasized that these spectral time series should not be over-interpreted in the absence of field measurements. It is important to note that manual validation was only available for ice-out because changes in greenness are too subtle to detect with the human eye using visual inspection. Future in situ validation of green-up is possible, for example by using continuous oxygen sensors and should be prioritized in field data collection. Field studies are needed to link and validate the remote sensing signatures we have evaluated here to physical and ecological properties on the ground in lakes.

Future Research Directions
New near-daily remote sensing imagery offers a new approach for monitoring at high spatial and temporal resolutions. This study demonstrates that these dense time series observations can be used to detect ice-out and peak green-up in water bodies as small as 10 ha (Table 1). When coupled with high resolution sampling methods, these observations represent a powerful new frontier for freshwater monitoring. For example, advances in continuous oxygen sensors [2] and high throughput imaging systems such as the Imaging FlowCytobot [3] have revolutionized our ability to track phytoplankton dynamics at finer temporal scales. Previously, the sheer number and remote nature of Arctic and boreal lakes has precluded the widespread characterization of lake ecological dynamics with these new techniques. Combining high resolution field techniques with high resolution remote sensing will revolutionize our ability to track lake ecological dynamics at unprecedented scales. Satellite remote sensing provides a powerful set of observations for insights into lake ecological trends.
Climate-driven temporal change in lake ecological patterns has only been recorded through paleolimnological records that involved site specific physical sampling of sediments, and significant investment in personnel [34,77]. Through high temporal frequency remote sensing, we provide a framework to better identify the influence of both management and changing climate on water resources for a subset of arctic-boreal lakes in western North America. Future studies leveraging satellites observations to examine the millions of arctic-boreal lakes could further illuminate major controls. Future in situ work to validate Cubesat-derived observations will also be essential. Remote sensing observations can inform field campaigns, which are then crucial for quantifying the ecological processes that are taking place at a mechanistic level. These two approaches should not be considered as separate efforts, but highly complementary and necessary companions for advancing our understanding of global change impacts on lake ecosystems.

Conclusions
High spatiotemporal resolution CubeSat sensors significantly improve our ability to monitor freshwater ecosystems from afar. A key need exists to demonstrate the utility of these emerging observations. In this study, we use time-series imagery to generate high resolution phenology of lake color. We are able to successfully predict lake ice-out and show links between climate and ice-out patterns. Finally, we compare CubeSat and Sentinel-2 surface reflectance products and show a strong correlation between the two sensors. Our findings show ice-out and green-up vary more over space than through time and that lakes with earlier ice-out have later green-up. This finding is significant as ice-out has advanced significantly with warming with likely consequences for the timing of productivity. The global coverage of PlanetScope imagery suggests the methods demonstrated here have the potential to be applied to lakes across northern latitudes to identify ice-out, changes in the growing season length and changes in green reflectance that have ecological significance.