An Automatic Method to Detect Lake Ice Phenology Using MODIS Daily Temperature Imagery

: Lake ice phenology is a climate-sensitive indicator. However, ground-based monitoring suffers from the limitations of human vision and the difﬁculty of its implementation in harsh environments. Remote sensing provides great potential to detect lake ice phenology. In this study, a new automated method was developed to extract lake ice phenology parameters by capturing the temporal pattern of the transitional water/ice phase using a parameterized time function. The method is based on Moderate-Resolution Imaging Spectroradiometer (MODIS) daily temperature products, which have unique potential for monitoring lake ice cover as a result of providing four observations per day at 1 km spatial resolution from 2002 to 2016. Three seasonally ice-covered lakes with different characteristics in different climate regions were selected to test the method during the period of 2002–2016. The temporal pattern of water/ice transition phase was determined on the basis of unfrozen water cover fraction extracted from the MODIS daily temperature data, and was compared with the MODIS snow and reﬂectance products and Landsat images. A good agreement with an R 2 of above 0.8 was found when compared with the MODIS snow product. The annual variation of extracted ice phenology dates showed good consistency with the MODIS reﬂectance and AMSR-E/2 products. The approach was then applied to nine seasonally ice-covered lakes in northern China from 2002 to 2016. The strongest tendency towards a later freeze-up start date was revealed in Lake Qinghai (6.31 days/10 yr) among the lakes in Tibetan plateau, and the break-up start and end dates rapidly shifted towards earlier dates in Lake Hulun ( − 3.73 days/10 yr; − 5.02 days/10 yr). The method is suitable for estimating and monitoring ice phenology on different types of lakes over large scales and has a strong potential to provide valuable information on the responses of ice processes to climate change.


Introduction
Lake ice phenology, which encompasses freeze-up and break-up periods/dates and ice cover duration, is important for the biological, chemical and physical processes of lakes, especially in cold regions [1,2]. Long-term records of ice phenology serve as sensitive indicators of climate variability [3]. Increases in the interannual variability in lake ice phenology have been found in recent studies [4][5][6]. Ice freeze-up occurs later and ice breakup earlier, with the result of ice cover being lost at an accelerating rate in the Northern Hemisphere [2,5,6]. These trends in ice phenology are closely related to increasing air temperatures [7][8][9] and to the integrated effects of changing wind speeds, snow cover, and solar radiation [10].
Despite their robustness, ground-based lake ice observations have experienced a global decrease since the 1980s [11]. In turn, satellite remote sensing, as an alternative tool, has played a growing role in lake ice monitoring [12]. The brightness temperatures obtained from microwave satellite sensors, such as the Advanced Microwave Scanning dimictic lake located on the northeastern Tibetan Plateau. It is the largest inland lake in China and is subject to a semiarid continental climate characterized by warm summers, cold winters, and more precipitation in summer than in other seasons [19]. Lake Ngoring is the highest (~4200 m asl) freshwater lake in China and is located in the source region of the Yellow River, surrounded by hills covered with alpine meadows. The south and east shorelines freeze and melt first [42]. Lake Hulun is a shallow monomictic lake located in the eastern part of the Mongolian Plateau in the northern part of China and is affected by semiarid continental and monsoon climates. The lake freezes from the complex shoreline area of the lake, and melting starts from the northwest coast and gradually toward the east coast [43].
Remote Sens. 2021, 13, x 3 of 23 chosen to develop, evaluate, and test the proposed approach ( Figure 1 and Table 1). Three of the lakes, Lake Qinghai, Lake Ngoring, and Lake Hulun, were selected as target lakes to develop and test the ice phenology extraction algorithm. Lake Qinghai is a large brackish dimictic lake located on the northeastern Tibetan Plateau. It is the largest inland lake in China and is subject to a semiarid continental climate characterized by warm summers, cold winters, and more precipitation in summer than in other seasons [19]. Lake Ngoring is the highest (~4200 m asl) freshwater lake in China and is located in the source region of the Yellow River, surrounded by hills covered with alpine meadows. The south and east shorelines freeze and melt first [42]. Lake Hulun is a shallow monomictic lake located in the eastern part of the Mongolian Plateau in the northern part of China and is affected by semiarid continental and monsoon climates. The lake freezes from the complex shoreline area of the lake, and melting starts from the northwest coast and gradually toward the east coast [43]. (D) Lake Namco; (E) Lake Zharinamco; (F) Lake Mapang; (G) Lake Bosten; (H) Lake Xinkai; and (I) Lake Hulun.  (D) Lake Namco; (E) Lake Zharinamco; (F) Lake Mapang; (G) Lake Bosten; (H) Lake Xinkai; and (I) Lake Hulun. The products use the split-window algorithm to retrieve surface temperatures from band 31 and band 32 [45]. The "land" surface temperature products include both land and inland water [46]. The daily temperature data are available four times per day: at approximately 10:30 and 22:30 from Terra and at approximately 1:30 and 13:30 from Aqua (local time).

Other MODIS Data and Landsat Data
The MODIS daily snow product (MOD10A1, version 6) with a spatial resolution of 500 m and daily temporal resolution from 2002 to 2016 was used to compare the temporal pattern of the lake ice cover derived from the MODIS daily temperature product. The MODIS daily snow cover product is available from the U.S. National Snow and Ice Data Center (https://nsidc.org/data/dataprod/mod10.php, accessed on 15 January 2021), and is derived from the SNOWMAP algorithm, which uses the normalized difference snow index (NDSI) and decision strategies to identify snow cover [30,47]. The pixel is labeled as snow when NDSI ≥ 0.4, as recommended by Riggs, G. A. et al. [48]. The land is classified into 6 types: snow, inland water, land, cloud, missing data, and other (unclassified). The assessment showed that MODIS snow cover products have an overall accuracy of~93% under clear-sky conditions [47]. In this study, we used the snow cover data from Terra only, since the cloud/snow discrimination has errors in Aqua cloud-mask algorithm [30,48]. Additionally, the snow cover data from Aqua in version 6 used band 6 to replace band 7 in the NDSI calculation which is different from the previous version. More details about the MODIS snow products can be found in [30,47,48].
The cloud-free images from the MODIS Daily Reflectance product and Landsat images from 2002 to 2016 were also selected to evaluate the lake water/ice cover information from MODIS temperature products. The MODIS Daily Reflectance product (MOD09GQ, version 6), has a spatial resolution of 500 m and daily temporal resolution, and is available at (https://modis.gsfc.nasa.gov/data/dataprod/mod09.php, accessed on 15 January 2021). The Landsat 5, 7, and 8 images have a high spatial resolution of 30 m and a temporal resolution of 16 days, and are available from the standard United States Geological Survey Landsat Surface Reflectance products. Due to the 16-day temporal resolution and the influence of cloud cover, the availability of Landsat images is limited, and cloud-free images were preferentially selected. In our study, 4 images from Landsat under clear sky condition were chosen in each lake during freeze-up and break-up periods, respectively.
The air temperature from 2002 to 2016 was adopted from the China Meteorological Administration (CMA) dataset which has been checked by thorough quality control [49].

Lake Ice Phenology Detection Method
The automated detection of lake ice phenology from MODIS daily temperature data was performed in three main steps ( Figure 2): (i) the MODIS data from Terra and Aqua were combined to increase data availability, (ii) the unfrozen water cover fraction was extracted to characterize the spatial-temporal change pattern of water/ice transition, and (iii) a curve-fitting model was applied to capture the ice phenology characteristics in transitional periods. The individual steps of the method are presented in the following sections.

Terra and Aqua Combination
The lake-extent mask and quality control were applied prior to the Terra and Aqua data combination. The lake mask excluded lake boundaries using a buffer zone (3 × 3pixel matrix [50,51]) to eliminate the lake-land mixed-pixel effect. Quality control was applied to identify "valid pixels" based on the quality flags from the MOD11A1 and MYDA1 quality control layers. When the pixel quality flags were "good quality", "average land surface temperature error ≤ 0.01", "average emissivity error ≤ 0.01", or "average emissivity error ≤ 0.02", the corresponding pixels were identified as "valid pixels". When the pixel quality flags were "cloud effects" or others, the corresponding pixels were identified as "invalid pixels".
The Terra and Aqua data combination (also referred to as the MODIS Daily combi-

Terra and Aqua Combination
The lake-extent mask and quality control were applied prior to the Terra and Aqua data combination. The lake mask excluded lake boundaries using a buffer zone (3 × 3-pixel matrix [50,51]) to eliminate the lake-land mixed-pixel effect. Quality control was applied to identify "valid pixels" based on the quality flags from the MOD11A1 and MYDA1 quality control layers. When the pixel quality flags were "good quality", "average land surface temperature error ≤ 0.01", "average emissivity error ≤ 0.01", or "average emissivity error ≤ 0.02", the corresponding pixels were identified as "valid pixels". When the pixel quality flags were "cloud effects" or others, the corresponding pixels were identified as "invalid pixels".
The Terra and Aqua data combination (also referred to as the MODIS Daily combination) referred to merge MOD11A1 and MYD11A1 data acquired on the same day on a pixel-by-pixel basis. Based on the quality flags, the valid pixels were always given priority over invalid pixels. The output temperature images named as MOD11_Merge and the combination was executed by the following step for a single day: (a) If the same pixel marked as a "valid pixel" was found in more than two images in a given day, the value of the corresponding pixel in MOD11_Merge was calculated by the average temperature from the "valid pixels" and also flagged as "valid"; (b) If the same pixel marked as a "valid pixel" was found in only one image in a given day, the corresponding pixel in MOD11_Merge was taken from the "valid pixel" and also flagged as "valid"; (c) If the same pixel marked as a "valid pixel" was not found in any images in a given day, the corresponding pixel was marked as "invalid".

Extraction of the Unfrozen Water Cover Fraction
The temporal pattern of water/ice transition phase was determined by unfrozen water cover fraction extracted from MOD11_Merge images during the period 2002-2016. The water/ice status in MOD11_Merge image was classified by a temperature threshold initially on a pixel-by-pixel basis. The temperature thresholds were ±0.5 • C, here, and their sensitivity analysis was performed as discussed in Section 5.1. The unfrozen water cover fraction was calculated by the number of pixels classified as "water", accounting for the lake area by the following steps: (a) If the temperature of the pixel was above 0.5 • C, we regarded the pixel as water and classified it as "Pixel water "; (b) If the temperature of the pixel was below −0.5 • C, we regarded the pixel as ice and classified it as "Pixel ice "; (c) If the temperature of the pixel was between −0.5 • C and 0.5 • C, we regarded the pixel as starting to freeze (melt) and classified it as "Pixel ice-water-mixed "; (d) The unfrozen water cover fraction was calculated by Equation (1).
Unfrozen water cover fraction = total number of Pixel water total number of Pixel valid (1) Afterwards, the adjacent temporal filter was used on the same pixel in the previous and subsequent days without reducing spatial and temporal resolution to deduce the continuous stable ice-covered or ice-free pixel from cloud-covered pixels. As the adjacent temporal filter was mainly used for stable-status pixel deduction in our study, two temporal windows of ±5 days and ±7 days were chosen, respectively.
Outliers in the unfrozen water cover fraction were mainly caused by the high percentage of cloud cover, as shown in Figure 3. The outliers were eliminated by the following steps. First, the unfrozen water cover fraction was set to 1 for dates before the mean air temperatures in autumn dropped below 0 • C. The mean air temperature values were preliminarily smoothed by a 31-day running average to remove synoptic variability [52]. Second, the unfrozen water cover fraction was removed if the percentage of "invalid" pixels in an image was larger than 80%. After outlier removal, the temporal pattern of the unfrozen water cover fraction (increase or decrease) was applied to describe the water/ice transition process and then extract lake ice phenology. Remote Sens. 2021, 13, x 7 of Figure 3. Extraction of the unfrozen water cover fraction values in Lake Hulun (a), Lake Ngoring (b), and Lake Qinghai (c). The gray columns represent the cloud cover. The blue dot line represents the unfrozen water cover fraction without outlier removal, and the red dot line represents the unfrozen water cover fraction after outliner removal.

Curve Fitting for Lake Ice Phenology Detection
We proposed a time-dependent logistic function to describe the temporal pattern the unfrozen water cover fraction spanning one hydrological year. The transition betwe open water and ice-covered states was parameterized by the sigmoidal growth cur characterized by an initial slow growth, but accelerating thereafter and decreasing in final phase before approaching the upper asymptote [53]. The growth or decay of the p tial ice cover was reflected by variation of the unfrozen water fraction between the limit values of 1 and 0, where 1 corresponds to open water. The logistic function model h been found to represent the sigmoid growth profile well [37,54] and have frequently be used to fit phenology time series [38,40,41,55]. The logistic function model used here w modified from [37] and is illustrated in Figure 4. The model equation is as follows: where is the unfrozen water cover fraction at time ; is the maximum value of unfrozen water cover fraction, which is equal to 1 in our study; is the rate of curvatu at the inflection point, referred as in the freeze-up process and in the breakprocess; is the days of the hydrological year, is the time at a given inflection po referred to as in the freeze-up process and in the break-up process; and and are the fitting parameters and represent the curvature and position of function, resp tively. The logistic function was similar to the cumulative normal distribution and has t key fitting parameters ( and ) representing the curvature and position of the funct on the time axis, respectively. and Lake Qinghai (c). The gray columns represent the cloud cover. The blue dot line represents the unfrozen water cover fraction without outlier removal, and the red dot line represents the unfrozen water cover fraction after outliner removal.

Curve Fitting for Lake Ice Phenology Detection
We proposed a time-dependent logistic function to describe the temporal pattern of the unfrozen water cover fraction spanning one hydrological year. The transition between open water and ice-covered states was parameterized by the sigmoidal growth curve, characterized by an initial slow growth, but accelerating thereafter and decreasing in the final phase before approaching the upper asymptote [53]. The growth or decay of the partial ice cover was reflected by variation of the unfrozen water fraction between the limiting values of 1 and 0, where 1 corresponds to open water. The logistic function model has been found to represent the sigmoid growth profile well [37,54] and have frequently been used to fit phenology time series [38,40,41,55]. The logistic function model used here was modified from [37] and is illustrated in Figure 4. The model equation is as follows: where y is the unfrozen water cover fraction at time x; a is the maximum value of the unfrozen water cover fraction, which is equal to 1 in our study; k is the rate of curvature at the inflection point, referred as k f in the freeze-up process and k b in the break-up process; x is the days of the hydrological year, x t is the time at a given inflection point, referred to as x f in the freeze-up process and x b in the break-up process; and k and x t are the fitting parameters and represent the curvature and position of function, respectively. The logistic function was similar to the cumulative normal distribution and has two key fitting parameters (k and x t ) representing the curvature and position of the function on the time axis, respectively.
where X is the start or end date of the corresponding transitional period (freeze-up or break-up) and is the open-water fraction threshold, which is set to 0.995 and (1-0.995) for completely ice-free and completely ice-covered lake surfaces, respectively.

Cloud Contamination Removal by Terra and Aqua Combination
The valid pixel fraction calculated by the number of valid pixels in the lake mask, and the valid day fraction calculated by the number of the day that image was covered by 80% valid pixels, were jointly analyzed to evaluate the reduction in cloud contamination ( Table 2). The valid pixel fraction was improved in the MOD11_Merge images to 77.2%, 70.7%, and 79.1% for Lakes Hulun, Ngoring, and Qinghai, respectively. The valid day fraction was improved in the MOD11_Merge images to 66.2% in Lake Hulun, to 54.5% in Lake Ngoring, and 69.1% in Lake Qinghai. These numbers were more than two times larger than those in the corresponding images from MOD11A1 and MYD11A1. Obvious increases were observed in the valid pixel fraction and the valid day fraction in MOD11_Merge images after the combination. Hence, we concluded that the proposed approach effectively increases data availability and increases the amount of good-quality data while retaining the essential information for lake ice extraction. The logistic function was applied by fitting x t and k separately to the freeze-up and break-up periods in each hydrologic year. With k and x t known, the specific ice phenology characteristics, freeze-up start and end dates (FUS and FUE), break-up start and end dates (BUS and BUE), frozen ice cover duration (days between freeze-up start date and break-up end date, FIC), and complete ice cover duration (days between freeze-up end date and break-up start date, CID)) were calculated as follows: where X is the start or end date of the corresponding transitional period (freeze-up or break-up) and y 0 is the open-water fraction threshold, which is set to 0.995 and (1-0.995) for completely ice-free and completely ice-covered lake surfaces, respectively.

Cloud Contamination Removal by Terra and Aqua Combination
The valid pixel fraction calculated by the number of valid pixels in the lake mask, and the valid day fraction calculated by the number of the day that image was covered by 80% valid pixels, were jointly analyzed to evaluate the reduction in cloud contamination ( Table 2). The valid pixel fraction was improved in the MOD11_Merge images to 77.2%, 70.7%, and 79.1% for Lakes Hulun, Ngoring, and Qinghai, respectively. The valid day fraction was improved in the MOD11_Merge images to 66.2% in Lake Hulun, to 54.5% in Lake Ngoring, and 69.1% in Lake Qinghai. These numbers were more than two times larger than those in the corresponding images from MOD11A1 and MYD11A1. Obvious increases were observed in the valid pixel fraction and the valid day fraction in MOD11_Merge images after the combination. Hence, we concluded that the proposed approach effectively increases data availability and increases the amount of good-quality data while retaining the essential information for lake ice extraction.

Comparison of the Unfrozen Water Cover Fraction with Other Datasets
The performance of the unfrozen water cover fraction obtained from MOD11_Merge daily images was compared with that from MODIS daily snow products, MODIS daily reflectance products, and Landsat images. These products have previously proven their ability to extract lake ice characteristics [18,21,23,56].
The unfrozen water cover fraction calculation and the outlier removal method (Section 3.2) were applied to the MODIS daily snow cover products. As shown in Figures 5 and 6, the unfrozen water cover fraction obtained from MOD11_Merge was overall consistent with that of the MODIS snow cover product, with values of R 2 larger than 0.70 during the freeze-up pattern and 0.89 during the break-up pattern (Table 3). Compared with the freeze-up patterns, the break-up patterns in MOD11_Merge and MODIS snow cover product showed higher agreement, with R 2 values of 0.912, 0.889, and 0.948, respectively. The freeze-up pattern in MOD11_Merge and MODIS snow cover in Lake Qinghai had a lower R 2 , bias, and RMSE than the other two lakes, with values of 0.72, 0.22, and −0.55, respectively. This may be due to the thin ice or initial ice information, which was not easily detectable, especially in large lakes [14,15].   Figure 6. Comparison of unfrozen water cover fractions in Lake Hulun, Lake Ngoring, and Lake Qinghai from MOD11_Merge temperature images and MODIS snow product data, taking the year 2011 as an example. The visual inspections of lake ice under clear-day conditions from the MODIS daily reflectance products and Landsat images were further used to evaluate the water/ice cover information. The unfrozen water cover fraction in the phases of break-up process was consistent with MODIS reflectance products, MODIS snow products, and MOD11_Merge temperature images, as shown in the Supplementary Materials, Figures S2-S4. This consistency was notably seen in the developed and final phases of ice cover than in the initial phase, especially in smaller Lake Ngoring (Supplementary Materials, Figure S4). When compared with Landsat images with a high resolution of 30 m, thin ice or fast ice at the lakeshores was difficult to detect from the MODIS products, either from snow products and temperature products, as shown in the Supplementary Materials, Figures S5-S7. The misclassification of thin ice was possibly caused by the lower temperature resulting from undetected cloud cover in the MODIS temperature products [23].

Comparison of Derived Lake Ice Phenology with other Lake Ice Datasets
The lake ice phenology was extracted from the temporal pattern of unfrozen water cover fraction by characterizing the water/ice transition process (as shown in Supplementary, Figure S1). The extracted ice phenology dates of Lake Qinghai, Lake Ngoring, and Lake Hulun were then compared with published lake ice datasets including microwave and optical satellite sensors, as listed in Table 4. As shown in Figure 7a-d, the annual variation of ice phenology in Lake Qinghai from the proposed method was consistent with other datasets, with an R 2 of 0.90 when compared with MODIS snow products [56], an R 2 of 0.87 when compared with MODIS reflectance data [19], and with an R 2 of 0.80 when compared with passive microwave data [15]. Table 4. The data sources and the durations of other publishing lake ice phenology datasets.

Lake Name Data Source and Reference Duration
Lake Hulun MODIS reflectance product [43] 2002-2016 Lake Ngoring MODIS snow product [23] 2002-2015 AMSR-E/2 [57] 2000-2016 Lake Qinghai SSM/I [15] 2002-2015 MODIS snow product [23] 1979-2016 AMSR-E/2 [57] 2000-2016 MODIS reflectance product [19] 2000-2016 Remote Sens. 2021, 13, x Figure 7. Comparison of lake ice phenology dates between our study (the right y-axis) and other studies (the l in Lake Qinghai (a-d), Lake Ngoring (e-f), and Lake Hulun (g-h). The lake ice phenology dates extracted in our shown as dotted lines. The y-axis is the day of the hydrological year (1 September is referred to as Day 1 in our  Comparison of lake ice phenology dates between our study (the right y-axis) and other studies (the left y-axis) in Lake Qinghai (a-d), Lake Ngoring (e-f), and Lake Hulun (g-h). The lake ice phenology dates extracted in our study are shown as dotted lines. The y-axis is the day of the hydrological year (1 September is referred to as Day 1 in our study).
The trends of break-up dates obtained from MOD11_Merge data had a high consistency with those of other MODIS products, less with passive microwave data, especially in Lake Qinghai, as shown in Table 5. However, the trends of freeze-up dates showed inconsistencies. For example, the trend of FUS in Lake Qinghai was −4.09 days/10 yr, while it was 4.00 days/10 yr from Cai, Y. et al. [15] and 6.31 days/10 yr in our study. This inconsistency could potentially be attributed to several factors, such as misclassification of the thin ice at the initial phase of freeze-up period, inhomogeneous definitions of lake ice phenology in different studies, and the coarser spatial resolution of the passive microwave data [13,14].

Interannual Variability in Lake Ice Phenology
The extraction approach was applied to other lakes in different climatic regions of northern China. The temporal patterns and trends of ice phenology from 2002 to 2016 are shown in Figure 8 and Table 6. In general, shallower lakes had earlier freeze-up dates than deeper lakes (Figure 8): FUS and FUE in Lake Hulun (mean depth 5 m) were 69 and 91 days earlier than in Lake Namco (mean depth 43 m). The trends in ice phenology variables among the Tibetan lakes during the period 2002-2016 revealed a tendency towards delayed FUS dates, with the highest value being 6.31 days/10 yr (p < 0.05) in Lake Qinghai and the lowest value being 0.09 days/10 yr in Lake Zharinamco. Trends of BUE dates were also positive in all lakes on the Tibetan Plateau except Lake Silingco and Lake Qinghai, with insignificant BUE trends of −1.00 days/10 yr and −4.69 days/10 yr, respectively. Compared to the lakes on the Tibetan Plateau, the break-up starts and end dates rapidly shifted to be earlier in Lake Hulun (−3.73 days/10 yr; −5.02 days/10 yr). The short length of the data record (15 years) may be one of reasons that trends have no significance. layed FUS dates, with the highest value being 6.31 days/10 yr (p < 0.05) in Lake Qinghai and the lowest value being 0.09 days/10 yr in Lake Zharinamco. Trends of BUE dates were also positive in all lakes on the Tibetan Plateau except Lake Silingco and Lake Qinghai, with insignificant BUE trends of −1.00 days/10 yr and −4.69 days/10 yr, respectively. Compared to the lakes on the Tibetan Plateau, the break-up starts and end dates rapidly shifted to be earlier in Lake Hulun (−3.73 days/10 yr; −5.02 days/10 yr). The short length of the data record (15 years) may be one of reasons that trends have no significance.

Factors Influencing Lake Ice Phenology
The growth and decay of lake ice cover are governed by thermodynamic interactions between atmosphere-water-ice interfaces and are affected by climatic and lake-specific factors. The major climatic factors controlling the heat exchange with the atmosphere are air temperature, solar radiation, precipitation, and wind speed [5,8,10]. Lake-specific factors, such as lake morphometry (depth, area, and volume) and location (longitude, latitude, and altitude), determine heat storage in lakes and additionally affect lake ice formation [58]. To assess how lake-specific and climatic factors contribute to the freeze-up and break-up patterns and the changes in the ice cover duration, the correlations between ice phenology time series extracted in our study and possible influencing factors including lake specific factors and major climatic factors were analyzed. The climatic variables, including air temperature, wind speed, solar radiation, and snow depth, were obtained from 2002 to 2016 from the China Meteorological Administration dataset based on the nearest weather station. The mean values of the lake surface water temperature (LSWT) data during the ice-off period used in our study were calculated based on the MOD11_merge dataset. Our LSWT results showed high consistency with publicly available satellite-derived LSWT products, such as the Arclake [59] and LSWT datasets from AVHRR [60], with R 2 values of 0.91 and 0.93, respectively. The other climatic factors were obtained from weather stations near the lakes.
As shown in Figure 9, the results suggest that climate factors strongly affected the break-up pattern, while lake-specific factors had a greater effect on the freeze-up pattern that agrees with previously published findings [61][62][63]. Deep high-altitude lakes showed later freeze-up dates and shorter ice cover durations, while shallower lakes located at high latitudes showed earlier freeze-ups and longer ice cover durations. The air temperature had a significant effect on the break-up date and the duration of the ice-covered period. The dates at which the air temperature reached 0 • C in the spring even had a high R 2 value (0.81) for the break-up periods. The LSWT had a more obvious impact on the freeze-up date (freeze-up day and rate) and CID than on the break-up pattern.
Remote Sens. 2021, 13, x 15 of 23 Figure 9. Correlation maps between lake ice phenology variables and climatic (a) and lake-specific factors (b). Red indicates a positive correlation, and blue indicates a negative correlation. The size of the circle represents the absolute value of the correlation coefficient. The * represents significance at a confidence level of 0.1. Ts (ice_off) is the lake water surface temperature during the ice-off period, and Ta (yearly) is the annual mean air temperature. The 0 °C dates (autumn/spring) represent the dates at when air temperature reached 0 °C in the autumn/spring. WD (yearly) and SWR (yearly) represent the annual mean wind speed and shortwave solar radiation, respectively.

Sensitivity Analysis of Classification for Water/ice Status Pixels
The sensitivity of the water/ice classification status to the threshold value was evaluated by a lake ice phenology extraction under different freezing and melting temperatures, which was varied in 0.25 °C intervals within the range from −1.25 °C to 1 °C. Figure  10 shows the temporal pattern of the ice fitting variables ( and , and and ) among the three lakes, under different freezing and melting temperature settings. Con- Figure 9. Correlation maps between lake ice phenology variables and climatic (a) and lake-specific factors (b). Red indicates a positive correlation, and blue indicates a negative correlation. The size of the circle represents the absolute value of the correlation coefficient. The * represents significance at a confidence level of 0.1. T s (ice_off) is the lake water surface temperature during the ice-off period, and T a (yearly) is the annual mean air temperature. The 0 • C dates (autumn/spring) represent the dates at when air temperature reached 0 • C in the autumn/spring. WD (yearly) and SWR (yearly) represent the annual mean wind speed and shortwave solar radiation, respectively.

Sensitivity Analysis of Classification for Water/ice Status Pixels
The sensitivity of the water/ice classification status to the threshold value was evaluated by a lake ice phenology extraction under different freezing and melting temperatures, which was varied in 0.25 • C intervals within the range from −1.25 • C to 1 • C. Figure 10 shows the temporal pattern of the ice fitting variables (x f and k f , and x b and k b ) among the three lakes, under different freezing and melting temperature settings. Consistent fluctuations were found in the interannual variabilities in x f and x b . The variabilities in k f and k b were relatively stable except when the temperature threshold was too high or too low ( Figure 11). As shown in Figure 10a-c, the freeze-up process occurred earlier when the freezing temperature was set higher. The temporal pattern of the unfrozen water cover fraction was stable and consistent when the freezing temperature was set between −0.75 • C and 0.5 • C. The break-up process showed relatively high variability in the initial phase, especially when the melting temperature was set lower than 0 • C (Figure 10d-f). This variability was possibly due to the melt-refreeze events captured when the temperature setting was too low, which were not representative of the starting point of the actual break-up process.
The trends of x f and k f were consistent when the temperature threshold varied between −0.75 • C and 0.5 • C (Figure 12a,c). The trends of x b performed consistently well under different temperature settings, while the trends of k b showed good consistency when the melting temperature was above 0 • C, possibly caused by the high fluctuations in the initial phase of the break-up process. These findings suggest that the trends of ice phenology are robust when the freezing temperature is set between −0.75 • C and 0.5 • C and the melting temperature is set above 0 • C.
Remote Sens. 2021, 13, x 16 of 2 in the initial phase of the break-up process. These findings suggest that the trends of ic phenology are robust when the freezing temperature is set between −0.75 °C and 0.5 °C and the melting temperature is set above 0 °C. The temperature classification applied for the water/ice status should consider th lake surface situation in the freeze-up pattern and break-up pattern, respectively. Th freezing temperature can be influenced by several factors. Water salinity can decrease th freezing temperature and delay ice formation [64]. Moreover, the ice pixel temperatur can be influenced by the bubble contents, impurities, and fracture patterns of ice surface [65]. Ice break-up starts after the ice/snow surface is warmed up to the melting point [63 Melting snow, snow-ice and black ice have different reflectance values [66,67]. Further more, the day and night melt-freeze events produce meltwater on the ice surface befor the ice cover break-up completely. To consider these uncertainties, and the physical basi of ice cover, we determined the classified pixel as "ice" when its temperature was lowe than −0.5 °C and classified the pixel as "water" when its temperature was higher tha 0.5°C.   (d−f), and the third row presents the time series of (g−i), and the fourth row presents the time series of (j-l). The temperature classification applied for the water/ice status should consider the lake surface situation in the freeze-up pattern and break-up pattern, respectively. The freezing temperature can be influenced by several factors. Water salinity can decrease the freezing temperature and delay ice formation [64]. Moreover, the ice pixel temperature can be influenced by the bubble contents, impurities, and fracture patterns of ice surfaces [65]. Ice break-up starts after the ice/snow surface is warmed up to the melting point [63]. Melting snow, snow-ice and black ice have different reflectance values [66,67]. Furthermore, the day and night melt-freeze events produce meltwater on the ice surface before the ice cover break-up completely. To consider these uncertainties, and the physical basis of ice cover, we determined the classified pixel as "ice" when its temperature was lower than −0.5 • C and classified the pixel as "water" when its temperature was higher than 0.5 • C. Figure 12. (a-d) The ice phenology variable trends under different freezing and melting temperatures of pixel settings. The x-axis is the lake, and the y-axis is the trend of the unfrozen water cover fraction. The red and blue tiles indicate positive trends and negative trends, respectively.

The Advantages and Limitations of the Method
The proposed algorithm of the ice phenology extraction possesses several distinct advantages for wide-scale application on seasonal ice-covered lakes. An important qualitative advance is to scale the transition period between the open water and the ice-covered stage by a parameterized function of time (Equation (3)). The scaling-based method has been proven to be robust to outliers and can be applied to lakes with varying morphometry. Additionally, the characterization of the seasonal ice phenology can be extended beyond the binary ice-on/ice-off events by the direct assessment of the "partial ice cover" duration. The latter characteristic is crucial for understanding the lake response to ice cover formation, especially in medium to large lakes. During ice formation in autumn, the length of the transitional period determines the rate of cooling across the entire water column: a faster formation of the complete ice cover allows more heat to be stored in the near-bottom part of the lake, affecting the water-ice heat fluxes throughout the entire icecovered period and the ice thickness [63,68,69]. In turn, the length of the transition period to the ice-free state in spring causes long-lasting effects on the vertical thermal stratification in the following summer and thereby affects crucial lake water quality characteristics, such as the dissolved oxygen content [63,69]. Both autumn and spring transitional periods are the results of interaction between the air-lake heat transport, the available solar radiation, and the heat storage in the lake water column. Therefore, they serve as sensitive indicators of the regional response to global climate change. When applied to long-term Figure 12. (a-d) The ice phenology variable trends under different freezing and melting temperatures of pixel settings. The x-axis is the lake, and the y-axis is the trend of the unfrozen water cover fraction. The red and blue tiles indicate positive trends and negative trends, respectively.

The Advantages and Limitations of the Method
The proposed algorithm of the ice phenology extraction possesses several distinct advantages for wide-scale application on seasonal ice-covered lakes. An important qualitative advance is to scale the transition period between the open water and the ice-covered stage by a parameterized function of time (Equation (3)). The scaling-based method has been proven to be robust to outliers and can be applied to lakes with varying morphometry. Additionally, the characterization of the seasonal ice phenology can be extended beyond the binary ice-on/ice-off events by the direct assessment of the "partial ice cover" duration. The latter characteristic is crucial for understanding the lake response to ice cover formation, especially in medium to large lakes. During ice formation in autumn, the length of the transitional period determines the rate of cooling across the entire water column: a faster formation of the complete ice cover allows more heat to be stored in the near-bottom part of the lake, affecting the water-ice heat fluxes throughout the entire ice-covered period and the ice thickness [63,68,69]. In turn, the length of the transition period to the ice-free state in spring causes long-lasting effects on the vertical thermal stratification in the following summer and thereby affects crucial lake water quality characteristics, such as the dissolved oxygen content [63,69]. Both autumn and spring transitional periods are the results of interaction between the air-lake heat transport, the available solar radiation, and the heat storage in the lake water column. Therefore, they serve as sensitive indicators of the regional response to global climate change. When applied to long-term remote sensing observations, the proposed approach would allow tracing the climatically driven effects on non-linear interactions between various components of the surface heat balance during autumn and spring. The proposed approach also builds on the strengths of the MODIS daily temperature products. The Terra and Aqua daily combination from four observations with a time step of 4 to 7 hours demonstrates the effectiveness in reducing the cloud contamination while retaining the essential spatial information. The automatic temperature-based algorithm was also shown the ability to correctly extract ice phenology on different lake types concerning potential influences caused by salinity, morphometry, and geographic location.
The limitations of the proposed algorithm are intrinsic to remote sensing products. Firstly, thin ice cover is formed when lake surface temperature approaches the freezing temperature, but it may be repeatedly destroyed during windy periods with prolonged fluctuation of the surface temperature around 0 • C [63,69]. Therefore, thin ice is difficult to detect using temperature-based methods. Secondly, the severe cloud contaminations caused inevitable data gaps, although the data availability was improved by Terra and Aqua daily combination in our study. Additionally, MODIS-derived LSWT was revealed to have a negative bias compared with in situ observations, especially during the daytime. The negative bias could potentially have been caused by several things, including undetected cloud cover [46,51], cool skin, and warm layer effects [70][71][72]. Furthermore, Terra and Aqua daily combination could cause additional uncertainties due to ignoring the diurnal variation of lake surface temperature. Finally, although the spatial resolution of 1 km from MODIS temperature products is relatively high compared to other remote sensing products, like AMSR-E, it is still relatively coarse for phenology studies on lakes with horizontal dimensions smaller than several kilometers, especially located at high elevations or in areas covered by persistent snow cover. Progress in phenology observations on small lakes could be made by obtaining sub-kilometer spatial resolutions while retaining the sub-daily temporal resolution of data.

Conclusions
In this study, we developed a new automatic method to characterize lake ice phenology using MODIS daily temperature products from 2002 to 2016 by parameterizing the temporal pattern of the water/ice transition process. Based on this method, the ice phenology of nine seasonal ice-covered lakes in different climate regions was evaluated. The results showed that Lake Qinghai had the most obvious trend towards delayed freeze-up start date (FUS) (6.31 days/10 yr) among the lakes of the Tibetan plateau, and Lake Hulun the most obvious tendency towards earlier break-up start and end dates (BUS and BUE), with change rates of -3.73 days/10 yr and −5.02 days/10 yr, respectively. The proposed method has two major advantages, making it robust and applicable: (i) it avoids any lake-specific or empirical thresholds, allowing the effective treatment of data gaps caused by atmospheric conditions, and (ii) it uses a time-dependent parameterized function to characterize the water/ice transition process, making it possible to trace the ice phenology on a more detailed level than traditional binary classification. The algorithm could provide data support for the correct interpretation of the global-scale climate change effects on lake ice phenology: apart from the general shortening of the ice-cover period, the interplay between atmospheric warming and seasonally available shortwave solar radiation may produce significant changes in the duration of transitional-partially ice-covered-periods with important consequences for air-lake interactions and internal lake mixing.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs13142711/s1. Figure S1. The process of curve fitting of unfroze water fraction in three lakes Hulun from 2002 to 2016. The performances of logistic function fitness of unfrozen water cover fraction are good with the R 2 large than 0.95 in three lakes. Figure S2. Comparison of unfrozen water fraction extraction in 2011 in Lake Hulun from MODIS reflectance image, MODIS snow product, and MOD11_Merge data. The images in the first column are from MOD09 daily reflectance product (RGB:341). Figure S3. The comparison of unfroze water fraction extraction in 2011 in Lake Qinghai from MODIS reflectance image, MODIS snow product, and MOD11_Merge data. The images in first column were from MOD09 daily reflectance production (RGB:341). Figure S4. The comparison of unfroze water fraction extraction in 2011 in Lake Ngoring from MODIS reflectance image, MODIS snow product, and MOD11_Merge data. The images in first column were from MOD09 daily reflectance production (RGB:341). Figure S5. The comparison of unfroze water fraction extraction in Lake Hulun from MOD11_Merge temperature data and Landsat images. The images in first column were the classification of water and ice based on the proposed method in this day. The images in second column represent the lake surface situation from Landsat images (RGB true color combination: band 532) at the same dates. The third column represent the green cycle area zoomed in Landsat images. The green cycles represent the difference between MOD11_Merge and Landsat images. Figure S6. The comparison of unfroze water fraction extraction in Lake Hulun from MOD11_Merge temperature data and Landsat images. The images in first column were the classification of water and ice based on the proposed method in this day. The images in second column represent the lake surface situation from Landsat images (RGB true color combination: band 532) at the same dates. The third column represent the green cycle area zoomed in Landsat images. The green cycles represent the difference between MOD11_Merge and Landsat images. Figure S7. The comparison of unfroze water fraction extraction in Lake Hulun from MOD11_Merge temperature data and Landsat images. The images in first column were the classification of water and ice based on the proposed method in this day. The images in second column represent the lake surface situation from Landsat images (RGB true color combination: band 532) at the same dates. The third column represent the green cycle area zoomed in Landsat images. The green cycles represent the difference between MOD11_Merge and Landsat images.
Funding: This study was funded by the National Key R&D Program of China (2017YFA0603601) and the National Natural Science Foundation of China (41930970). The study was performed within the project "Lakes of Tibet in the Climate System" (LaTiCS) funded by the Sino-German Centre for Research Support (CDZ Project GZ1259) and German Research Foundation (DFG Project KI 853/13-1). GK was additionally supported by the DFG Project KI 853/16-1.