Analysis of Ice Phenology of Middle and Large Lakes on the Tibetan Plateau

Considered as a sensitive indicator of climate change, lake ice phenology can have significant influences on regional climate by affecting lake-atmosphere energy and water exchange. However, in situ measurements of ice phenology events are quite limited over high-elevation lakes on the Tibetan Plateau, where satellite monitoring can make up such deficiency. In this study, by a combination of AMSR-E (2002–2011) and AMSR-2 (2012–2021) passive microwave data, MODIS optimal products and in situ measurements of temperature profiles in four lakes, the ice phenology events of 40 high-elevation large lakes were derived and their inter-annual trends and influencing factors were analyzed. The freeze-up start date (FUS) mainly occurs in November-December with an average date of 9 December and the break-up end date (BUE) is concentrated in April-May with a multi-year average of 5 May. Under climate warming, 24 of the 34 (70.6%) lakes show delayed FUS at an average trend of 0.35 days/year, and 7 (20.6%) lakes show advanced BUE (rate of change CR = −0.17 days/year). The average ice coverage duration (ID) was 147 days, and 13 (38.2%) lakes shortened ID at an average rate of −0.33 days/year. By synthesizing other ice phenology products, we obtained the assembled products of lake ice phenology, and found that air temperature dominates during the freeze-thaw process, with a higher dependence of BUE than that of FUS on air temperature.


Introduction
Lake ice, as an important component of the cryosphere, has important implications for local and global climate, energy balances, and hydrological cycles. For example, longer ice-off conditions increase the heat transfer, the annual evaporation, and the greenhouse gases emissions from lakes due to chemical and biological processes [1]. In addition, the cooling effect of latent heat flux through evaporation can affects regional climate, such as lake-effect clouds and precipitation [2,3]. At the same time, lake ice is sensitive to changes in climate [4]. Some studies have suggested that lake ice phenology is a good proxy for climate change, and both freeze-up and break-up dates can represent air temperature change [5][6][7]. Moreover, in some case, lake ice records may be a more reliable climate indicators than air temperature records [8,9]. However, it is very difficult to conduct high-density and continuous artificial lake ice observations. The development of remote sensing technology has made it possible to monitor lake ice continuously on a regional and global scale [1,10,11] surface and the atmosphere, and then cause the change of the regional or even the whole TP heating [39,40]. Studies have shown that the thermal effect of the TP is an important reason for its influence on the atmospheric circulation and the climate of Asia and even the world [41][42][43][44]. At the same time, the TP was identified as one of the most sensitive regions in the world to changes in climate [45]. However, due to the lack of ground observation, studies on the impact of global climate change on the TP is still limited. Lake ice phenology retrieved from remote sensing data, as an effective indicator of climate change, can fill a gap in our knowledge about the impact of climate change in this remote region. Many lake ice phenology datasets on the TP have been developed [18,[30][31][32][46][47][48][49][50][51][52]; however, the uncertainties in multiple remote sensing data, multiple methods, and multiple definitions of lake ice phenology parameters in previous studies, make the consistency between lake ice phenology datasets very weak. In this study, based on the passive microwave remote sensing data and in situ measurements, we summarize and analyze the lake ice phenology of the Tibetan Plateau. The contents are as follows: Section 2 introduces the study area and methods, including the introduction of data, and the extraction method of ice phenology events, etc.; Section 3 introduces the integrated analysis of eight lake ice phenology datasets, the influencing factor analysis, the cluster analysis, and the trends of lake ice phenology events during 2002-2021; Section 4 discusses the influences of snow cover on the lake ice phenology events and their trends, the effect of pixel position, and the source of errors existing in this study; Section 5 documents the conclusions.

Study Area and Methods
Tibetan Plateau (26 • 00 12 N-39 • 46 50 N, 73 • 18 52 E-104 • 46 59 E), known as "the Third Pole of the earth", is located in the southwest of China, and has an area of 2572.4 × 10 3 km 2 [53]. The Tibetan Plateau is rich in water resources, and lakes with an area of 1 km 2 and above account for 57.2% of the total lake area in China [54]. In order to adapt to the spatial resolution of passive microwave remote sensing data, this study selected 40 medium and large lakes for ice phenology research ( Figure 1, Table 1). The ASTER Global Digital Elevation Model V003, which has a horizontal resolution of 30 m and can be obtained from NASA EARTHDATA (https://www.earthdata.nasa.gov/, accessed on 8 January 2022), was used to obtain the attribute information (area, altitude, latitude, and longitude) of the targeted lakes (Table 1). Table 1. Lake attributes and average lake ice phenology parameters. The average lake ice phenology parameters are the integration results of eight datasets (this study and seven datasets [18,[30][31][32]47,50,51]). Acronyms denote start of ice break-up (BUS), end of ice break-up (BUE), start of freeze-up (FUS), end of freeze-up (FUE) dates, and ice coverage duration (ID).    Table 1.   are 1000 m [55]. This study selected MODIS 1B products (MOD02HKM and MYD02HKM) (https://www.earthdata.nasa.gov/, accessed on 25 March 2022) with a resolution of 500 m after radiometric calibration, and the lake ice was visually discriminated according to a false-color composite image formed by the red band (B1), the near-infrared band (B2), and the blue band (B3).
The passive microwave imager AMSR-E (The Advanced Microwave Scanning Radiometer for EOS) and AMSR-2 (The Advanced Microwave Scanning Radiometer 2) were launched with the Aqua satellite and the Global Change Observation Mission 1st-Water satellite respectively, which provide vertical and horizontally polarized brightness temperature observations with spatial resolutions ranging from 3 × 5 km (89.0 GHz) to 75 × 43 km (6.925 GHz). The G-Portal (https://gportal.jaxa.jp/gpr/, accessed on 4 January 2022) provides Level-3 brightness temperature products of AMSR-E and AMSR-2. In this study, the level-3 horizontal polarization brightness temperature data of 18 GHz was selected. Daily snow thickness distribution data in China stretch from 1 January 1979 to 31 December 2020, with a spatial resolution of 25 km [56][57][58][59]. This dataset was used to calculate the average snow days with a snow depth of ≥5 cm from December to May in five regions of the TP from 2002 to 2020 (Figure 1), focusing on the impacts of snow cover on the break-up end date.

ERA5-Land
ERA5-Land is a reanalysis dataset provided by the European Centre for Medium-Range Weather Forecasts. It provides data on water and energy cycles from 2.89 m below ground to 2 m above the surface since 1950, with a horizontal resolution of 0.1 • (9 km) and a temporal resolution up to 1 h. This study uses the ERA5-Land monthly average data provided by the Copernicus Climate Change Service Data Platform from 1950 to the present (https://cds.climate.copernicus.eu/, accessed on 1 April 2022), and selects the monthly average data of 2 m air temperature, 10 m wind speed, surface downward shortwave radiation and snow depth from 2000 to 2021 to analyze the influencing factors of ice phenology events.

Lake Temperature Observation
In Nam Co [60][61][62], Peiku Co [63], Bangong Co, and Daze Co [64], the vertical temperature gradient measurements have been collected. The observation period of Nam Co is from November 2011 to June 2014 (daily), the observation period of Peiku Co is from June 2016 to May 2017 (hourly), and the observation period of Bangong Co is July 2012 to August 2013 (hourly), the observation period of Daze Co is from August 2012 to August 2013 (hourly). The lake temperature observations can provide in situ datasets for judging the freezing and thawing dates of lakes (as shown in Supplementary Figure S1): take the first minimum point of water temperature lower than the maximum density temperature as the freeze-up start date, the minimum water temperature point as the freeze-up end date, the point where the water temperature starts to rise rapidly and the stratification is obvious as the break-up start date, and the point where the water temperature first peaked after the break-up start date is the break-up end date.

Extraction of Lake Boundaries
After the DEM was projected, the continuous area with slope 0 was identified as a lake to determine the lake boundary. The lake boundaries of HydroLAKES [65] were used to verify the extraction results, and lakes with an area >200 km 2 were selected as the target lakes ( Figure 1 and Table 1).

Extraction of Lake Ice Phenology Events
Phenology is the field of science concerned with cyclic and seasonal natural phenomena, especially in relation to plant and animal life and climate. In glaciology, freshwater ice scientists refer to freeze-up/break-up and duration of ice on lakes and rivers as ice phenology [66]. Lake ice freeze-up start (FUS) corresponds to the date when detectable ice appears. Freeze-up end (FUE) is defined as the date when the lake is fully ice-covered. Break-up start (BUS) is the first day of the year on which a detectable ice-free water surface appears. Break-up end (BUE) is defined as the date on which a lake is completely clear of ice. Time between FUS and FUE determines the freeze-up duration (FUD), while time between BUS and BUE determines the break-up duration (BUD). Lake ice cover duration (ID) is defined as the period between FUS and BUE, and complete ice cover duration (CID) is defined as the period between FUE and BUS.
In this paper, the lake ice phenological parameters were determined by the brightness temperature curve of the central lake pixel. The linear interpolation method was used to complete the missing brightness temperature data, and the median filtering method was used to denoise the brightness temperature sequence. The brightness temperature difference search method [24] combined with temperature threshold was used to automatically extract FUE and BUS dates. The specific discrimination algorithm is shown in Equations (1) and (2). (2) where, TB i is the filtered brightness temperature on the ith day, i = 1, 2, . . . , 365 or 366, and min () and max () are the minimum and maximum functions. The maximum and minimum values of the time series of the brightness temperature difference were taken as the abrupt change points of the brightness temperature curve.
The positive and negative curve intersection method [25] combined with the brightness temperature threshold was used to determine FUS and BUE. The positive and negative curves are shown in Equations (3)- (6). Curves Y1 and Y4 were used to determine FUS, and curves Y2 and Y3 were used to determine BUE; the method is described in more detail in the Supplementary Materials.
After lake ice phenology parameters were automatically extracted, the visual interpretation of the brightness temperature curve was used to adjust the unreliable results.

Statistical Analysis
The Mann-Kendall trend analysis [67,68]  (negative) statistic Z is used to represent the increasing(decreasing) trend, and the size of |Z| represents the significance of the trend change. If |Z| ≥ 1.28, the trend characteristics are significant, and the confidence level is ≥ 90%. The trend values are characterized by the slope β. In this paper, the statistic Z was used to analyze the trend change, and the slope β was used to analyze the rate of change.
In the integrated analysis of the eight lake ice phenology datasets, the deviation of each dataset was measured by the root mean square error (RMSE, Equation (7)), where the subscript i is the ith lake (i = 1,2, . . . , 40), n is the total number of lakes (n = 40), j is the jth lake ice phenology dataset (j = 1,2, . . . ,8), y ij is the ice phenology dates of the ith lake in the jth dataset, and y i is the average ice phenology of the ith lake calculated from multiple datasets. RMSE j is RMSE of the jth dataset.
In order to find the dependence of the lake average freezing and thawing dates (Table 1) on climatic factors (temperature, wind speed, snow depth, radiation) and local factors (latitude and longitude, area), correlation analysis was carried out, and a t-test was used to test the significance of the correlation coefficient (r). Finally, cluster analysis on the lakes was conducted based on the most important factors.

Verification and Assessment of Lake Ice Phenology Datasets
The extracted lake ice phenology events indicate that some lakes were not frozen/ completely frozen during 2002-2021 (as shown in Supplementary Table S1). To illustrate the reliability of the results, MODIS false-color composite images during ice covered seasons were selected to verify the unfrozen conditions. The MODIS false-color image show consistency with the results obtained from the AMSR-E/2 data (as shown in Supplementary Figure Figure S2f). Xijir Ulan Lake is the only extra, where the lake was defined as not completely frozen by the MODIS composite image, but as not frozen by AMSR-E/2. This may be probably related with the position of the selected microwave pixel, check Section 4.2 for more content.
The estimated ice phenology events from the satellite datasets and the lake temperature profiles (as shown in Supplementary Figure S1) are summarized in Table 2; the FUE and BUE dates obtained by satellite products are close to the results obtained by the lake temperature profiles, with high consistency regardless of the size of the lake. In contrast, the FUS and BUS dates obtained by the lake temperature profile are earlier than the satellite products, where ice forming or melt can happen at shallower depths, but satellite products can only sense the status of the lake surface at deeper depths because of its coarse resolutions. For example, from 2012 to 2013 in Dagze Co, the biases of FUE and BUE dates in most of studies were within 1 and 2 days respectively, while that of FUS and BUS dates ranged from 16 to 25 days and 37 to 47 days, respectively. This similar phenomenon with relatively small (large) biases in FUE and BUE (FUS and BUS) dates also occurred in Bangong Co and Nam Co. The observation data of Paiku Co from 2016 to 2017 show that the lake was in a relatively uniform state of mixing throughout the winter and no frozen ice appeared (as shown in Supplementary Figure S1). The MODIS image on 19 Mar 2017 proves this situation (as shown in Supplementary Figure S2c); however, the dates of freezing and thawing were determined by Guo et al. [51], who used model simulations.
Further, the multi-source average of each lake ice phenology parameter was used as the standard, and the root mean square error (RMSE) was used to evaluate the deviation of each dataset. Overall, the results from Qiu [47] and Cai et al. [50] are closer to the average lake ice phenology (RMSE are both within 6 days). Although this study was not the best when compared with the average results, our estimation shows some advantages when compared with that of lake temperature profiles ( Table 2). For example, the estimated FUS, FUE (BUE) dates in Nam Co during 2012-2013 (2013-2014) in this study were the closest to the estimation from the lake temperature profiles. The errors of FUE and BUE dates between our estimation and those estimated through lake temperature profiles were all within 5 days. These indicate that the lake ice phenology parameters determined in this study are reliable.

Comprehensive Integrated Analysis of Lake Ice Phenology Datasets
The derived lake ice phenology events by our method were integrated with the results of seven other studies [18,[30][31][32]47,50,51] (Figure 2a-d), and the average freeze-thaw dates of the 40 lakes in the TP region were obtained ( Table 1).
The boxplot shows the magnitude of the difference between the datasets. Through the different remote sensing data types, the FUS determined by passive microwave is generally later than the results obtained by MODIS (Figure 2a), and this can be attributed to the freezing and thawing patterns of lakes. Many lakes begin to freeze in the shallow water near the lakeshore. The MODIS pixels with higher spatial resolution cover the entire lake surface, and this change is more easily observed, while passive microwave pixels, with a coarse resolution and locating at the center over deep water, have less interference from the surrounding land surface. Differences in the spatial distribution of pixels lead to differences in the observed timing of lake ice formation, which is similar to the results of Cai et al. [19,50]. In addition, the initially formed thin ice may not be easily detected by the brightness temperature threshold [11], resulting in a later FUS determined by passive microwave data. Meanwhile, most of the FUE (BUS) dates obtained by passive microwave are earlier (later) than MODIS (Figure 2b,c), which may be since the lake ice ratio corresponding to the brightness temperature threshold is smaller than that used by MODIS to determine lake ice phenology parameters. Moreover, compared with the FUS, FUE, and BUS dates, the consistency of the BUE date is better, and there is no particularly large difference (Figure 2d), which is consistent with the conjecture of Cai et al. that the impact of different remote sensing data on the freezing date may be greater than that on melting dates [50]. However, in general, the BUE date extracted by passive microwave is Sensors 2023, 23, 1661 9 of 20 earlier than that of MODIS data, which may be related to the earlier ice melting in deep water. The calculation of duration also reflects the differences in remote sensing data: the FUD date, BUD date, CID, and ID obtained by passive microwave data are mostly shorter than those obtained by MODIS data (as shown in Supplementary Table S2).  The boxplot shows the magnitude of the difference between the datasets. Through the different remote sensing data types, the FUS determined by passive microwave is generally later than the results obtained by MODIS (Figure 2a), and this can be attributed to the freezing and thawing patterns of lakes. Many lakes begin to freeze in the shallow water near the lakeshore. The MODIS pixels with higher spatial resolution cover the entire lake surface, and this change is more easily observed, while passive microwave pixels, with a coarse resolution and locating at the center over deep water, have less interference from the surrounding land surface. Differences in the spatial distribution of pixels lead to differences in the observed timing of lake ice formation, which is similar to the results of Cai et al. [19,50]. In addition, the initially formed thin ice may not be easily detected by the brightness temperature threshold [11], resulting in a later FUS determined by passive microwave data. Meanwhile, most of the FUE (BUS) dates obtained by passive microwave are earlier (later) than MODIS (Figure 2b,c), which may be since the lake ice ratio corresponding to the brightness temperature threshold is smaller than that used by MODIS to determine lake ice phenology parameters. Moreover, compared with the FUS, FUE, and BUS dates, the consistency of the BUE date is better, and there is no particularly large difference (Figure 2d), which is consistent with the conjecture of Cai et al. that the impact of different remote sensing data on the freezing date may be greater than that on melting dates [50]. However, in general, the BUE date extracted by passive microwave is earlier than that of MODIS data, which may be related to the earlier ice melting in deep water. The calculation of duration also reflects the differences in remote sensing data: the FUD date, BUD date, CID, and ID obtained by passive microwave data are mostly shorter than those obtained by MODIS data (as shown in Supplementary Table S2).  Table 1 [18,[30][31][32]47,50,51].
The above analysis shows that there are differences in lake ice monitoring by using different remote sensing data, and the multi-sources average can balance the differences and obtain general lake ice phenological results ( Table 1). The FUS of lakes on the Tibetan Plateau is mainly concentrated in November and December. Among them, Hoh Xil Lake and Ulan Ul Lake began to freeze at the earliest (27 October), while Xuru Co began to freeze at the latest (29 January), and arrived in late January of the following year, all lakes were in a completely frozen state (Nam Co was completely frozen at the latest, 29 January). In late February of the following year, the lakes began to melt one after another. Among them, Taro Co began to melt at the earliest time (19 February). By April and May, the lake ice of most lakes disappeared completely, but the ice in a very few lakes could last until June and July. The ice coverage duration varies widely, from 79 days (Xuru Co) to 239 days (Hoh Xil Lake).

Analysis of Influencing Factors of Lake Ice Phenology Parameters
We carried out a statistical analysis of the dependence of the average ice phenology parameters (Table 1) on climatic factors and local factors (Table 3) in 40 lakes on the TP. The analysis results showed that temperature had the most significant effect on the freeze-thaw dates. The correlation coefficient between "average temperature in December" and FUS can reach 0.77, and the correlation between "average temperature in June" and "BUE" was higher (r = −0.79). This confirms that lake phenology was highly determined by thermodynamic factors as reported by Kouraev et al. [33], with BUE more thermally determined [66]. There were high correlations between freeze-thaw dates and latitude (r = −0.71 for FUS and r = 0.41 for BUE), solar radiation (r = 0.54 for FUS and r = −0.16 for BUE), and altitude (r = −0.1 for FUS and r = 0.38 for BUE), which is attributed to the autocorrelation between these factors and temperature. In general, low altitudes, low latitudes and strong solar radiation correspond to high air temperatures. The correlation between snow cover and freeze-thaw dates were significant (r = −0.41 for FUS and r = 0.29 for BUE), because snow accumulation can affect the growth and thaw rates of the ice cover by adding insulation and albedo. Snowfall is beneficial to the initial formation of lake ice and slow down the melting process of lake ice, which is consistent with previous studies [35,[69][70][71][72][73]. The correlation between wind speed (in December) and freeze-thaw dates was small (r = −0.21 for FUS), mainly because the freezing and thawing of lakes is a fast-changing process, which can be completely frozen or thawed in about half a month, while the monthly wind speed contains too much information and does not reflect the characteristics of wind speed changes in the days before and after FUS and BUE. Overall, air temperature dominated during the freeze-thaw process, and BUE had a stronger dependence on air temperature than FUS, and snow was important in the freeze-thaw process. Through the above correlation analysis, the most important factors for lake ice phenology (elevation, latitude, area, average temperature in December, average temperature in June, shortwave radiation, and snow depth) were selected to perform cluster analysis on the lakes. The ice phenology dates of each group of lakes are a reflection of climate and lake attribute characteristics (as shown in Supplementary Figure S3, Table S3).

Trend Analysis of the Lake Ice Phenology Extracted by AMSR-E/2 from 2002 to 2021
Among the 40 lakes, six (Dogai Coring, Xijir Ulan Lake, Taro Co, Paiku Co, Tangra Yumco, Xuru Co) have unfrozen/incomplete freezing conditions, thus they are not considered in trend analysis. Among the remaining 34 lakes, 24 (70.6%) had later FUS and 12 were statistically significant at the 90% confidence level. The average rate of change in later FUS was 0.35 days/year, while only 7 (20.6%) lakes had earlier BUE at an average trend of −0.17 days/year. Ice cover duration (ID) was affected by both FUS and BUE. Among them, 16 (47.1%) lakes extended the ID at an average rate of 1.17 days/year, 13 (38.2%) lakes shortened ID at an average rate of −0.33 days/year, and 5 (14.7%) lakes had no obvious trend in ID ( Table 4). The above results do not seem to be consistent with the changing characteristics of the earlier BUE and the shorter ID under the background of global warming, and further discussion can be found in Section 4.1.

Effects of Snow Cover on Lake Ice Phenology and Trend Change
In order to understand the reasons for the later BUE of most lakes, the BUE variation curves of lakes from 2002 to 2021 were drawn. Taking Qinghai Lake as an example (Figure 3), the BUE of the last few years (2018-2021) was later than other years, which obviously increases the positive trend of BUE, and this phenomenon generally exists in Har Lake, Kusai Lake, Zhuonai Lake, Serling Co and other lakes. To illustrate the influence of the BUE anomaly from 2018 to 2021 on the overall change trend, the change trend of BUE of each lake from 2002 to 2018 was calculated (Table 5). Compared with 2002-2021 (Table 4), more lakes had later BUE and shorter ID from 2002 to 2018 (Table 5), and the proportion of earlier BUE increased from 20.6% to 41.2%. The later FUS is still common (76.5%), and correspondingly, more lakes had shorter ID with proportion increased from 38.2% to 64.7% ( Table 6). The trend analysis results (2002-2018) are consistent with Cai et al. They used MODIS snow products to study the phenological changes of lake ice in 58 lakes in the TP region from 2000 to 2017, and found that 81.03% of the lakes have later FUS, 50% of the lakes have earlier BUE, and 68.79% of lakes have shorter ID [50]. Liu et al. [74] investigated the changing characteristics of lake ice phenology over the TP during 2002-2015 and found that the lakes in the southern TP had delayed break-up dates and prolonged ice durations, which was attributed to the impact of the winter North Atlantic Oscillation (NAO). Both studies indicate that the trend toward earlier BUE need not be dominant.
vestigated the changing characteristics of lake ice phenology over the TP during 2015 and found that the lakes in the southern TP had delayed break-up dates an longed ice durations, which was attributed to the impact of the winter North Atlan cillation (NAO). Both studies indicate that the trend toward earlier BUE need not b inant.     The delay of BUE is likely to be related to the increase of snow cover [69,71]. To confirm this idea, a long-term series of daily snow depth data in China (1979-2021) was selected for verification. From 2002 to 2020, the average snow days with a snow depth ≥ 5 cm from December to May of the next year in each region were drawn (Figures 1 and 4). It can be clearly seen that 2018 and 2019 were snowy years, which correspond to the later BUE in 2018 and 2019. Therefore, the trend towards later BUE in most lakes from 2002 to 2021 is likely to be caused by snow cover: from 2018 to 2020, the abnormal increase in snow cover led to abnormally late BUE, which increased the trend of delaying BUE, thus covering up the trend of advancing BUE.

Effect of Pixel Position on Lake Ice Phenology Parameters
It is worth noting that the difference in pixel position also affects the results of lake ice phenology. Taking Nam Co as an example, two pixels were selected in the east and west of the lake (the distance between the two is about 30 km), and the lake ice phenology parameters from 2002 to 2021 (except 2011 to 2012) were extracted respectively (as shown in Supplementary Tables S4 and S5). The pixel position will affect the judgment of whether the lake is completely frozen. The difference of break-up dates caused by the difference of pixel position was also obvious, which can reach up to 20 days, and the break-up dates of eastern pixel were later than that of the western pixel ( Figure 5).

Effect of Pixel Position on Lake Ice Phenology Parameters
It is worth noting that the difference in pixel position also affects the results of lake ice phenology. Taking Nam Co as an example, two pixels were selected in the east and west of the lake (the distance between the two is about 30 km), and the lake ice phenology parameters from 2002 to 2021 (except 2011 to 2012) were extracted respectively (as shown in Supplementary Tables S4 and S5). The pixel position will affect the judgment of whether the lake is completely frozen. The difference of break-up dates caused by the difference of pixel position was also obvious, which can reach up to 20 days, and the break-up dates of eastern pixel were later than that of the western pixel ( Figure 5).
December to May of the following year. The division of regions is shown in Figure 1.

Effect of Pixel Position on Lake Ice Phenology Parameters
It is worth noting that the difference in pixel position also affects the results of lake ice phenology. Taking Nam Co as an example, two pixels were selected in the east and west of the lake (the distance between the two is about 30 km), and the lake ice phenology parameters from 2002 to 2021 (except 2011 to 2012) were extracted respectively (as shown in Supplementary Tables S4 and S5). The pixel position will affect the judgment of whether the lake is completely frozen. The difference of break-up dates caused by the difference of pixel position was also obvious, which can reach up to 20 days, and the break-up dates of eastern pixel were later than that of the western pixel ( Figure 5).

Uncertainty Analysis
Insufficient spatial and temporal resolution of remote sensing data and the limitations of the methods brought uncertainties in the extraction of lake ice phenology parameters. In this study, the absence of remote sensing data will directly affect the extraction results. In addition, the coarse spatial resolution (0.1 • ) will result in the brightness temperature curve being contaminated by lakeshore radiation, while the degree of influence of radiation pollution on the results cannot be quantified. The limitations of the brightness temperature difference search method have been pointed out in some other studies [26,27]; when the jumping phenomenon of high platform values on the brightness temperature curve is not obvious, the results obtained through the method of searching for the maximum and minimum are most likely not located at the initial position of high platform values, which may also appear in our estimation. As for some unfrozen/not completely frozen lakes, they cannot be automatically identified by the program, and can only be judged by the method of visual interpretation of the brightness temperature curve.
The available time series of AMSR-E/2 data allowed us to derive trends in lake phenology parameters for the studied lakes, however, the results of trend analysis are easily affected by fluctuations in individual years due to short time series. For example, advanced BUE events appear more when heavy snow years are removed. Therefore, a combination of satellite products (optical and passive microwave data) and numerical simulations is suggested for long-term trend analysis.
In the process of merging the eight lake ice phenology datasets, the results of lake ice phenology derived from different remote sensing data were very different. For example, the average FUS date of Dogai Coring determined by optical remote sensing and passive microwave remote sensing can differ by 36 days [30,47]. The average FUE (BUS) date of Migriggyangzham Co (Xijir Ulan Lake) can differ by 30 (74) days [18,31]. As five of the eight datasets came from MODIS products, the average results will contain more optical information and the RMSEs of passive microwave data are larger. In addition, the study period, and the number of lakes in each dataset were different. Consequently, lakes with more datasets got more representative averages, while lakes with fewer datasets were more susceptible to individual results.
The frozen-thaw process will be affected by atmospheric heat exchange, water heat storage, solar radiation, ice and snow conditions, wind, inflow, etc. [69]. However, water depth as the function of lake heat storage was not considered in this study. Although the effects of wind and radiation were considered, the data used do not reflect changes in wind speed around freeze-thaw dates and the inter-annual variations in radiation.

Conclusions
This paper uses the AMSR-E/2 passive microwave brightness temperature data from 2002 to 2021 (except 2011 to 2012) to extract lake ice phenology parameters for lakes in the TP region, and the MODIS 1, 2 and 3 band false color composite images were used to verify the unfrozen/incompletely frozen lakes, and the two had good consistency. An analysis of the lake ice phenology trend of 34 lakes that were completely frozen in each of year from 2002 to 2021 found that 70.6% of the lakes had a later FUS, and only 20.6% of the lakes had earlier BUE. Considering that snow cover may be the reason for the abnormal trend of BUE, the last few snowy years (2018-2020) were eliminated and the trend analysis was carried out, and it was found that the average FUS of 34 lakes in 2002-2018 was 4 December, the trend towards later FUS was still common (76.5% of the lakes, average change rate = 0.59 days/year), and a few lakes (23.5%) have a trend towards earlier FUS (average change rate = −0.68 days/year). The average BUE was 7 May, and the proportion of lakes with earlier BUE increased to 41.2% (average change rate = −0.38 days/year), but lakes with a lagging trend in BUE still dominated (52.9%, average change rate = 0.79 days/year). The average ID affected by FUS and BUE was 153 days, 64.7% of lakes had shorter ID (average change rate = −0.76 days/year), and 35.3% of lakes had longer ID (average change rate = 1.24 days/year) ( Table 6).
Comparing the results of this study with the other seven lake ice phenology datasets, it was found that the type of remote sensing data has an influence on the results of lake ice phenology parameters. In general, the FUS, BUS (FUE, BUE) dates extracted from passive microwave data were later (earlier) than the extraction results of MODIS data. The FUD date, BUD date, and CID and ID values of lakes determined by passive microwave data were almost shorter than the results from the MODIS data, but the specific differences will also be affected by the algorithm, lake properties, etc.
The lake temperature profile measurements provide a new validation dataset for the determination of lake ice phenological parameters. The FUE and BUE dates determined by lake temperature were in good agreement with remote sensing results, but FUS and BUS dates were significantly earlier than those derived from remote sensing. In order to balance the influence of remote sensing data types and algorithms on the extraction results of lake ice phenology parameters, the results of each dataset were synthesized to obtain the average ice phenology results of lakes in the TP area. It was found that the TP lakes mainly began to freeze in November and December, and the lake ice completely melted in April and May of the following year. However, the freezing and thawing dates of lakes vary greatly, mainly due to different meteorological conditions and local factors. Through the correlation analysis of lake ice phenology influencing factors, it is found that temperature is dominant, and the dependence of BUE date on temperature is higher than that of the FUS date. Considering several factors (elevation, latitude, area, average temperature in December, average temperature in June, shortwave radiation, and snow depth) that have a greater impact on lake ice phenology, cluster analysis of lakes found that the lakes in high altitude, high latitude and low temperature areas have the characteristics of early freezing, late melting, and long lake ice duration.
Although the application of deep learning and machine learning algorithms in lake ice research is very limited, its advantages have been shown. For example, the accuracy performance of convolutional neural network (CNN) is greater than 92% [75]. These automatic classification methods are more suitable for lake ice identification than previous methods (such as the threshold method). Therefore, advanced data processing methods should be introduced in the future study of lake ice phenology on the TP, such as, k-means algorithm, Iterative Region Growing and Semantics (IRGS) algorithm, support vector machine (SVM) algorithm, the genetic algorithm (GA), and convolutional neural network (CNN), etc.    Table S1: Unfrozen/incompletely frozen lakes during 2002 to 2021. NF is the abbreviation for not frozen, and NCF is the abbreviation for not completely frozen. Blank means frozen. Table S2: Average freeze-up duration (FUD), break-up duration (BUD), lake ice cover duration (ID), and completely ice cover duration (CID) of each lake ice phenology dataset. Time between freeze-up start date and freeze-up end date determines the FUD, while time between break-up start date and break-up end date determines the BUD. ID is defined as the period between freeze-up start date and break-up end date, and CID is defined as the period between freeze-up end date and break-up start date. All variables are in "days". FID numbers indicate the lakes, with lake information shown in Table 1. Table S3: Average values of basic characteristics and ice phenology for each lake group. The number of lakes in each group (N), altitude (Alt), latitude (Lat), lake area (Area), average temperature in December (T-12), average temperature in June (T-6), annual shortwave radiation (SW), and snow depth (SD) are considered. Lake ice phenology parameters include freeze-up start date (FUS), freeze-up end date (FUE), break-up start date (BUS), break-up end date (BUE), and the period between freeze-up start date and break-up end date (ID). Table S4: Lake ice phenology parameters during 2002 to 2021 in Nam Co extracted by the western pixel. Table S5: Lake ice phenology parameters during 2002 to 2021 in Nam Co extracted by the eastern pixel. (References [18,25,[30][31][32]47,50,51]