Phenological Changes and Driving Forces of Lake Ice in Central Asia from 2002 to 2020

Abstract

Lake ice phenology is an indicator of past and present climate, it is sensitive to regional and global climate change. In the past few decades, the climate of Central Asia has changed significantly due to global warming and anthropogenic activities. However, there are few studies on the lake ice phenology in Central Asia. In this study, the lake ice phenology of 53 lakes in Central Asia were extracted using MODIS daily LST products from 2002 to 2020. The results show that MODIS-extracted lake ice phenology is generally consistent with Landsat-extracted and AVHRR-extracted lake ice phenology. Generally, lakes in Central Asia start to freeze from October to December. The trends in the lake ice phenology show strong regional differences. Lakes distributed along the Kunlun Mountains show overall delayed trends in all lake ice phenology variables, while lakes located in southwestern Central Asia show clear advancing trends in the freeze-up start dates (7.06 days) and breakup end dates (6.81 days). Correlations between the phenology of lake ice and local and climatic factors suggest that the ice breakup process and the duration of its complete coverage depend more on heat, while precipitation mainly affects the freezing time of the ice. Wind speed mainly affects the time of completely frozen of ice. In general, the breakup process is more susceptible to climatic factors, while local factors have strong influences on the freeze-up process.

1. Introduction

Global climate change is profoundly affecting human survival and development and is one of the major challenges faced by the international community today [1]. Historical records show that the warming trend began in the 1950s and changed abruptly in the mid to late 1980s [2]. In recent years, the global warming trend has become more pronounced. From 2011 to 2020, the average temperature has increased by 1.09 ${}^{\circ }$C over the pre-industrial period (1850–1900), while the 2001–2010 period was 0.99 ${}^{\circ }$C warmer than the pre-industrial period [3]. This abnormal warming trend has increased evaporation from the land surface and accelerated the melting of glaciers. The increased evaporation from the land surface has reduced the areas of lakes at lower elevations and plains, while the melting of glaciers has led to the formation of new lakes at higher elevations [4]. In the past decades, a widespread water loss occurred in the global endorheic system [5] and the total lake area in Central Asia has shown a decreasing trend [6,7,8]. Previous studies have shown that nine lakes in Central Asia have decreased in area by about 45,352 km${}^2$, accounting for approximately 49.62% of their total surface areas [6] and the number of lakes decreased from 7309 in 2000 to 6585 in 2015, a decrease of 85 lakes per year [7]. Thus, global climate change has greatly influenced the evolution of lakes.

Global climate change also has a significant impact on lake processes, such as changes in lake size, physical and chemical properties, and lake ice phenology. Among the above processes, lake ice phenology is the most important factor, which directly reflects climate change, thus making it a sensitive indicator of climate change [1,9,10,11,12,13,14]. Lake ice phenology describes the freeze–thaw cycles of lakes and the durations of their ice covers, regulating the chemical, biological, and physical processes in lakes [13,15,16,17,18]. Lakes are particularly important in arid and semi-arid regions. Lakes and their ice covers can provide a number of ecosystem services for people living in Central Asia, such as water supply, agricultural irrigation, and some recreational activities [19]. In addition, lakes play an important role in social and ecological systems in these inland areas where surface flow is scarce [5].

Prior to the widespread use of remote sensing, studies on lake ice phenology relied on traditional methods, including field observations and hydrological studies. However, data collections using these traditional methods were laborious and the data were spatially sparse [20]. Remote sensing data with different spatial and temporal resolutions can fill the data gaps of traditional methods over long time periods. Moderate-resolution imaging spectroradiometer (MODIS) data have temporal resolution of one day and can be important for the long-term monitoring of lake ice. Several MODIS products have been used to detect lake ice phenology, including the MODIS snow product [9,19] and reflectance product [21,22]. In this study, we used the MODIS Land Surface Temperature (LST) product because it contains both daytime and nighttime images and presents more information than other products.

Therefore, we extracted the lake ice phenology from 2002 to 2020 via Google Earth Engine (GEE) and geemap [23]. This study aimed to (1) quantify the changes in lake ice phenology over the past 18 years and (2) analyze the driving factors of changes in lake ice phenology.

2. Materials and Methods

2.1. Study Area and Lake Selection

The study area is located in the mid-latitudes of the Northern Hemisphere, with a spatial extent of 28° N, 61° N, 44° E, and 97° E (Figure 1). The weather in Central Asia is controlled by a continental climate with cold winters and severe droughts in the summer. The overall topography of Central Asia is high in the southeast and low in the northwest, and the numerous inland lakes are unevenly distributed. Annual precipitation in Central Asia is higher in mountainous areas (about 1000 mm) and lower in desert areas (less than 100 mm) [6]. However, in this region, potential evaporation is much higher than annual precipitation [24], making the region one of the world’s most important arid zones. Due to the unique climate, most of the lakes in Central Asia rely on glacial melt, mountain precipitation, and river runoff for their water supply.

2.2. Datasets

2.2.1. MODIS Daily Land Surface Temperature (LST) Products

MODIS is a key instrument installed on the Terra and Aqua satellites; it has provided data since 2000 and 2002, respectively. The MODIS sensor has 36 spectral bands covering wavelengths from 0.405 $\mathsf{\mu }$m to 14.385 $\mathsf{\mu }$m. Our primary data are daily land surface temperature products from MODIS (MOD11A1, MYD11A1, version 6) for the period from 2002 to 2020, obtained from NASA’s Earth Observing System Data and Information System website via the Google Earth Engine (GEE) platform (https://modis.gsfc.nasa.gov/data/dataprod/mod11.php (accessed on 15 January 2022)). The daily LST product retrieves land surface temperatures using band 31 and band 32 through the split-windows algorithm [25], including daytime LST and nighttime LST of the land surface.

Among the various sensors and products, we chose the MODIS sensor and the MODIS LST daily product because the MODIS sensor has a temporal resolution of 1 day and its LST daily product contains a daytime LST and a nighttime LST with a spatial resolution of 1 km, which can be directly used to create continuous time series for further analysis. The MODIS LST product has the unique advantage of providing more surface information (four images per day) compared to other MODIS products, which increases its priority in this study.

2.2.2. Lake Information

The HydroLAKES database contains freshwater lakes, salt lakes, and artificial reservoirs, and was developed by Global HydroLAB [26]. The main processes used to produce HydroLAKES include manual identification and removal of polygons from rivers and wetlands, removal of duplicate and overlapping polygons, and the correction of damaged or incomplete polygon geometry. Surface area, depth, volume, shoreline length, elevation, and other valuable information for lakes described in the HydroLAKES database are available on the HydroSHEDS website (https://www.hydrosheds.org/page/hydrolakes (accessed on 3 January 2022)).

2.3. Methods

The algorithm for calculating ice–water fraction mainly contains four steps. Firstly, input MODIS datasets and HydroSHEDS datasets; apply quality control to MODIS datasets to remove cloudy pixels and pixels with high LST error values and select lakes. Secondly, merge the lake ice time series data retrieved from the LST_Day_1km band and LST_Night_1km from MOD11A1 and MYD11A1 products from 2002 to 2020. Thirdly, classify lake ice through a LST threshold value and calculate the ice fraction; export the lake ice fraction time series in a csv file. Lastly, use the merged ice fraction time series to extract ice phenology variables, including freeze-up start (FUS), freeze-up end (FUE), breakup start (BUS), breakup end (BUE), freeze-up duration (FUD), ice completely covered duration (ICD), and the breakup duration (BUD).

In this study, FUS represents the first day when the ice fraction is above 0.2, FUE represents the first day when the ice fraction is above 0.8, BUS represents the last day when the ice fraction is above 0.8, and BUE represents the first day when the ice fraction is above 0.2 after BUS. FUD represents the day between FUS and FUE, ICD represents the day between FUE and BUS, and BUD represents the day between the BUS and BUE day of the year.

2.3.1. Quality Control and Lake Selection

Using optical sensors to retrieve lake ice phenology can be largely affected by the clouds. We removed cloud-covered pixels and high LST error pixels through the quality control in the QC_Day and QC_Night bands of MOD11A1 and MYD11A1 datasets. In this study, we set pixels with the quality flag of “LST produced”, “good data quality”, and “average LST error ≤ 1 K” as valid pixels, and classified other pixels as invalid pixels. To simplify the following calculations, we set the value of valid pixels to 1 and the value of invalid pixels to 0. The formula is given as follows in Equations (1) and (2):

$$Mas{k}_{\left(t\right)}=\left\{\begin{array}{ccc}\hfill Vali{d}_{Pix}& \hfill =\hfill & \hfill 1\\ \hfill Invali{d}_{Pix}& \hfill =\hfill & \hfill 0\end{array}\right.$$

(1)

$$Imag{e}_{\left(t\right)}^{(c,S)}=Imag{e}_{\left(t\right)}^{(r,S)}\times Mas{k}_{\left(t\right)}$$

(2)

where t is the day of the year of the image; S is the name of the satellite (Terra and Aqua); Image${}^r$ is the LST daily product of the Aqua and Terra satellite; Image${}^c$ is the image after the quality control operation. By applying Equations (1) and (2) to the whole image collection, we removed most of the low-quality pixels of the whole image collection.

To avoid noise introduced by the mixed pixels, we only chose lakes that had surface areas larger than 30 km${}^2$. After the selection process, there remained 53 lakes with surface areas ranging from 32.43 to 23,865.91 km${}^2$.

2.3.2. MODIS Daily LST Combination

Since clouds are constantly moving in the atmosphere, cloud coverage varies with time, making it possible to fill in missing values due to quality control operations by combining images from different times of the day. If a pixel is a cloudy pixel in one satellite observation and a cloudless pixel in another observation, the pixel value of the cloudless pixel is used to fill the cloudy pixel as the final result. Considering that the air temperature varies over time, which means that a valid pixel may have two LST values, we used the largest value as the pixel value. In addition, choosing the largest LST value as the pixel value also reduced the error of water–ice classification, because water freezes when the temperature is below 0 degrees Celsius. Therefore, we used the maximum value combination method to combine the data of MOD11A1 and MYD11A1 (of the same day) after quality control processing. The formula is listed as follows in Equation (3):

$$Merg{e}_{(x,y,t)}=Max(Imag{e}_{(x,y,t)}^{(c,T)},Imag{e}_{(x,y,t)}^{(c,A)})$$

(3)

where x is the index of column (horizontal axis); y is the index of row (vertical axis); T stands for the Terra satellite and A stands for the Aqua satellite. Equation (3) was applied to data from MOD11A1 and MYD11A1 products on the same day. Maximum LST values for daytime and nighttime were obtained from two satellite products.

2.3.3. Lake Ice Penology Variable Extraction

The time series of water–ice transitioning was extracted from the merged image in Section 2.3.2 for the period from 2002 to 2020. The water–ice phases in the MODIS daytime and nighttime merged images were classified by applying a temperature threshold. Most of the lakes in the CA region are saline, which has a lower freezing point. Thus, we selected 272 K as the temperature threshold. If the LST value of the pixel was lower than 272 K, we set the pixel to ice and classified it as $Pixe{l}_{ice}$; if the LST value was greater than 272 K, we set the pixel to water and classified it as $Pixe{l}_{water}$. The water–ice fraction was calculated by Equation (4):

$$Fractio{n}_t=\frac{N_{ice}}{N_{valid}}$$

(4)

where t is the index of time and includes daytime and nighttime; $N_{ice}$ is the total number of $Pixe{l}_{ice}$ and $N_{valid}$ is the total number of valid pixels. Equation (4) was applied to data after operation from Section 2.3.2. Daytime and nighttime lake ice fractions were obtained. To further reduce the error of the water–ice fraction, we combined two lake ice fraction time series using Equation (3).

3. Results

3.1. Validation of Ice Phenology Curve

The accuracy of derived ice–water fraction of a lake from the MODIS daily LST combined image was estimated by comparing the ice–water fraction derived from AVHRR daily SST product and Landsat 8. The overall accuracy was calculated as the root mean square error (RMSE) of the time difference and R${}^2$.

The lake ice fraction time series of three large lakes derived from the MODIS daily LST combined image, AVHRR daily SST image, and Landsat 8 image were compared to access the performance of the method. The lake ice fraction time series derived from the MODIS daily LST combined image was consistent with the AVHRR daily SST image and Landsat8 image (Figure 2 and Figure 3). The R${}^2$ of the comparison ranged from 0.68 to 0.81 and the root mean square error (RMSE) ranged from 0.17 to 0.25 (

Table 1. The R${}^2$ and RMSE for the comparison between the MODIS daily LST combined image, AVHRR and Landsat8.

	MODIS vs. AVHRR			MODIS vs. Landsat8
	Large Aral Sea	Lake Balkhash	Lake Zaysan	Large Aral Sea	Lake Balkhash	Lake Zaysan
R${}^2$	0.80	0.81	0.71	0.77	0.79	0.68
RMSE	0.17	0.19	0.23	0.18	0.19	0.25

), which show that the water–ice fraction derived from the MODIS daily LST combined image had good quality.

3.2. Lake Ice Phenology Characteristics in Central Asia

Remarkable differences in lake ice phenology were found in Central Asia lakes. Some lakes, such as Xiaoxihaiz Shuiku, Kara-Bogaz-Gol, Lake Issyk Kul, Lake Sarygamysh and Lake Aydar, were rarely fully ice-covered in the winter, while some other lakes were frozen all year around. For some lakes, such as Xiaoxihaiz Shuiku, compared to other time periods, the lake ice covers for the period from 2007 to 2008 were significantly increased. This may be due to the occurrence of the La Nina event in 2007. The mean ’day of year’ (doy) of FUS was 106.85 in all lakes; the mean doy of FUS was 137.02 in not completely ice-covered lakes, and the mean doy of FUS was 94.95 in completely ice-covered lakes (

Table 2. Lake ice phenology characteristics in Central Asia.

Lake Type	FUS (doy)	FUE (doy)	BUS (doy)	BUE (doy)
Not completely frozen in winter	137.02	NA	NA	180.32
Completely frozen in winter	94.95	125.26	227.35	252.20

NA means no data, notice that 8 August is referred to day 1 of the hydrological year.

); the mean doy of FUE was 125.26 in completely ice-covered lakes (Table 2); the mean doy of BUS was 227.35 in completely ice-covered lakes (Table 2); the mean doy of BUE was 231.86 in all lakes, the mean doy of BUE was 180.32 in not completely ice-covered lakes, and the mean doy of BUE was 252.20 in completely ice-covered lakes (Table 2).

Lake Seletyteniz had the mean earliest doy of FUS (64.50), Lake Issyk Kul had the mean latest doy of FUS (165.17), Lake Seletyteniz had the mean earliest doy of FUE (89.89), Lake Sarezskoye had the mean latest doy of FUE (165.67), Lake Kaydak had the mean earliest doy of BUS (179.61), Lake Arkatag had the mean latest doy of BUS (271.5), Lake Saryyazin Skoye had the mean earliest doy of BUE (153.33), and Lake Arkatag had the mean latest doy of BUE (310.17).

Previous studies have revealed that the lakes in the Northern Hemisphere experienced earlier breakups and later freeze-ups. Not all lakes in Central Asia had the same trend as previous studies [15,27]. The trends of the ice phenology variables are presented in Figure 4. From 2002 to 2020, 26 lakes (49.06%) had advanced FUS dates, and 21 lakes (39.62%) had advanced BUS dates. Due to some lakes rarely being completely frozen in the winter and the ’day of year’ (in FUE and BUS) was too noisy for analysis, we set the trend value to zero in FUE and BUS. Lakes, located in the southwest part of Central Asia had obvious advancing trends in FUS dates. For the BUS dates, lakes in the north and west parts of Central Asia had obvious advancing trends.

3.3. Ice Phenology in Different Elevation Lakes

Lakes were clustered into three groups according to the elevation using K-means (Group 1: elevation higher than 3000 m, Group 2: elevation from 1000 to 2000 m, Group 3: elevation lower than 1000 m). The mean temporal patterns and trends of ice phenology of three groups of lakes from 2002 to 2020 are presented in Figure 5 and

Table 3. Phenology variable change and trend.

Group	FUS		FUE		BUS		BUE		Elevation (m)
	Mean Doy	Trend (Day/Year)	Mean Doy	Trend (Day/Year)	Mean Doy	Trend (Day/Year)	Mean Doy	Trend (Day/Year)
1	87.56	0.13	113.67	1.04 **	254.09	0.55	284.74	0.67 **	3925
2	125.22	0.08	135.08	0.71 **	220.32	−0.65 **	231.20	0.15	1535
3	108.14	−0.35	137.95	0.40	210.86	−0.82 **	217.60	−0.31	207

** represents the statistical significance at the 0.05 level.

. Generally, higher altitude lakes had earlier FUS, FUE dates and later BUS, BUE dates. Compared to Lake Kaydak (elevation −29 m), in group 3, the FUS and FUE in Lake Aksayquin (elevation 4844 m) in group 1 were 27 and 41 days earlier, while the BUS and BUE were 80 and 74 days later. Furthermore, the trends of each lake’s ice phenology variable in the three groups were consistent overall with the trends of earlier breakups and later freeze-ups in the Northern Hemisphere, as revealed by previous studies [12,15,27,28]. The trend of ice phenology in group 2 correlated with previous studies, while the trends in group 1 (an earlier breakup trend) and group 3 (a later freeze-up trend) existed some differences (see Table 3). The FUD, ICD, and BUD are listed in

Table 4. Lake ice process duration change and trend.

Group	FUD		ICD		BUD
	Mean Duration	Trend (Day/Year)	Mean Duration	Trend (Day/Year)	Mean Duration	Trend (Day/Year)
1	26.74	0.84 **	140.42	−0.49	31.34	0.05
2	23.46	−0.02	109.42	−0.71 *	21.57	0.08
3	32.90	0.02	89.18	−0.41	24.71	0.31

** represents the statistical significance at the 0.05 level and * represents the statistical significance at the 0.1 level.

.

3.3.1. Group 1: Elevation Higher than 3000 m

Lake Karakul in group 1 had the longest average FUD of 33.06 days (see Appendix A

Table A1. Mean Lake ice phenology of 53 lakes in Central Asia.

Lake Name	Mean FUS	Mean FUE	Mean BUS	Mean BUE	Mean FUD	Mean ICD	Mean BUD
Akkol	88.67	129	213.44	242.16	40.33	84.44	28.72
Aksayquin	82.06	99.61	259.39	286.89	17.56	159.78	27.5
Alakol	125.39	144.89	219.06	236.28	19.5	74.17	17.22
Aqqikkol	88.39	101.72	242.33	280.33	13.33	140.61	38
Arkatag	77.56	92.67	271.5	310.17	15.11	178.83	38.67
Ayakkum	94.44	105.44	234.5	256.39	11	129.06	21.89
Aydar	146.78	NA	NA	164.28	NA	NA	NA
Balkhash	107.78	124.06	229.06	243.11	16.28	105	14.06
Barkol	67.5	90.89	218.72	248.61	23.39	127.83	29.89
Biylikol	132.17	NA	NA	178.06	NA	NA	NA
Bosten	122.72	143.27	220.89	234.5	20.56	77.62	13.61
Chardarinskoye	127.72	NA	NA	172.89	NA	NA	NA
Chatyr-kel’	71.72	104.39	261.72	294.22	32.67	157.33	32.5
Ebi	102.89	124.06	209.78	231.33	21.17	85.72	21.56
Inder	102.39	NA	NA	224.78	NA	NA	NA
Issyk Kul	165.17	NA	NA	190.33	NA	NA	NA
Jili	110.83	136.56	221.94	250.56	25.72	85.38	28.61
Kamyshlybas	111.17	144.67	196.22	231.61	33.5	51.55	35.39
Kara-Bogaz-Gol	144.67	NA	NA	170.39	NA	NA	NA
Karakul	93	126.06	266.06	289.5	33.06	140	23.44
Karasor	84.06	120.06	228.67	254.28	36	108.61	25.61
Karatomarskoye	88.39	139.89	230.78	257.06	47.5	90.89	26.28
Kaydak	109.67	141.89	179.61	212.72	32.22	37.72	33.11
Kayrakkumskoye	134.39	NA	NA	159.83	NA	NA	NA
Kayrakkumskoye1	139.11	NA	NA	173.25	NA	NA	NA
Kilibek	71.22	121.39	207.94	244.39	50.17	86.55	36.44
Kishi’karoy	76.67	118.22	219.67	245.39	41.56	101.45	25.72
Kushmurun	73.5	113.78	236.33	254.61	40.28	122.55	18.28
Kuyumazarskoye	149	NA	NA	166	NA	NA	NA
Large Aral Sea	108.56	129.83	211.56	227.22	21.28	81.73	15.67
Manas	93.78	119.28	209.33	230.83	25.5	90.05	21.5
Markakol’	97.83	119.89	264.44	283.33	22.06	144.55	18.89
Qapshaghay Bogeni Reservoir	133.33	164.22	193.06	219.11	30.89	28.84	26.06
Sailimu	129.56	157.39	245.06	268.94	27.83	87.67	23.89
Sarezskoye	99.06	165.67	244.78	269.33	66.61	79.11	24.56
Sarygamysh	150.33	NA	NA	194.28	NA	NA	NA
Saryyazin Skoye	136.06	NA	NA	153.33	NA	NA	NA
Sasykkol	105.78	126.44	221.06	240.94	20.67	94.62	19.89
Seletyteniz	64.5	89.89	235.17	251.89	25.39	145.28	16.72
Shalkar	101.06	157.11	208.78	236.89	56.05	51.67	28.11
Shalkarteniz	72.61	123.44	212.11	237.28	50.83	79.94	25.17
Small Aral Sea	110.78	132.17	215.72	240.39	21.39	83.55	24.67
Song-kel	94.22	113.78	252.44	291.11	19.56	138.66	38.67
Sor Khaki	97.5	NA	NA	225.28	NA	NA	NA
Tengiz	74.61	96.67	234	252.28	22.06	137.33	18.28
Toktogul’skoye	154.28	NA	NA	197.11	NA	NA	NA
Tudakol	144.56	NA	NA	167.89	NA	NA	NA
Ul’ken-karoy	72.28	131.06	224.94	244.89	58.78	93.88	19.95
Ulungur	108.28	131.33	230.94	252.56	23.06	99.61	21.62
Uyaly	113.39	148.44	208.72	232.17	355.06	60.28	23.45
Xiaoxihaiz Shuiku	139.5	NA	NA	195.56	NA	NA	NA
Zaysan	97.89	116.56	241.17	255.22	18.67	124.61	14.05
Zhalauly	82.83	118.17	218.56	245.06	35.33	100.39	26.5

) and the only negative trend of FUD, this might be because Lake Karakul had the second largest area and depth (except for Lake Sarezskoye because the duration of Lake Sarezskoye was too noisy to analyze). Lake Aqqikkol had the shortest average FUD of 13.33 days (see Appendix A Table A1). Lake Arkatag had the longest average ICD of 178.83 days (see Appendix A Table A1), this might be because Lake Arkatag had the shallower lake depth in group 1. Lake Sarezskoye had the shortest average ICD of 79.11 days (see Appendix A Table A1) and the most negative trend of ICD in group 1. Lake Ayakkum had the shortest average BUD of 21.89 days, while Lake Arkatag had the longest average BUD of 38.67 days (see Appendix A Table A1).

3.3.2. Group 2: Elevation from 1000 m to 2000 m

Bosten Lake had the shortest average FUD of 20.56 days (see Appendix A Table A1) and Lake Sailimu had the longest average FUD of 27.83 days (see Appendix A Table A1). For ICD, Bosten Lake had the shortest average ICD of 77.62 days (see Appendix A Table A1), while Markakol had the longest average ICD of 144.55 days (see Appendix A Table A1) and the only positive trend of ICD. For BUD, Bosten Lake had the shortest average BUD of 13.61 days (see Appendix A Table A1), while Lake Barkol had the longest average BUD of 29.89 days (see Appendix A Table A1).

3.3.3. Group 3: Elevation Lower than 1000 m

Lake Balkhash had the shortest average FUD of 16.28 days (see Appendix A Table A1), while Lake Ul’ken-karoy had the longest average FUD of 58.78 days (see Appendix A Table A1). For ICD, the Qapshaghay Bogeni Reservoir had the shortest average ICD of 28.84 days (see Appendix A Table A1), while Lake Seletyteniz had the longest average ICD of 145.28 days (see Appendix A Table A1). For BUD, Lake Balkhash had the shortest average BUD of 14.06 days (see Appendix A Table A1), while Lake Shalkar had the longest average BUD of 28.11 days (see Appendix A Table A1).

4. Discussion

4.1. Correlation of Lake Ice Phenology Variables with Influencing Factors

Both climatic conditions and local factors affect lake ice phenology to varying degrees [15]. Climatic conditions control the heat exchange between lake and air [29], while local factors influence the heat capacity, which ultimately affects the phenology of lake ice. For example, climatic conditions, such as wind speed, can destroy thin ice when a lake with thin ice encounters strong winds [30]. In addition, lakes with larger volumes and areas require more heat to freeze and breakup during a hydrological year [19,31]. Major climatic factors include a mean 2 m of air temperature, mean temperature during deglaciation, wind speed, short-wave radiation, and precipitation [29,30,32,33]. Local factors include latitude, longitude, elevation, depth, and area [19,30,33]. Climatic factors including air temperature, wind speed, and precipitation, are derived from monthly average products of ERA5 and shortwave radiation from TerraClimate.

Figure 6a shows that temperature-related climatic factors mainly affect the ice breakup period, which is consistent with previous studies [19,29,34], and the completely frozen period of ice. Precipitation is negatively correlated with FUS and positively correlated with FUE, which has a strong influence on the completely frozen period of ice. This may be due to the fact that precipitation not only decreases the temperature [35,36], but also destroys the thin ice layer on the lake surface, which makes the freezing start earlier and the freezing end later, and finally prolongs the freezing process [30]. Short-wave radiation has a strong effect on FUS and FUD. The higher the wind speed, the earlier the freezing event occurs, the later the breakup event occurs, and the longer the ice completely covers the duration. This may be because wind accelerates convection over the lake surface, bringing cold air above the lake surface, accelerating the freeze-up process, delaying the breakup process [30], and prolonging the ICD. Strong winds may destroy the thin ice layer on the lake surface and prolong the freezing process [34,35,37].

Figure 6b,c, show the correlation between lake ice phenology variables and lake-specific factors. The results in Figure 6b shows that longitude and altitude have a strong influence on BUE, BUS, and FUE. This may be due to the fact that the higher the altitude, the lower the air temperature [38,39], and the altitude of Central Asia is higher in the east and lower in the west, which delays BUS, BUE, and advances FUE. Compared to BUE, BUS, and FUE, FUS is more dependent on latitude, lake area, and lake depth. This is because the temperature at high latitudes is lower than that at low latitudes, making FUS occur earlier. In addition, the larger the lake area and the deeper the average depth of the lake means that more heat needs to be released during the freezing period, thus affecting the date of FUS [9,30,34,40]. The result in Figure 6c shows that the lake area and the length of the lake shore have strong influences on BUD. This may be due to the fact that lakes with large areas and long shores have thinner ice compared to other lakes, which may be the reason for the negative correlation between BUD and lake area and shore length. FUD is negatively correlated with longitude and altitude. This is because the altitude in Central Asia is lower in the west and higher in the east, and the temperature is lower at higher altitudes than at lower altitudes. However, the anomalous relationship between latitude, lake area, and FUD is difficult to explain, probably because of the high salinity of lakes in Central Asia. ICD has a negative correlation with latitude and a positive correlation with longitude and altitude. This may be because the elevation in Central Asia is higher in the south and lower in the north, and elevation may be a more important factor than latitude.

4.2. Limitation and Prospection

The limitations of the study are listed below. First, the daily LST combination process (see Section 2.3.2) was unable to fill the data gaps caused by severe cloud contamination. Second, due to the lack of lake salinity data, using a threshold of 272 K to determine the freezing point of all lakes, instead of using a threshold based on lake salinity for each lake, may lead to some uncertainties in extracting lake ice phenology. In addition, the merging of daytime LST and nighttime LST ignores the temperature variation and may introduce some additional uncertainties [29]. Moreover, lake salinity is an important factor affecting lake ice phenology [33]. Due to the lack of lake salinity, the analysis between lake ice phenology and incluence factors (e.g., lake-specific factors) may not be comprehensive. Furthermore, anthropogenic activities affecting lake ice phenology were not addressed in this study. However, anthropogenic activities are also important factors affecting climate change, as anthropogenic global warming is best estimated at 1.07 ${}^{\circ }$C [3], which has a significant impact on terrestrial systems [41,42]. Finally, the anomalous relationship between latitude, lake area, and FUD is difficult to explain, probably because of the high salinity of lakes in Central Asia.

As a sensitive indicator of regional climate change, the phenology of lake ice is still inadequately studied in some regions, especially in Central Asia. In addition, satellite products with long-term and temporal continuities have relatively low spatial resolutions and are weak in monitoring small lakes, while satellites with higher spatial resolutions are not continuous in time. Further studies can be conducted in the following aspects:

(1)

Using multiple satellite data to retrieve lake ice phenology for all lakes, including small lakes.

(2)

Quantifying and evaluating how anthropogenic activities affect lake ice phenology variables.

(3)

Projecting future lake ice variations in different representative concentration pathways (RCPs).

5. Conclusions

In this study, we used a combined method applied to MOD11A1 and MYD11A1 products to obtain daily LST and classified water and ice pixels using temperature thresholds. Based on this approach, we extracted lake ice phenology for 53 lakes in Central Asia for 2002–2020, including FUS and BUE for 53 lakes and FUE and BUS for 38 lakes. The results show that the time series of water–ice components derived from the MODIS daily-combined LST dataset and the AVHRR and Landsat8 datasets are generally highly consistent.

In this study, lakes in Central Asia start to freeze from October to December. The mean FUS and BUE of the Central Asia lakes were 106.85 and 231.86, respectively. Trends in ice phenology for each lake showed strong regional differences. Lakes distributed along the mountains show overall delayed trends among all lake ice phenology variables. In contrast, lakes in southwestern Central Asia showed clear trends of advancing FUS and BUE. From 2002 to 2020, the average advancement rate of FUS was −0.29 days/year for the 53 lakes in Central Asia, the average delay rate of FUE was 0.24 days/year for the 38 lakes in Central Asia, the average advancement rate of BUS was −0.22 days/year for the 38 lakes in Central Asia and the average advancement rate of BUE was −0.25 days/year for the 53 lakes in Central Asia.

Lake ice phenology in Central Asia is influenced by climatic factors (temperature, precipitation, ice age temperature, wind speed, and shortwave radiation) and local factors (latitude, longitude, lake area, mean lake depth, shore length, and altitude). In general, temperature and wind speed have the most direct effects on lake ice phenology, and precipitation mainly affects FUD. Lakes at high elevations tend to have longer ICD. Lake areas mainly affect FUD and BUD.