Analysis of ice phenology of lakes on the Tibetan Plateau from MODIS data

The Tibetan Plateau includes a large system of endorheic (closed basin) lakes. Lake ice phenology, i.e. the timing of freeze-up and break-up and the duration of the ice cover may provide valuable information about climate variations in this region. The ice phenology of 59 large lakes on the Tibetan Plateau was derived from Moderate Resolution Imaging Spectroradiometer (MODIS) 8-day composite data for the period from 2001 to 2010. Ice cover duration appears to have a high variability in the studied region due to both climatic and local factors. Mean values for the duration of ice cover were calculated for three groups of lakes defined by clustering, resulting in relatively compact geographic regions. In each group several lakes showed anomalies in ice cover duration in the studied period. Possible reasons for such anomalous behaviour are discussed. Furthermore, many lakes do not freeze up completely during some seasons. This was confirmed by inspection of high resolution optical data. Mild winter seasons, large water volume and/or high salinity are the most likely explanations. Trends in the ice cover duration derived by linear regression for all the studied lakes show a high variation in space. A correlation of ice phenology variables with parameters describing climatic and local conditions showed a high thermal dependency of the ice regime. It appears that the freeze-up tends to be more thermally determined than break-up for the studied lakes.


Introduction
The Tibetan Plateau (TP) was identified as one of the most sensitive regions in the world to changes in climate (Liu and Chen, 2000).At the same time, because of its unique properties, the TP has a strong impact on the global climate system, acting as a large heat source and moisture sink during the summer (Ueno et al., 2001;Hsu and Liu, 2003;Sato and Kimura, 2007).It impacts the monsoon circulation as a barrier to zonal and meridional air motion (e.g.Kurzbach et al., 1993;Barry, 2008).Recently, many studies have been focusing on various aspects of climate change impact on the TP, yet our knowledge of its mechanisms and regional variations is still limited.An unbiased analysis of the impact of the changing climate on the TP is often hindered by lack of ground measurements.Sparseness of meteorological stations, rough conditions, and difficult accessibility call for the use of remote sensing methods.Variations in climate in this area are reflected by several phenomena that are well documented by satellite records such as snow distribution (Shreve et al., 2009;Xu et al., 2009;Kropáček et al., 2010), changes in glacier extent (Ye et al., 2006, Bolch et al., 2010), and fluctuations in lake area (Bianduo et al., 2009), which can therefore be seen as climate indicators.
An important indicator of changes in climate derived from satellite data is lake ice phenology.This has been proved by a number of recent studies on regional (Ruosteenoja, 1986;Barry and Maslanik, 1993;Duguay et al., 2006;Che et al., 2009) and global scales (Walsh et al., 1998;Magnuson at al., 2000).Interannual variation in lake ice phenology allows estimates of local climatic variability (Walsh et al., 1998).The lakes on the TP are distributed mainly in the central and western part, while meteorological stations are concentrated in the eastern part of the plateau (Liu and Chen, 2000).Derivation of lake ice phenology of Tibetan lakes could therefore fill a gap in our knowledge about the impact of climate change in this remote region.

Lake ice phenology
Ice phenology deals with periodic formation and decay of ice cover over water bodies and changes in its timing as a result of seasonal and interannual variations in climate.Formation of ice cover affects water bodies in high latitudes and altitudes, where the temperature during cold season falls below 0 • C for a sufficiently long time period.Appearance of ice is initiated when the lake surface water temperature falls below 0 • C. The date when detectable ice appears is referred to as the Freeze Onset (FO) in this study.The ice formation usually begins along the shore of shallow bays.Ice growth, especially during the formation of a new thin ice, can be interrupted by strong wind events.The date of the end of the freeze-up period, when the lake is fully ice-covered, is denoted as Freeze-up date (FU).An increase in temperature of the air above and the water below the ice in spring leads to a thinning of ice, eventually resulting in the ice breaking up.Appearance of a detectable ice-free water surface, after which no more total freeze-up occurs, is called Break-up date (BU).Further deterioration of ice leads to an eventual disappearance of ice, denoted as Water Clean of Ice (WCI).
Since there are discrepancies in the definition of variables describing the duration of lake ice cover, a comparison of ice phenology records for different regions is often biased (Brown and Duguay, 2010).Here we use two variables describing the duration of ice cover: Duration of Ice cover (DI), which results from subtraction of the WCI and FO, and Duration of Complete Ice cover (DCI), which is defined as the period between BU and FU.Many lakes on the plateau freeze up as late as January or February.A freeze-up in December or November is counted as ice cover for the following year for sake of simplicity, since it belongs to the same ice cycle.The choice of lake-ice phenology variables as a climate indicator varies among authors.Liston and Hall (1995) suggest that FU and BU, together with maximum ice thickness, are useful indicators of regional climate change, while Latifovic and Pouliot (2007) use the end of freeze-up and the end of break-up in their analysis of lake ice phenology in Canada instead; however, their choice was influenced by the availability of comparable in situ data.

Factors influencing the lake ice phenology
Sensitivity of lake phenology to climate variables has been investigated by several authors, where it was proved that BU/FU are highly sensitive to changes in air temperature.Using the results of a numerical simulation, a delay of as much as 4 weeks in freeze-up and break-up dates has been found by Liston and Hall (1995) for St Mary Lake in Canada in reaction to an increase in air temperature of 4 • C.They also attribute a sensitivity of ice phenology to a decrease in wind speed, which amounts to 2 weeks delay for Great Slave Lake in Canada.Air temperature has been suggested to be the most significant factor influencing the ice phenology, for instance by Livingstone (1997) and Kouraev et al. (2007).Analysis of ice phenology of numerous lakes and rivers in the Northern Hemisphere showed that a 1 • C change in mean air temperature corresponds to a change of about 5 days in the duration of the ice cover season (Skinner, 1993;Magnuson et al., 2000).However Weyhenmeyer et al. (2004) suggested that the relationship between temperature and break-up date is non-linear, leading to more sensitive reactions of ice breakup date, and to temperature changes in cold regions; however, it is not only temperature that determines the ice regime of lakes.There are several factors that affect both the duration of ice cover and the timing of ice phenology events, which can be divided into local and climatic factors.Apart from temperature, climate-related factors include wind speed, radiation and snow cover, the latter acting as an insulating layer.The absence of snow cover thus results in thicker ice cover, and leads in turn to later break-up (Kouraev et al., 2007).Local factors that are related to the ice regime of lakes are: lake bathymetry, salinity, altitude, shape of the shoreline, and lake area.The relative importance of particular factors for the lake-ice regime can be assessed using numerical models.For instance, Ménard et al. (2002) showed by using a onedimensional thermodynamic lake ice model that lake depth is a determinant factor of freeze-up date for Canadian lakes.

Effects of lake ice phenology on local or regional climate
While the lake ice phenology is impacted by climate change, the lake ice regime can in turn influence the local or even the regional climate.Freeze-up and break-up processes lead to an abrupt change of lake surface properties (e.g.albedo, roughness), affecting mass and energy exchange between the lake water and atmosphere.The lakes basically act as heat sinks in summer and heat sources in autumn and winter (Wilson, 1977).The relatively cool/warm air above the lake will inhibit/enhance convection in summer/winter.The variations in ice phenology affect the amount of on-lake evaporation, since the lake-ice blocks the water-air interface, and sublimation does not compensate this effect.A deep lake, such as Nam Co, can supply a large amount of heat energy to the atmosphere in the post-monsoon period (Haginoya et al., 2009).Air masses passing over a large, relatively warmer open-water area during winter can lead to heavy snow fall on the downwind shore.This socalled lake effect is well-known in the case of the large Lurentinan Great Lakes (Eichenlaub, 1970).For the TP, lake effect has been described in the case of Nam Co, which is a large saline lake in the central part of the plateau (Li et al., 2009), and has been analysed by remote sensing data by Kropáček et al. (2010).This exchange process is obviously interrupted by lake freeze-up when no more warmer openwater is available.The relationship between lake freeze-up and lake-effect-determined snow cover in Nam Co Basin has been compared on time series of optical satellite and passive microwave data by Kropáček et al. (2011a).

Objectives
So far, little attention has been paid to the ice regime of lakes on the TP in comparison with high latitude lakes in other regions.The comprehensive analysis of Northern Hemisphere lake ice cover as a response to changing temperature presented by both Magnuson et al. (2000) and Dibike et al. (2011) (based on ground observations and modelling approaches, respectively), do not include the TP.Analysis of the ice regime of two large lakes on the TP has been carried out so far: the Qinghai Lake on the NE margin of the TP by Che et al. (2009) using data from a passive microwave satellite instrument, and Nam Co in the central TP by Ye et al. (2011), using a combination of passive microwave and optical satellite data.An overview of lake ice phenology of at least the largest lakes on the TP is, to our knowledge, still missing.
Ground observations of lake ice phenology on the remote TP with harsh winter conditions are not available.Such costintensive observations are furthermore complicated in the case of large lakes, where the surface of the whole lake cannot be observed from a shore station (Walker and Davey, 1993).The number of such stations is globally decreasing, for instance in Canada, decreasing from 140 sites in mid-1970's to less than 20 sites by the end of the 1990's.Obviously the analysis of lake ice phenology on the TP has to rely on satellite observations.This study is focused on the derivation of the time series of lake ice phenology events from MODIS 8-day snow composites, for a representative group of lakes on the TP, allowing for the calculation of the duration of the ice cover period.The relation of ice phenology to (i) local factors and to (ii) climate variables is analysed in order to assess the potential of ice phenology as a proxy for climate variability.Additionally, trends in duration of the ice cover in the period from 2001 to 2010 for the studied lakes are derived.

Study area
The TP has a mean elevation exceeding 4500 m a.s.l. and it is bordered by the high mountain ranges Karakorum and the Himalayas in the south and by Kunlun Shan in the north, extending to a length of several hundreds of kilometers.The climate of the study region is under the influence of the westerlies in winter and the Asian monsoons in summer, and is characterized by wet summers and cold, dry winters.The most relevant climate parameters for this study are given in Fig. 1, which shows air temperature and wind speed/direction on the TP.The mean annual air temperature for 2001-2010 (Fig. 1a) features an evident dependency on orography, but also a negative gradient from south-east to north-west.The central and southern part of the TP have a mean annual temperature between 0 and 5 • C, dropping to far below 0 • C in the north-western part of the plateau.Only in the south, around Lhasa, does the mean annual temperature exceed 5 • C. The climate of the TP is typical with a pronounced air circulation pattern driven by both the westerlies and the monsoon.Figure 1b, c and d shows, respectively, the mean 10 m wind speed and direction for 2001-2010.Wind speed is also strongly controlled by orography, with significantly higher wind velocities on mountain ridges.The dominant surface wind direction is west-east for most parts of the TP, which is owing to the prevailing influence of the westerlies in winter, with high velocities and a constant flow pattern.During the summer months the picture is different, with lower wind velocity and a cyclonic circulation forming in the low levels of the atmosphere and at the surface.This well-known low-pressure system is linked to the summer plateau heating and associated air updraft (Gao et al., 1981).
Analysis of meteorological records over the last decades reveals that both temperature and precipitation show an increasing trend (Liu and Chen, 2000;Xu et al., 2008).This general climate trend appears to feature some regional differences across the TP.The south-eastern part of TP becomes warmer and wetter, the south-western part follows the same trend but with smaller amplitude, while the north-eastern part of the TP turns warmer and drier (Niu et al., 2004).
There are a large number of lakes on the TP differing in size, salinity, altitude and shape.Putting the limit of the minimum area to 0.1 km 2 , their number reaches almost 6880, covering a total area of around 43 000 km 2 , out of which 1260 lakes exceed 1 km 2 (Jiang et al., 2008).This makes the TP one of the largest lake systems in the world, unique in the high-elevation conditions.There are 7 lakes larger than 500 km 2 distributed mainly in the SW part of the plateau, out of which 4 exceed 1000 km 2 .So far, there is little known about the bathymetry of even the largest lakes (Wang et al., 2010).
Most of the lakes on the TP are brackish or saline.An overview of salinity in three classes is presented by the "Map of the Lakes and Glaciers on the Tibetan Plateau" (Yao, 2007), while zoning of Tibetan lakes according to their hydrochemical type is described in Zheng and Liu (2009).Additionally, several figures on salinity of the Tibetan lakes are dispersed in literature.Nevertheless, a comprehensive data set describing the salinity of lakes on the TP is not available.Salinity, in the context of limnology, represents the total sum of ions dissolved in water.Most of the lakes on the TP represent a termini of inland drainage basins and are therefore salty.Only a few lakes, e.g.Taro Co in the SW part of the plateau, are through-flow lakes containing freshwater (salinity < 3 g L −1 ).Lakes in highly arid basins with little surface inflow may become hypersaline, which is, for instance, the case of Chabyer Caka (Zabuye Lake) with salinity of its surface brine reaching up to 425 g L −1 (Zheng and Liu, 2009), or Lagkor Co in the SW part of the plateau.Average salinity of five hydrochemical zones defined by Zheng and Liu (2009) ranges from 135 to 352 g L −1 .Such high salinity indeed influences the freezing temperature of the lake water.The magnitude of the freezing point depression decreases with an increase of salinity.This in turn leads to a shortening of the ice cover period of salty lakes in comparison with freshwater lakes under the same conditions, i.e. climate, elevation, etc.
Levels and extent of lakes on the TP experienced high oscillations as a reaction to changes in climate conditions in the past (Zhang et al., 2011;Phan et al., 2012;Kropáček et al., 2011b).Furthermore, the lake shore may have shifted by even hundreds of meters during the last 40 yr (Bianduo et al., 2009;Kropáček et al., 2011b).These oscillations were related to an increased runoff, owing to the retreat of mountain glaciers, according to Yao et al. (2007), who analysed several catchments in Tibet.There is a lot of confusion in the terminology, caused by parallel usage of several versions of names for the same lake.Often more than two versions are used in Tibetan, while different transcriptions from Tibetan (often via Chinese) increase the confusion.Here, we decided to adhere to a convention introduced by the Map of the lakes and glaciers on the TP (Yao, 2007) compiled by the Institute of TP Research, Chinese Academy of Science (ITP, CAS).

MODIS instrument and data description
The study area features favourably low cloud cover during winter and spring periods, and the data acquisitions by optical instruments are not hindered by polar darkness, as is the case for instances of high-latitude Canadian lakes (Howell et al., 2009).Furthermore, the medium resolution optical data, which have spatial resolution in the order of several hundred meters, allow for derivation of ice phenology for tens of lakes on the TP.This is in contrast to passive microwave data, with spatial resolutions in the order of tens of kilometers, which are useful only for the largest lakes (Che et al., 2009).In this study, we used MODIS data for the derivation of lake ice phenology of Tibetan lakes.
The Moderate Resolution Imaging Spectroradiometer (MODIS) instrument is carried by two polar orbiting satellites, Terra and Aqua.The system is able to acquire data of the entire Earth's surface every 1 to 2 days in 36 spectral bands by both Terra (in the morning) and Aqua satellites (in the afternoon).In order to keep the data amounts at a reasonable level in this study we used the MODIS snow product MYD10A2 (Hall et al., 2007) instead of primary spectral bands.This product has a resolution of 500 m and is a composite classification image over an 8-day period, based on a Normalized Difference Snow Index (NDSI) and several other criteria tests that minimize the impact of cloud cover.The classification scheme matches the needs of this study as it contains the following five classes: Snow, No snow, Lake, Ice and Cloud.

Derivation of ice phenology dates from the MODIS 8-day composite data
The frozen lake surface is represented by classes Ice and Snow in the MODIS 8-day composites.We adopted the percentage of pixels classified as class water as a measure for the frozen status of the lake.This approach provides more accurate estimates than querying the number of pixels belonging to classes Ice and Snow, since it accounts better for missregistrations between MODIS scenes.For practical reasons, a threshold for the ice/water extent had to be set in order to detect the ice phenology dates.Factors 5 % and 95 % for the area of open water were chosen as thresholds for the identification of the ice phenology events in order to exclude the influence of noise.Such absolute thresholds would be sensitive to various disturbing effects.Instead we used an adaptive approach.The data set was first divided by a median into two parts: one below median corresponding to frozen lake, and the other one with the values above median corresponding to unfrozen conditions.In the next step, mean values were calculated from both data sets.The thresholds 5 % and 95 % were then applied to the range limited by these calculated mean values, instead of the whole range 0-100 %.This helped us to suppress the influence of mixed pixels and remaining discrepancies between the used lake outlines and the real shorelines.
For the extraction of ice phenology events from MODIS data, all lakes larger than 100 km 2 were selected.This resulted in 65 lakes that are well-distributed over the TP (Fig. 2).Out of this group, 6 lakes (Yinmar, Yamzhog Yumco, Nau Co, Chabyer Caka and Ringinyububu Co) had to be sorted out, since their histograms were too noisy for a reliable extraction of ice phenology.For instance, the Yamzhog Yumco has a highly complex tentacle-like shape that results in an extremely noisy signal.In order to account for regional differences in the ice regime of lakes in such a large area, the lakes were grouped by k-means clustering in a multidimensional space, defined by six variables (area, elevation, latitude, longitude, shape index, and mean temperature) into three groups denoted by abbreviations A, B and C respectively (see Fig. 2 and Table 1).These groups correspond to relatively geographically compact regions located in the north-west, north-east and south, respectively.This approach aims to select homogenous groups in terms of factors affect-Table 1.Average values of basic characteristics and ice phenology for each lake group.Duration of complete ice cover (DCI) and duration of ice cover (DI) in days calculated as a mean value for each of the three groups over the period from 2001 to 2010 and standard deviations (σ ) of both variables, which is a measure of variation of ice cover duration amongst the lakes in each group.In order to get a regional overview of the ice regime on the TP, mean values of DI and DCI, phenology dates and trends for these groups were calculated and compared.
For definition of the spatial extents of the lakes we used Global Lakes and Wetlands Database -GLWD (Lehner and Döll, 2004), which is a global data set created as a combination of various available sources.The shorelines from the GLWD were checked against the MODIS images.In some instances, a certain amount of manual editing of the shoreline had to be involved, i.e. shifts or shape adjustments.In order to account for small variations of the lake shapes caused by mixed pixels, a mean image was calculated from MODIS images covering the autumn period after the monsoon and before the lake freeze-up.This image, which shows in one layer spatial variations of shape between single acquisitions, was used as a base layer for adjustments of the vector data.Another source of the size variation is probably a shoreline displacement connected to the lake volume/size oscillations (Bianduo et al., 2009;Wu and Zhu, 2008;Kropáček et al., 2011b) during the study period.These variations usually do not exceed the pixel size of MODIS.
In the next step, histograms were extracted from all MODIS images in the time series separately for each lake.Graphs of the percentage of open water area were constructed for the whole time period.The ice phenology events were then identified automatically as intersections of the curve representing open water with thresholds 5 % and 95 %.It was observed that multiple intersections with thresholds occur during the freeze-up and break-up period.This can occur owing to ice retreat caused by wind events and refreezing during the break-up period, or by cloud cover not accounted for by the 8-day compositing approach.Another reason is the presence of noise in the time series induced by classification uncertainty.The magnitude of the noise is given by variation of the water curve during the summer period, when it is not due to cloud cover.
In order to identify correctly the FU in case of multiple intersections of the open water curve with the threshold, which is mainly owing to noise, a point that has been set well-before possible early freeze up was used.The first threshold crossing was searched for from this point.Analogously, all the other ice phenology dates were identified (Fig. 4).Since the points of the open water have 8-day sampling, the exact day of each ice phenology date was obtained as a linear interpolation of the first point below and above the threshold.All identified points representing the ice phenology dates were plotted together with the open water curve for a visual check.

Validation of open water area derived from MODIS 8-days composites
The accuracy of the derivation of the open water area for a lake from the MODIS 8-day composites was estimated by a comparison with a number of reference satellite images (Fig. 5).These images cover various phases of freeze-up and break-up processes and are from various satellite sensors, including Envisat/ASAR (Advanced Synthetic Aperture Radar, wide swath mode data), Landsat, TerraSAR-X (Synthetic Aperture Radar) and SPOT (System for Earth Observation).For each validation image, a percentage of the open water area was derived and compared with the water curve from MODIS 8-day data, which was interpolated for each day.A time difference in days between each point representing the validation image and a corresponding point on the interpolated curve with the same percentage of open water was derived.An overall accuracy was then calculated as a root mean square error of these time differences in days.An RMS error equal to 9.6 days was obtained from twenty validation images for Nam Co.A mean time difference as low as 1.2 days proves that there is no systematic error in the open water area estimation.

Correlation of the lake phenology with the local and climate factors
In an attempt to identify the main drivers of the ice regime of the study area, the extracted ice phenology data was compared with data describing local and climate factors.The local factors, which were possible to describe by available data, were: latitude, longitude, elevation, area and shape.Latitude and longitude were retrieved from the GLWD using automatically generated centroids of lake polygons.The centroids were manually shifted inside the lake in several cases where they fell outside the lake polygons, owing to the concave shape of the shore.Lake altitude was identified in the Shuttle Radar Topography Mission (SRTM) digital elevation model using the edited centroids.A simple shape index was defined as a ratio between perimeter and area, which can be easily The identified lake ice phenology events are marked with symbols: "x" for FO, "+" for FU, "*" for BU and "♦" for WCI.The noise in the total ice cover is often caused by confusion with clouds, which are actually rare in the winter period on the TP.retrieved from the GIS layer of lakes.This index describes shape complexity of the shore.Lakes with a complex shape are likely to have a lower depth and, as a consequence, a lower thermal inertia in comparison with a rounded lake of the same size.Lakes with high shape index are likely to provide more shallow bays where the freeze-up usually starts.In addition, the parameter lake area itself is a proxy of the lake volume and is related to the thermal inertia of the lake.For characterization of the climate for the selected lakes it is not possible to simply rely on measurements of nearby ground meteorological stations, because of their sparseness.The density of stations decreases towards the west, leaving the whole western part of the plateau with some five stations.To our knowledge, only Nam Co and Gyaring Co have a ground meteorological station located in the immediate vicinity of its lake shore.Instead we rely on a climate data set generated with the Advanced Weather Research and Forecasting (WRF) numerical atmospheric model (Skamarock and Klemp, 2008).Two-way grid nesting in the parent domain (30 km spatial resolution) yields to a regionalscale spatial resolution over the TP of 10 km (Fig. 1).All details of the model configuration are presented in Maussion et al. (2011).The 10 yr data set (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) is composed by consecutive model runs, reinitialized every day using data from the Global Forecast System's final analysis.The model is strongly constrained by the observed synoptic weather patterns and thus represents a regional atmospheric reanalysis.
Other studies successfully made use of this modelling technique for other regions of comparable characteristics (complex terrain, scarce observations), for example Iceland (Bromwich et al., 2005) or Greenland (Box et al., 2006).Producing a high-resolution reanalysis data set is part of a larger project (Maussion et al., 2012).For the present study we use time series of daily and monthly averages of hourly values for two variables (2 m air temperature and 10 m wind speed).A validation effort using 84 weather stations located within the 10 km domain, showed a good match between the measured and the simulated daily mean temperature (mean bias of −0.91 K, Root Mean Square Deviation (RMSD) of 2.74 K and coefficient of determination (r 2 ) of 0.92).Wind speed was validated for monthly means owing to its high variability.Direct comparison with observations is difficult because of the altitude difference between station and model grid points, but the inter-annual and annual variability of wind speed are well-reproduced by the model with an r 2 of 0.44.Contours of mean temperature, wind speed and wind direction on the TP are shown in Fig. 1.Based on the reanalysed data, the mean annual temperature was identified for each lake.In addition we derived a thawing index as the cumulative temperatures above 0 • C, and the freezing index as below 0 • C days.These two variables characterize the temperature in periods that affect ice thawing and ice melting respectively.Correlation coefficients were calculated between variables describing ice phenology and those describing the local and climate factors (see Table 3).

Results and discussion
The lake ice phenology for 59 lakes on the TP has been extracted from the time series of MODIS 8-day composites for the period from 2001 to 2010.This allowed us to get an overview of the regional differences in ice regime over the TP (Fig. 6).There is a remarkable difference in the ice cover duration in terms of both DCI and DI amongst the studied lake groups (Table 1).The groups resulting from k-means clustering form well-defined geographical regions (Fig. 1).Group C in the south features the lowest values for both DCI and DI, while group A, containing mainly lakes in the north-west of the plateau, has the highest values for DCI and DI.The ice cover duration for group A is roughly two-times greater than group C. Owing to these high differences in ice cover duration amongst the groups, any general figures for the whole TP would not be representative.Statistical significance of the differences in DCI and DI amongst all three groups was checked by ANOVA (Analysis of variance between groups) test (F-value = 13.8,P-value = 0.001 for DCI and F-value = 9.2, P-value = 0.001 for DI).There is a high variation in both the DCI and DI in each group (Table 1).Grouping of the lakes into a higher number of clusters (four and five) using K-means clustering resulted in spatially heterogeneous regions and moreover did not lead to a decreased variation in ice phenology variables in the groups.
For several lakes it was not possible to extract the ice phenology because of a high level of noise in the signal obtained from the MODIS 8-day composites.This was most likely to be owing/attributed to the effect of mixed pixels.This effect impacts to a greater degree smaller lakes and lakes with complex shape.Since the lower limit of the used method, in terms of lake size, is shape dependant, it would be difficult to set a threshold of certain minimal lake size as a clear limitation of the method.In the case of some lakes, open water was misclassified as class land in the MODIS 8-day composites, which can be probably explained by a high sediment load in the water, for instance in the Ulan Ul Lake.

Validation of the no complete freeze-up events
It appears that there was no complete freeze-up in some years for several lakes.This may be surprising in an area with a harsh climate and with a long period of temperatures below 0 • C.This occurs also in the case of some large lakes like, for instance, Taro Co.In order to exclude the possibility that the phenomena is only random noise in the data, quick-looks of Landsat scenes were checked for the presence of open water in the relevant periods.The presence of open water during usual ice-on season was confirmed in this way for a number of lakes (see Table 2).

Group A
The main characteristics of this group are high elevation, low mean annual temperature, and relatively small size.This cluster contains only twelve lakes, which are located mainly in the very western part of the TP bordered from the north by the Kunlun Mts. (Fig. 2, Appendix A).Three lakes are located further in the east, somewhat separated from the main group surrounded by lakes of the group B. The reason why these lakes belong to the group A instead of the group B is a low mean-annual temperature (Lixi'Oidaim Co: −7.4 • C, Taiyang: −7.6 • C and Hoh Xil: −8.3 • C) when comparing to the mean value for the group B: −5.4 • C. The ice regime of lakes in this group is determined by cold continental climate that is reflected by a long ice cover period: DCI 149 and DI 209 days (Table 1).Longmu Co has a remarkably short duration of ice cover, both in terms of DI and DCI, and in some seasons it did not freeze-up completely, while the longest ice The Cryosphere, 7, 287-301, 2013 www.the-cryosphere.net/7/287/2013/cover was observed by the Hoh Xil Lake.The anomaly of Longmu Co is not easy to explain since it is one of the highest lakes of the group and its mean temperature is not exceptionally high.The mean trend in the ice cover duration is 2.6 and 1.4 days yr −1 for DCI and DI, respectively (Table 3).Nevertheless, there is a high variation of the trend in this group (Appendix B).The mean trend in DCI is also strongly influenced by the high trend of the two lakes Aksayqin and North Heishi, with extremely low values of DCI in 2001 and 2003.

Group B
This group contains in total eighteen lakes, out of which five are larger than 300 km 2 (Fig. 2, Appendix A).The lakes are mostly located in the northern part of the plateau.This group includes also three lakes in the west and another three lakes in the south-west surrounded by group A and C, respectively.This inclusion is due to their high elevation and low mean annual temperature.On average, lakes in this group feature a relatively long ice cover period (DCI 109 days and DI 159 days).There is a high variation of both DCI and DI (Table 1), which is influenced by a large variation of elevation (Appendix B) and a large range of latitude.Lake Ayakkuh features remarkably short mean DI of 88.3 days in comparison with the average DI of the group, mainly due to its low elevation (4161 m a.s.l.).Lakes Xijir Ulan, Yaggain Canco and Dogai Coring did not freeze up completely during the seasons 2001 and 2002.The mean trend in duration of lake ice period is negative in terms of DCI, but positive in terms of DI (Table 3), and amounts for −0.2 and 1.7 days yr −1 , respectively; however, there is again a high variation in both the DI and DCI in this group (Table 3).Four lakes in this region appear to have a rapid increase in the ice cover duration in terms of the DI: Dogai Coring, Xijir Ulan, Ayakkuh Lake and Yaggain Canco.All of these lakes have also a high variation in ice cover duration (Appendix B).In the case of Dogai Coring, Xijir Ulan and Yaggain Canco, both the increasing trend and the high variation are strongly impacted by the occurrence of no freeze-up events in the studied period.

Group C
Lakes of this group, which can be characterized by a low elevation and a high mean annual temperature (Table 1), are distributed approximately along the parallel 31 • N in the south of the Tibetan Plateau (Fig. 2).This cluster contains twenty nine lakes, which is the highest number of all the three groups (Appendix C), including many of the largest lakes on the TP, e.g.Serling Co, Nam Co, Zhari Namco, Ngangla Ringco, Tangra Yumco and Taro Co. Except for Mapam Yumco, La'nga Co, Paiku Co and Puma Yumco, all lakes are in the central part of the plateau north from the Gangdise and Nyainqentanglha Mts., which act as a barrier to the monsoon circulation.The lakes Bangong Co and Geyze Caka are located in the western part of the plateau.The lake ice cover in this region, with an average DCI of 72 days and a DI of 126 days, is the shortest of the three regions.A number www.the-cryosphere.net/7/287/2013/The Cryosphere, 7, 287-301, 2013 of lakes in this group, e.g.Taro Co, Tangra Yumco, Paiku Co, Xuru Co, Monco Bunyi, Dajia Co, Bam Co, Npen Co in and Ringco Ogma, did not freeze up completely in some years.Gyeze Caka did not freeze-up completely during the whole period from 2001 to 2010, and there are no signs of it freezing up before, which could be due to its hypersalinity ("caka" means "salty lake" in Tibetan).Another interesting lake in this group is Tangra Yumco, which is regarded as the deepest lake on the TP (Wang et al., 2010).It is divided into a northern and a southern part by a narrow strait.The northern part freezes over later or not at all, as seen on Landsat images (for instance on 31 March 2007, 27 February 2008and 15 March 2010).This is obviously caused by heat accumulation in the northern part that reaches a depth of 214 m, which is approximately 100 m deeper than the southern part (Wang et al., 2010).In the case of Taro Co, the influence of high salinity on this effect can be excluded since this lake, while having an outlet towards Chabyer Caka, is a fresh water lake (Wang et al., 2002).The DI of lakes in this region varies from 59.2 days (Paiku Co) to 204.8 days (Pangkog Co).In the case of Paiku Co it was only in 2008 that it froze over completely.Pangkog Co, which is located immediately to the east of Serling Co, is, with an altitude of 4529 m a.s.l., the third lowest lake in the group.A possible explanation of its long ice cover would be its relatively small size of 100.7 km 2 combined with its low depth.Mean trend in ice cover duration of this group is −1.1 and −1.6 days yr −1 in DCI and DI, respectively, with a relatively low variation (Table 3).

Trends in duration of ice cover for the three groups
Trends in DCI and DI for all regions were calculated as mean values of trends in each group.Statistical significance of trends in both DCI and DI are given in Appendices A, B and C. The mean values for each group are listed together with the corresponding standard deviations in Table 3.The statistical significance test by means of ANOVA (F-value = 3.45, P-value = 0.038 for trend in DCI and Fvalue = 0.05, P-value = 0.95 for DI) revealed that the groups do not differ in terms of mean trend.The standard deviations reveal the magnitude of differences in the behaviour of lakes in each group.The especially high standard deviation in the case of group B suggests that there is a strong influence from the local factors on the ice regime of these lakes.There are differences between the trend in DCI and the trend in DI for the same region, which reflects variation in duration of the freeze-up and break-up periods and it is also highly influenced by no complete freeze-up events.The differences in duration of the freeze-up and break-up periods may be caused by snow fall and wind events, insulating effect of snow cover slowing up the break-up, changes in albedo, etc.It has to be noted that the variations of ice cover duration of all groups exceed the calculated trend (Table 3), meaning that even if the lakes were grouped with respect to the factors influencing their ice regime, there are large differences in the behaviour of single lakes in each group.It is not easy to identify the reason for this; however, this may be owing to some local factors, for instance a feedback effect of snow cover.

Correlation of ice phenology variables with climatic and local parameters
Statistical analysis of the dependency of ice phenology variables on climatic and local parameters was carried out for all 59 studied lakes (Table 4).In order to assess the sensitivity of the analysis to lake selection and to suppress the influence of noise introduced by mixed pixels, which obviously play a more important role in the case of smaller lakes, the same analysis was repeated for lakes larger than 300 km 2 (18 lakes in total, Table 5).Results of the analysis showed that there is a high correlation between ice cover duration (DCI and DI) and temperature and both the thawing and freezing indices.The correlation coefficient, which is 0.62 to 0.63 for all lakes, increases up to 0.82 when selecting lakes larger than 300 km 2 .This confirms a high thermal determination of lake phenology reported by many studies for a number of lakes worldwide (Ruosteenoja, 1986;Manguson et al., 2000;Latifovic and Puoliot, 2007).It is likely that the dependency of ice phenology on thermal variables would be even higher if the freezing point depression caused by salinity could have been taken into account.It appears that the correlation with the three variables based on temperature is higher for FO and WCI than for FU and BU.This suggests that the FO and WCI are more thermally determined.Surprisingly there is only a little difference in dependency of freeze-up and breakup events on the freezing and thawing indices.It appears that for the large lakes there is a high dependency of ice phenology on the shape index.The strikingly lower value for all lakes, including many smaller lakes around 100 km 2 , is likely to be due to the generalized lake outlines in the lake database, which affects the reliability of the calculated shape index.High correlation of ice cover duration (r = 0.46 for DCI and r = 0.45 for DI of the group of the large lakes) with latitude is in fact related to a decrease in temperature and irradiation towards the north.There is probably no strong driver While there is no significant correlation of ice phenology variables with the perimeter/area index for all lakes, there is a strong correlation for the selected large lakes achieving r = 0.79 in the case of duration of complete ice cover (DCI).This shows that the shape complexity significantly affects the ice phenology, in particular, the freeze-up process providing shallow bays for the development of early ice.On the contrary, the break-up process seems to be largely unaffected by the shape complexity, even in the case of the large lakes.For the parameter area, a statistically significant correlation was found only in the case of DCI and FO (for the 18 largest lakes).This suggests that the studied lakes with a large area do not always have a large volume and thermal capacity.
The analysis revealed much higher dependency of FO and WCI to the climate and local parameters than in the cases of FU and BU.This can be caused by two factors.Either wind or the local conditions can lead to a higher variation in FU and BU, but do not affect the FO and WCI in the same way; alternatively, the spatial variation of BU and FU is purely caused by systematically higher inaccuracy in BU and FU estimation.Irrespective of the reason, the FO and WCI (respective DI as their difference) derived from MODIS 8-day snow composites appears to be a better indicator for changes in climate on the TP.While comparing the correlation of BU and FU with temperature, it is the FU that appears to be more thermally determined, which contradicts findings of Duguay et al. (2006), who found the BU to be a more robust indicator of climate oscillations than FU in the case of Canadian lakes.It seems that apart from temperature, there are some impor-tant factors that influence the break-up on the TP.Wind disturbances during the freeze-up period, marked by the maximum wind speed for all three groups in February (Fig. 3), are thus likely to have a lower impact on the freeze-up process than the insulating effect of snow cover on the break-up process in the case of Tibetan lakes.

Conclusions
Ice phenology was derived for 59 lakes on the TP from MODIS 8-day composite data.The obtained results imply that despite some influence of local settings, e.g.salinity, shape or wind exposure, the ice phenology of the Tibetan lakes is to a high degree determined by climate, especially the air temperature, which is in agreement with results from other cold regions in the world.This suggests that changes in the ice phenology of Tibetan lakes can indicate variations in regional climate.This is highly valuable since the area is only sparsely covered by ground observations, and at the same time plays an important role in regional and global climate.The ice phenology of the studied lakes, averaged over the acquisition period, features a high variation both in the lake groups and when comparing them, which reflects differences in geographical position (altitude and latitude) and local settings, e.g.shape complexity and salinity.This study showed that a number of Tibetan lakes do not freeze up completely during some winter seasons.This interesting phenomenon and its consequences for local climate, including snow cover regime, should have proper attention paid to it in further research.Both of the variables, used as a measure of duration of the ice cover, DI and DCI, can be used as climate indicators, although DI appeared to be more stable in terms of variation in two of the three regions than DCI.Additionally, the estimation of DCI is hindered when the lakes do not freeze up completely.It can be thus concluded that DI allows for a more reliable characterization of the ice regime from MODIS 8-days composites for lakes on the TP.Furthermore, it appears that FO is more thermally determined than BU in the case of the studied lakes.
The available time series of MODIS data allowed us to derive a trend in ice cover duration for the studied lakes, which appeared to be highly variable over the TP.The high variation remains even after grouping the lakes by clustering in a multi-dimensional space defined by relevant climatic and local variables.These facts, together with the relatively short evaluation period of ten years, would make any general conclusion about a trend in ice phenology on the TP dubious.Even though this study confirmed the usefulness of lake ice phenology as a climate indicator, the extraction of a significant trend for the assessment of changes in climate on the TP requires a longer time series of observations.The high variation of the trend amongst the lakes also points out the fact that whenever ice phenology will be used as a proxy of climate change on the TP, proper care has to be paid to the selection of a representative sample of lakes.

Fig. 1 .
Fig. 1.Climate properties of the Tibetan Plateau for the period 2001-2010: (a) annual mean 2 m temperature, (b) annual mean 10 m wind speed, (c) mean 10 m wind speed for December, January and February, (d) mean 10 m wind speed for June, July and August.Data points below 3000 m a.s.l.have been masked out for clarity.Flow arrows are shown for every 7-th model grid point and indicate average wind direction and wind speed.See Sect.3.4 for details on the reanalysed climate data set.

Fig. 3 .
Fig. 3. Mean monthly values of 2 m air temperature and 10 m wind speed averaged for each lake group based on reanalysed climate data (see Sect. 3.4 for details) for the period 2001-2010.Both variables have a marked seasonality.

Fig. 4 .
Fig. 4.An example of curves for open water (below), ice and snow (middle) and clouds (above), derived from MODIS 8-days composite data for Dorsoidong Co.The identified lake ice phenology events are marked with symbols: "x" for FO, "+" for FU, "*" for BU and "♦" for WCI.The noise in the total ice cover is often caused by confusion with clouds, which are actually rare in the winter period on the TP.

Fig. 5 .
Fig. 5. Validation of the time series representing the area of open water for Nam Co, derived from MODIS 8-day composite data by comparing the open water curve with points representing the area of open water extracted from high-resolution satellite images (Landsat, Envisat/ASAR, TerraSAR-X and SPOT).

Fig. 6 .
Fig. 6.Ice cover duration for the studied lakes on the TP showing duration of freeze-up period, duration of complete ice cover (DCI) and duration of break-up period as heights of the respectively coloured sections of bar symbols, whereas the total size of the bar represents the duration of ice cover (DI).

Table 2 .
Validation of the no complete freeze-up events detected by MODIS data for selected lakes by Landsat data.Percentage of ice cover was determined from available Landsat images falling to the period of usual ice cover determined for the particular region.

Table 3 .
Trends in days yr −1 in duration of both DCI and freezing period DI in the period from 2001 to 2010 for each group.The corresponding standard deviations of trends gives an idea about the variation of trends in each group.

Table 4 .
Correlation of ice cover duration (DCI -duration of complete ice cover and DI -duration of ice cover) and ice phenology dates with climate and local parameters calculated for all 59 studied lakes.Statistically significant values at the level of significance 0.02 are bold.

Table 5 .
Correlation of ice cover duration (DCI -duration of complete ice cover and DI -duration of ice cover) and ice phenology dates with climate and local parameters calculated for lakes larger than 300 km 2 (18 largest lakes).Statistically significant values at the level of significance 0.02 are bold. of the ice regime, which would lead to an east-west gradient in ice phenology, since only a weak correlation was revealed between the ice phenology variables and longitude.The ice cover duration appears to change with altitude, owing to a decrease in temperature connected to a decrease in air pressure, which can be expressed by temperature lapse rate.The selection of the large lakes leads to a significant increase in dependency of ice cover duration to altitude.This could have been likely caused by an inaccuracy in determination of lake ice phenology dates for smaller lakes with noisy time series of the satellite data classifications.