Remote sensing of lake ice phenology in Alaska

The timing of lake ice breakup and freezeup are important indicators of climate change in Arctic and boreal regions because they respond rapidly and directly to variations in climate conditions. Despite its importance, lake ice phenology remains poorly documented in most lakes of Alaska. To fill this data gap, we constructed a remote sensing-derived lake ice phenology database covering all lakes in Alaska larger than 1 km2 (n = 4241) over the period 2000–2019. This dataset, which includes lake ice on/off dates and lake ice duration, was based on an automatic method using daily moderate resolution imaging spectroradiomenter (MODIS) imagery to measure lake ice fraction. This method extracts lake ice pixels from MODIS images using a dynamic threshold that was calibrated against Landsat Fmask. Different filters that account for clouds, polar night, and other sources of error were applied to increase the accuracy of lake ice phenology estimation. Trend analysis shows earlier breakup (−5.5 d decade−1) for 440 lakes and later breakup (7.5 d decade−1) for four lakes (p < 0.05). A total of 289 lakes had significant trends toward later freezeup (2.9 d decade−1) and 11 lakes towards earlier freezeup (−3.3 d decade−1). Most lakes with significant trends are north of the Brooks Range. This dataset can contribute to increased understanding of interactions between lake processes and climate change, and it supports the study of biogeochemical, limnological and ecological regimes in Alaska and pan-Arctic regions.


Introduction
Lakes cover approximately 2% of the Earth's land surface globally (Adrian et al 2009), and as much as 8% in Arctic and boreal regions (Watts et al 2012). Lake properties such as area variation, water temperature, dissolved organic carbon and plankton composition, which respond rapidly and directly to environmental and climate conditions, can serve as important indicators of climate and environmental changes (Adrian et al 2009). This is particularly true in Arctic and boreal regions, where climate is changing rapidly (Vincent et al 2013). Among these properties, lake ice phenology is a particularly robust proxy for climate variability. Climate and environmental changes (e.g. air temperature) contribute to the variability of lake ice formation and decay (Palecki andBarry 1986, Robertson et al 1992). Since temperature is the primary driver of lake ice variations (Šmejkalová et al 2016), earlier breakup and later freezeup would be expected under a warming climate (Sharma et al 2020). In addition to its use as a proxy for climate variability, the development and loss of ice is important because of its substantial impact on lakes themselves and nearby terrestrial systems. Many processes happening in or near lakes are affected by the phenology of lake ice. For example, an increased ice-off period would result in a larger direct heat transfer between lakes and the atmosphere, export of more greenhouse gases due to chemical and biological processes, and increased annual evaporation from the lake surface (Latifovic andPouliot 2007, Adrian et al 2009). The presence of ice in lakes also alters the specific light condition and the turbulence in the water column, both critical to lacustrine plankton communities (Adrian et al 1999).
Studies focusing on the role of lake ice phenology in climatology and ecology suggest that lake ice conditions, especially as measured by freezeup and breakup dates, reflected changes to climate (Sharma et al 2019(Sharma et al , 2020. Magnuson et al (2000), for example, found significant trends towards an earlier breakup and later freezeup. From 1846 to 1995 the mean shift in freezeup dates was 0.58 d decade −1 later and changes in breakup dates averaged 0.65 d decade −1 earlier based on the historical ground observations for 39 lakes and rivers in the northern hemisphere. Latifovic and Pouliot (2007) analyzed 54 year trends of breakup and freezeup dates in 42 large Canadian lakes and found that the majority showed earlier breakup and later freezeup dates. The most significant trend was detected for lakes in the far north for the last decade of the study period, with an earlier breakup trend of 9.9 d decade −1 and later freezeup trend of 7.5 d decade −1 . Šmejkalová et al (2016) established a lake ice breakup phenology dataset for ∼13300 Arctic lakes via remote sensing. The changes in breakup dates ranged from −1.0 to −10.5 d decade −1 , depending on the region. All these studies agree well with the common expectation that there would be a shift towards earlier breakup and later freezeup dates in the context of global warming, which indicates that timing of lake ice can be used to infer environmental and climate changes.
Over the past century, an increase of 1.4 • C (compared to 0.8 • C worldwide) was detected for Interior Alaska by analyzing a 100 year climate record from Fairbanks (Wendler and Shulski 2009). Because boreal ecosystems are highly sensitive to temperature, this significant warming trend has important implications for both biological and physical processes in the region (Stone et al 2002). There has been considerable work on understanding lake ice conditions and variability in Alaska. For example, modeling approaches based on energy budget have been used to simulate the growth and decay of lake ice in Alaska (Jeffries et al 1996, Zhang and Jeffries 2000, Morris et al 2005. However, lake ice models usually simplify the physical processes governing lake ice formation and melt, and large uncertainty may exist in climate reanalysis and outputs of climate models that are used to force the lake ice model (Xue et al 2017). Therefore, in situ data is still important for calibrating and validating such models (Zhang and Jeffries 2000, MacKay et al 2017, Xue et al 2017. Even though in situ observations are critical for monitoring the impacts and spatial patterns of climate change, they remain sparse in remote Arctic and boreal regions. Remote sensing offers an alternative for directly observing lake ice phenology in Alaska. Commonly used remote sensing approaches for detecting lake ice breakup and freezeup can be divided into categories of passive microwave, active microwave and optical, based on the sensors used. Passive microwave sensors, which are not susceptible to cloud cover, can be used for mapping ice extent or measuring ice phenology over large rivers or lakes (Brakenridge et al 2007, Che et al 2009). However, the coarse spatial resolution of derived ice products (∼5-75 km) limits their use to the largest water bodies (Beitsch et al 2014, Zhang and Gao 2016). Active sensors such as synthetic-aperture radar (SAR) can also be used to monitor lake ice dynamics such as breakup and freezeup (Floyd et al 2014) at high spatial resolution (1-30 m) (Zhu and Bamler 2010). But the low temporal resolution (5-6 d) of SAR sensors can result in a large uncertainty, which limits its application for climate studies. Additionally, the intensity of backscatter in SAR imagery, which is used for differentiating ice from water, can be affected by the roughness of the surface. Therefore, it may be difficult for SAR to distinguish between ice and open water when the ice is wet and smooth (Wakabayashi et al 1993, Lundhaug 2000. Optical satellite imagery is also commonly used for detecting lake ice extent. Landsat images (30 m spatial resolution) can be used for monitoring ice effectively based on the different reflectivity characteristics of ice and water (Nolan et al 2002, Duguay et al 2003, Korzeniowska and Korup 2017. Despite the advantage of Landsat's high spatial resolution, the low effective temporal resolution, especially in cloudy regions, limits its monitoring capability. Moderate spatial resolution sensors such as moderate resolution imaging spectroradiometer (MODIS) have been broadly applied to detecting lake ice phenology because of their high temporal resolution (e.g. twice daily for MODIS). Arp et al (2013) (Šmejkalová et al 2016), rather than the full spatiotemoporal dynamics of ice extent. To our best knowledge, no comprehensive lake ice phenology dataset which includes all the lake ice phenology timing (i.e. breakup, freezeup, ice duration) covers the entirety of Alaska. To fill the gap, we develop an approach to estimate lake ice phenology in Alaska that contains all Alaskan lakes larger than 1 km 2 . By applying a polar-night-correction operation, outliers caused by low solar radiation are removed, which results in an enhanced accuracy for the estimates. This database, which includes lake ice breakup date, freezeup date and ice duration, is validated against independent lake ice phenology data. Using our lake ice dataset, we analyze spatial pattern of lake ice phenology and quantify corresponding temporal trends. We also seek to robustly quantify the uncertainty associated with these trends by propagating error regarding the impact of clouds on remotely sensed breakup/ freezeup detection into the trend analysis.

Material and methods
The first aim of this study is to develop a remotely sensed lake ice phenology dataset covering nearly all lakes in Alaska larger than 1 km 2 . For this purpose, we applied a method developed by Zhang and Pavelsky 2019 to calculate lake ice fraction. This method classifies ice pixels based on the red reflectance band of MODIS imagery, with an ice detection threshold calibrated against ice fraction from Landsat Fmask (Zhu and Woodcock 2012) over each lake. The time series of lake ice fraction is then used to identify lake ice breakup and freezeup dates. Multiple filters that account for clouds, polar night, and other sources of error were applied to increase the accuracy of lake ice phenology estimation.

Selected lakes
We obtained the boundaries of all lakes in Alaska from the Alaskan Lake Database Mapped from Landsat Imagery (Wang 2011). This Landsat-based dataset contains lakes in Alaska larger than 0.1 km 2 . Considering the moderate spatial resolution of MODIS imagery (250 m), only lakes with surface area >1 km 2 were selected in this study to increase the accuracy of the lake ice phenology detection. Smaller lakes, which are more likely to be primarily covered by mixed pixels, were not considered. In total, 4241 of 294 274 lakes were chosen, accounting for 61% of lake surface area in Alaska. Figure 1 shows the geographical locations of the lakes analyzed.

Satellite data 2.2.1. MODIS reflectance data
In this study, red reflectance data (620-670 nm, band 1) from the MODIS/Terra Surface Reflectance Daily product (MOD09GQ, collection 6) (Vermote and Wolfe 2015) collected between 2000 and 2019 was used to distinguish ice pixels from water over the study lakes. MOD09GQ provides two bands (red and near infrared) at a 250 m resolution in the Sinusoidal projection. It is a daily estimate of the surface spectral reflectance after atmospheric corrections. We used MOD09GQ band 1 (red) for two reasons. First, it has the finest spatial resolution among all MODIS products (along with band 2). Also, as vegetation has a lower reflectivity in the red band than the NIR band, using the red band to detect lake ice status can better remove the impact of changing vegetation. Second, the high temporal resolution of MOD09GQ makes it suitable for identifying lake ice breakup or freezeup events at a daily timescale.

MODIS cloud mask
To identify and remove clouds, MODIS cloud flags were extracted from the MOD09GA data product. MOD09GA contains several data quality control measurements that can be used in masking pixels with quality issues or assessing the uncertainties. Among these quality indicators, the state_1km_1 band provides the standard cloud mask, which includes cloud shadows, aerosol quantity, and cirrus detection at 1 km spatial resolution. The state_1km_1 is stored in an efficient bit-encoded manner with 16 bits. Each bit indicates a different meaning of cloud related status, such as cloud cover or cloud shadow. The cloud flags were applied to filter out MOD09GQ redimages with heavy cloud cover.

Landsat fmask
We use Landsat Fmask (Zhu and Woodcock 2012) to calibrate the reflectance threshold of MODIS images for ice classifications. The Fmask algorithm was originally developed for masking cloud, cloud shadow, and snow for Landsat TM and ETM+. Fmask classifies pixels into different groups for each individual image using Landsat top of atmosphere reflectance and brightness temperature as inputs based on an object-based cloud and cloud shadow matching algorithm. It has two advantages as reference data to calibrate the reflectance threshold for classifying MODIS images. First, compared with MODIS imagery, Landsat has much higher spatial resolution (30 m), which allows it to provide monitoring results with higher confidence under cloud free conditions. Second, the reliability of Fmask has been verified in multiple studies (Shao et   Note that multi-source remote sensing and reanalysis data with different spatial resolutions were used in this study. For processes that combine different spatial resolutions, we downscaled images with lower resolutions to higher resolutions, by assigning the same value to each subpixel.
We identify all Landsat images overlapping with the study period and compared lake ice fraction values based on Landsat to same-day values based on MODIS.

Surface air temperature
To inspect the relation between surface air temperature and lake ice phenology, we used temperature data from the NCEP/NCAR Reanalysis Project, which is a joint project between the National Centers for Environmental Prediction (NCEP, formerly 'NMC') and the National Center for Atmospheric Research (NCAR). The NCEP/NCAR Reanalysis project is using a state-of-the-art analysis/forecast system to perform data assimilation using past data from 1948 to the present. The data have 6 h temporal resolution and 2.5 degree spatial resolution.

Multi-error-removed improved-terrain (MERIT) DEM
The surface elevations of lakes were extracted from MERIT DEM, which is a high-accuracy global DEM at 3 arc second resolution (∼90 m at the equator) produced by eliminating major error components from existing DEMs. MERIT DEM separates absolute bias, stripe noise, speckle noise and tree height bias using multiple satellite datasets and filtering techniques. After the error removal, land areas mapped with 2 m or better vertical accuracy were increased from 39% to 58%.

In situ lake ice phenology data
Most of the in situ lake ice phenology data were obtained from studies led by the USGS and the University of Alaska Fairbanks (UAF). Surface and bed temperature data collected from autonomous thermistors typically moored near lakes centers were used to indicate timing of ice-cover formation (freezeup) and ice-cover loss (breakup) as described by (Arp et al 2013(Arp et al , 2016. In situ sensor freezeup and breakup data used for validation in this study are reported to the nearest day and the exact sensor model and configuration vary among separate studies (

Lake ice phenology calculation
The method for identifying lake ice breakup and freezeup timing generally contains five major steps (figure 2). At the very beginning, because cloud obscuration is the main limitation that hampers the use of optical sensors for ice cover detection and mapping, we removed cloudy pixels in MOD09GQ reflectance images using the cloud mask band 'state_1km_1' of MOD09GA. Second, we applied the method developed by Zhang and Pavelsky (Zhang and Pavelsky 2019) to classify lake ice pixels. Pixels with red band reflectance values greater than the dynamic threshold were classified as ice covered. The threshold for classifying ice vs. water was selected by calibrating the lake ice fraction against that from Landsat Fmask because we assumed that Landsat Fmask more accurately classifies lake ice cover due to its finer spatial resolution (figures S1 and S2 (available online at stacks.iop.org/ERL/16/064007/ mmedia)). After extracting lake ice pixels, ice fraction over lake surface was calculated in the third step.
Fourth, we applied a polar-night-effect correction operation to reduce the influence from lower solar radiation during fall/winter over lakes in highlatitude regions. Lower solar radiation during the polar night period is likely to be associated with lower reflectance values in red band satellite imagery. This reduced reflectance results in a misclassification between ice and water because the calculation of reflectance from radiance is less accurate when solar radiation is low. For a given location at the high latitudes, images incorporated into the daily MOD09GQ composite might be assigned from the observations of different orbit passes, which means not all images are acquired at the same time of day. This inconsistency can result in different illumination conditions from one daily composite to the next. An example of misclassification caused by polar night effect is presented in figure 3. At point 2, the remotely sensed ice fraction  dropped from 100% to 0% due to weak illumination conditions. To address this problem, we assume that lakes ice fraction does not decrease during winter when radiation is low. We adjust the ice fraction of cases like that shown in figure 3(b) to 100% if two criteria are met. First, the time span from sun rise to sun set is smaller than 7 h. Second, the previous ice fraction observation was larger than 80%, indicating an already frozen lake. Note that even though freezeup and breakup happen during days with sufficient daylight, correcting the low lake ice fraction during polar night is still necessary. This is because erroneous breakup and freezeup dates are sometimes detected by remote sensing algorithms when there is insufficient solar radiation (like polar night or influence by shadow), or the reflectance is extremely high due to clouds that are not detected by the cloud filter. The outlier control method in this study aims to select the longest ice period, which means the falsely detected breakup and freezeup dates can affect the performance of quality control operations. Figure S3 shows the comparison of lake ice fraction over different lakes after applying each filter, from which we can see the outliers present in the initial time series have been substantially reduced after applying the quality control and polar night effect corrections.
Last, we used the time-series of lake ice fraction to identify the timing of freezeup and breakup. We applied a threshold of 20% lake ice fraction to identify the lake freezeup and breakup dates. When the lake ice fraction first reached 20% each fall, we identified the freezeup date. The breakup date was defined as the day each spring when the lake ice fraction dropped below 20%. The definition of lake ice breakup and freezeup varies in different studies. For instance, Duguay et al (2006) used the mean reflectance value over lake surface to identify lake breakup/ freezeup dates. For Spencer et al (2008), freezeup was defined as when the lake ice fraction first reaches 90%. In this study, we chose 20% as the threshold because for most of the lake ice classifications (99.2%), the errors in ice detection, mainly caused by misclassification of cloud as ice, range from 0% to 20%. Lake ice fraction larger than 20% is more likely to be associated with real ice events. We calculate the lake ice duration for each year using the time span between the freezeup date and its following breakup date. Note that though the MODIS Cloud Mask product can be used effectively to remove cloud-affected pixels, in some cases clouds are misclassified as ice, which can result in ice detection in the middle of summer. We used additional two-step quality control operation to address this cloud residual problem (supplemental).

Estimation of freezeup and breakup date uncertainty
Breakup/freezeup uncertainties here were the measurement of missing data caused by clouds. In many cases, breakup or freezeup occured during a cloudy period, and we were not certain of the exact date. We only knew the general range, bounded by the dates of clear-sky images, when the breakup/freezeup processes have already occurred. In these cases, we considered the ice breakup or freezeup event to have happened on the middle day of the cloudy period between these two images. The uncertainty in each breakup or freezeup date was then the number of days between this middle day and the first day of this missing data period.

Trend analysis
Breakup and freezeup dates with uncertainties (calculated in section 3.4.2) larger than 10 d were removed when performing trend analysis. When calculating the trends, lake ice breakup/freezeup dates in certain years were also ignored if multiple breakup/ freezeup occurred. A lake having more than 2 years of such multiple breakup/freezeup are removed for trend analysis. In total, 343 lakes are ignored for the trend analysis. Trends of freezeup and breakup dates for the lakes are calculated using the remote sensing observations for the period 2000-2019. The magnitude of slope is estimated using Sen's slope (1968), which is a non-parametric test. By pairing all the ordered lake ice breakup or freezeup dates within the study period, multiple slopes can be calculated. Sen's slope is defined as the median value of all possible slopes. The significance of each trend is examined using the Mann-Kendall test (p < 0.05), a rank-based non-parametric test for trend (Mann 1945, Kendall 1975).

Validation of remotely sensed lake ice phenology
The comparison of breakup dates with in-situ data (figure 4) shows that the breakup dates were detected with relatively high accuracy. The overall mean absolute error (MAE), root mean squared error (RMSE), correlation coefficient (r) and bias were 5.0 d, 7.49 d, 0.84 and −1.13 d, respectively for the proposed algorithm. The negative biases indicate that the breakup dates detected by remote sensing methods were slightly earlier than the in-situ observations. This might be caused by different definitions of breakup dates between the remote sensing algorithm and ground-based observer. In addition, thin ice or ice with low albedo due to melt season processes might be detected as water by MODIS during the breakup period, leading to a negative bias of breakup estimation. When detecting freezeup, the overall MAE, RMSE, r and bias of the algorithm used in the study were 10.72 d, 13.74 d, 0.59 and −0.04 d, respectively (figure 4). Compared to breakup, freezeup detection resulted in a larger MAE and RMSE. This may be caused by heaver cloud coverage during the freezeup period. An outlier in our validation was found in breakup for Lake Minchumina. The breakup date from remote sensing in 2007 was the 177th day while the in-situ breakup date was the 140th day. This large difference was because the cloud filter failed to remove cloudy pixels during the ice melt period. Therefore, the remote sensing algorithm mistook the date when the cloud disappeared as the breakup date. In addition, we found that the bias was more negative for lakes in southwest Alaska with later freezeup dates. This was because the in-situ data in southwest lakes were collected from a different source. The inconsistent approaches of obtaining in-situ data might be a source of apparent differences in performance when evaluating the algorithm. To better understand the robustness of the results, such as how cloud cover influences the calculation of breakup/freezeup dates and the corresponding trends, visual inspection and uncertainty analysis are provided in supplemental materials (figures S4-S7). Shown as figure S4, the algorithm performs even better for smaller lakes than large lakes in terms of MAE and NRSME. This is because the threshold for ice classification in this algorithm was calibrated with Landsat Fmask which constrained the effect of mixed pixels. Furthermore, the in situ data of large lakes may be biased due to limitations of sensor detection area relative to much larger lake surface area. For example, the portion represented by the in situ sensor of a lake may be ice free even if centrally located, while more distant portions may be still ice-covered. This situation would be more common in large lakes with slower ice melt, which are also subject to wind redistribution of ice pans to the wind-ward shoreline.

Spatial pattern of lake ice phenology
Examination of patterns in lake ice breakup and freezeup timing across Alaska clearly shows individual and regional differences, with breakup and freezeup dates gradually changing from south to north (figure 5). The average value of lake ice breakup dates for all lakes was 28 May (day 148). To measure the variability of breakup/freezeup over each lake, we used the 5th and 95th percentiles instead of the earliest and latest breakup/freezeup. This is because some of the earliest/latest lake ice phenology was probably caused by outliers due to pervasive cloud effects. To identify the 5th and 95th percentiles, we sort the mean breakup dates of all lakes in ascending order. The 5th percentile mean lake ice breakup date was 14 April (day 104) in the southwest coastal region and the 95th percentile was 24 June (day 175) in the northa 71 d period separating lake breakup dates. Similar spatial patterns were detected for lake freezeup dates and ice durations. Lakes located in southern Alaska were more likely to be associated with later freezeup dates and short ice durations. The 5th and 95th percentile mean freezeup dates were 26 September and 23 November (day 269 to day 327). For lake ice durations, the 5th and 95th percentiles were 155 and 271 d. An uncertainty analysis to quantify the influence of missing data caused by cloud cover showed that the average uncertainties of lake ice breakup and freezeup dates were 1.9 and 3.0 d, respectively ( figure S6).
To illustrate what drivers might affect lake ice phenology, we compared different factors with mean lake ice breakup/freezeup dates. As expected, mean observed lake ice phenology was strongly correlated with latitude ( figure 6). The coefficient of determination (R 2 ) values between latitude and lake freezeup timing variation was 0.88. The R 2 value was 0.75 between latitude and breakup timing. The mean annual NCEP temperature were also highly correlated with lake ice phenology, with R 2 values of 0.67 and 0.82 for breakup and freezeup. However, we hardly found any significant relationship with altitude and lake surface area. The R 2 for altitude and lake surface area ranged from 0.001 to 0.014, indicating that the latitude and mean air temperature are more dominant on controlling lake ice phenology at large spatial scale in Alaska.

Trend analysis 3.3.1. Spatial pattern of lake ice phenology trends
The average rate of change in breakup timing was −2.2 d decade −1 . A total of 440 lakes (10.37% of all lakes) had a statistically significant earlier trend (−5.5 d decade −1 on average) at the p value of 0.05 (figure 7). Only four lakes show a statistically significant later breakup trend. If we only inspect the direction of the trends (both significant and insignificant), 88.78% lakes show earlier breakup trends. In terms of trend magnitude and trend significance, there was substantial region-to-region variability within Alaska. The significant earlier breakup trend was more evident in northern Alaska. A few lakes with significantly later breakup were also detected in the north-central region. Results reported by Šmejkalová et al (2016) for 2000-2013 using remote sensing observations also show a strong trend towards earlier breakup in far northern Alaska (∼−9 d decade −1 ). Analysis of trend against spatial domain revealed different trend behaviors. From the central to the north part of Alaska, trends generally increased, but  trends are most apparent for lakes north of the Brooks Range. On average, trends in Northern Alaska (−3.8 d decade −1 ) were 4.1 times larger than in the Yukon River Basin (−0.93 d decade −1 ), where little change was found. A total of 289 lakes had significant trends toward later freezeup (2.9 d decade −1 ) and 11 lakes towards earlier freezeup (3.3 d decade −1 ). Ignoring statistical significance, 79.8% lakes exhibit later freezeup trends. The smaller number of lakes with statistically significant trends in freezeup compared to breakup may reflect the fact that our data are noisier during the freezeup period. The average duration trend of all lakes was −0.8 d decade −1 (figure 6). A total of 195 lakes had significant trends (p < 0.05), with 189 lakes exhibiting negative trends (average −2.1 d decade −1 ) and six lakes showed positive trends (1.5 d decade −1 ). The uncertainties of breakup and freezeup trends for all lakes caused by clouds were 1.5 and 1.2 d decade −1 (figure S7).

Trend analysis by ecoregion
We examined the percentage of lakes that have significant lake ice timing trend at the scale of the ecoregion. An ecoregion is defined as an area of land and water containing vegetation communities that have the same ecological characteristics, such as species and ecological dynamics, environmental conditions, and interactions (Adrian et al 1999, Surdu et al 2014. These ecological features are critical for the long-term persistence of vegetation communities. Climate change is expected to affect variation in the seasonality of lake ice which has import and wide physical and biological ecosystem effects in the region (Reed et al 2009). Analyzing the variation in lake ice phenology across multiple lakes within an ecoregion can provide broader insights into the linkage between lake ice phenology and local ecology. For instance, the variations of lake ice would alter light availability in the water column at high latitudes, which has important consequences in the carbon balance and ecosystem behavior (Cory et al 2014). Knowing lake ice variability by ecoregions would be beneficial to understand the changes of ecosystems in the region and identify the ecoregions which are most vulnerable to climate change. Twenty ecoregions have previously been identified by synthesizing information on  the geographic distribution of environmental factors such as climate, physiography, geology, permafrost, soils, and vegetation (Gallant 1995).
The percentage of lakes with significant breakup trends across ecoregions suggest wide regional differences (figure 8). The maximum percentage value of lakes with significant breakup trends was detected in the interior highlands region (28%, n = 6/21). The interior highlands ecoregion is composed of rounded, low mountains, often surmounted by rugged peaks. The minimum percentage was detected in the Bristol Bay-Nushagak Lowlands, located in southwestern Alaska off Bristol Bay, and the Wrangell Mountains, which had no lakes with significant trends. The highest percentage of lakes with significant freezeup trends was detected in the Cook Inlet region. Out of 53 lakes in Cook Inlet, 18 lakes had significant trends towards later freezeup.

Discussion and conclusion
Variations of lake ice phenology can serve as a robust climate change indicator and are essential to understand interactions between lake ice and regional climate, environmental conditions and ecosystems (Rouse et al 2005, Adrian et al 2009. The changes of lake ice phenology may reflect the magnitude of climate change occurring in Alaska and be used for analyzing the effect of lake ice phenology on local ecosystems. Earlier breakup dates could increase the absorption of solar radiation and therefore accelerate the release of carbon from water to the atmosphere (Cory et al 2014). A comprehensive ice phenology dataset covering lakes in all sizes would be beneficial to better understand the fate of carbon in the Arctic. In addition, changes in ice duration may impact lake evaporation in high-latitude regions (Adrian et al 2009, Zhao andGao 2019).
In this study, we developed a first lake ice phenology dataset covering all lakes larger than 1 km 2 in Alaska. This dataset includes lake ice cover fraction and lake ice breakup and freezeup dates over 2000-2019. We explored spatial patterns in lake ice phenology, along with corresponding temporal trends. Unlike previous studies of lake ice phenology that examined only a few lakes (Duguay et al 2003, 2006, Arp et al 2013, Cai et al 2017 or focused only on ice breakup timing (Šmejkalová et al 2016), we used a well-validated method to measure breakup and freezeup timing, as well as ice duration, over thousands of lakes. As such, we provide a more complete picture of Alaskan ice phenology than do previous studies.
The spatial pattern of lake ice phenology in Alaska suggests that the mean lake freezeup timing, breakup timing and ice duration were highly correlated with latitude and temperature. The timing of lake ice phenology is governed by the surplus or deficit in the energy balance of the ice cover. The amount of solar radiation received at the land surface varies by latitude, the timing of ice phenology therefore depends on the latitudinal location of the lake (Morris et al 2005, Adrian et al 2009, Brown and Duguay 2010. For every degree rise in latitude, the average breakup date during the study period increased by about 4.1 d. The relationships with temperature are not surprising, though the fact that mean annual temperature explains such a large fraction of variability in freezeup timing (82%) shows its very strong importance. It is possible that a more detailed analysis of fall and spring temperatures would explain slightly more variability in freezeup and breakup dates, respectively. Meanwhile, there was a large variability over the entire domain in the spatial pattern of ice phenology timing. The difference of lake ice breakup timing between the 5th and 95th percentiles was more than 2 months over lakes in southwest and northern coastal regions, consistent with study by Arp et al (2013) that also found a 2 month lag in ice breakup timing across Alaskan lakes. The variability was primarily driven by temperatures, and might also be caused by some other factors such as elevation, lake surface area, and watershed area of individual lakes (Arp et al 2013). However, in this study, we found no strong relationship between lake ice phenology and lake elevation or lake surface area. Analysis of trends in lake ice timing shows that though ice breakup was occurring earlier and freezeup later for most lakes in Alaska, only a relatively small fraction of trends for individual lakes were significant. There were 440 lakes having significantly earlier breakup trends and 289 lakes having significant trends toward later freezeup. The correlation analysis shown in figure 6, and consideration of all the lakes, may limit detection of important patterns because lakes are not evenly distributed in Alaska. For example, while most lakes experiencing earlier breakup are located north of the Brooks Range, the region with the largest fraction of such lakes is the interior highlands, which only has 21 lakes (figure 8). It is also interesting to note the patterns suggesting that central Alaska may be experiencing the most robust trends in freezeup based on regional analysis, a pattern not as apparent in analysis of the unaggregated dataset.
This dataset still has some limitations. First, for regions with heavy clouds (e.g. the southern coast of Alaska), remotely sensed lake ice phenology have high uncertainties and may be biased due to the missing data. Besides this, due to the high reflectance of turbid water, the estimation of lake ice may be biased due to high turbidity . Second, the length of the data record is only 19 years, which is relatively short for examining trends. The short time series may be one of the reasons that only a small number of lakes have significant trends (Šmejkalová et al 2016). The fact that a large majority of lakes exhibit trends towards earlier breakup (∼88% of lakes) and later freezeup (∼79%), whether significant or not, suggests that a longer time series may, in the future, show significant trends for more lakes.
The dataset developed here can be used to support the study of biogeochemical, limnological and ecological regimes over large spatial scales. Given the high performance computing capabilities of Google Earth Engine (Gorelick et al 2017), an extended dataset including all lakes in the boreal region can be developed. In addition to direct scientific studies, such a dataset would also be useful for identifying ice covered lakes in other satellite datasets. For example, the Surface Water Ocean Topography (SWOT) mission will measure variations in area and water surface elevation of lakes as small as 0.0625 km 2 globally (Biancamaria et al 2016). As automated detection of ice with SWOT's Ka-band radar may be impossible, development of a global ice cover dataset from optical sources would help ensure the consistency of SWOT datasets.

Data availability statement
The data that support the findings of this study are available at https://github.com/tamuzhang/lake_ ice_AK.