Use of Landsat Satellite Images in the Assessment of the Variability in Ice Cover on Polish Lakes

: Despite several decades of observations of ice cover in Polish lakes, researchers have not broadly applied satellite images to date. This paper presents a temporal and spatial analysis of the variability in the occurrence of ice cover on lakes in the Drawskie Lakeland in the hydrological years 1984–2022 based on satellite data from Landsat missions 4, 5, 7, 8


Introduction
Elements of the hydrosphere of moderate latitudes undergo seasonal changes determined by climatic conditions. Parameters of rivers and lakes change broadly from high water temperatures in the summer season to a change in the physical state in the winter season, where the appearance of ice cover completely changes their functioning in comparison to other seasons. Due to this, the issue of ice cover is a frequently discussed subject that refers to both the physicochemical and biological conditions. Ice cover shapes the distribution of dissolved oxygen [1,2], light dispersion [3,4], metabolic dynamics [5], and heat content [6]. Research conducted on Lake Ulansuhai (China) showed that the presence of ice in the winter season may cause a rapid increase in nutrient concentration in water and sediments [7]. Moreover, a decline in ice cover favours wave growth in winter, as evidenced by the analysis of Lake Michigan [8]. The presence of ice affects the composition and activity of plankton [9,10] and waterfowl [11].
In the case of Poland, more than 7000 natural lakes exist across the country [12], constituting an important element of the natural environment, particularly in the north of Poland where they are the most abundant [13]. Numerous analyses of ice cover on lakes have been conducted to date [14][15][16][17][18][19], but the information is still limited to a relatively small population (several tens of lakes) under systematic monitoring. They are lakes with surface areas from 100 to 11,000 ha. They are therefore relatively large for Poland, where the average surface area of all lakes is 39.6 ha. [12]. This background lacks a reference to information on the remaining lakes, those with a smaller surface area that are equally important in environmental and economic terms. Importantly, study results for 220 lakes in the northern hemisphere showed that despite common warming patterns, individual lake parameters may determine the rate of loss of ice cover, and this knowledge is necessary for future activities in the scope of environmental modelling and protection [20]. Expanding knowledge in the scope of ice cover on lakes is possible owing to the application of remote sensing methods that have been dynamically developing at least for the last several decades.
Technical progress supports an increase in the accuracy and possibilities of satellite methods as well as the use of various combinations of data. Heinilä et al. [21] applied the method of accessing data regarding ice cover range in lakes based on multidimensional Gaussian distributions calculated for training data with the application of several reflectance and thermal bands and related indices. In the case of large lakes in the northern hemisphere, Wu et al. [22] assessed the ability of four machine learning classifiers to map ice cover of lakes, water, and cloudiness during periods of ice disintegration and freezing. It was determined that all classifiers using a combination of seven bands proved able to generate total classification accuracy above 94%. Wang et al. [23] investigated the efficiency of the semiautomatic segmentation "glocal" Iterative Region Growing with Semantics (IRGS) algorithm for the classification of lake ice by means of dual polarity RADARSAT-2 images obtained over Lake Erie. The analysis of Lake Erie showed that the "glocal" IRGS algorithm can provide credible ice-water classification by means of dually polarised images featuring high total accuracy. Satellite images permitted monitoring surface water/ice dynamics both in space and in time, although using them involves many challenges [24]. According to Herrick et al. [25] in the case of teledetection data, the dataset is still mainly available to experts in the field, and requires advanced skills. The method proposed herein applies a relatively simple approach illustrating the so-far-unavailable information regarding ice cover in Poland, based on remote sensing data. This offers new possibilities of investigation of lake ecosystems in reference to the phases of appearance and disappearance of ice, zones of its presence, and probability of occurrence. This information is key for fields such as limnology, hydrobiology, and ecology.
The study objective is the determination of the variability in occurrence of ice cover on selected lakes in Poland, and the determination of the effect of selected morphometric parameters on the occurrence of ice.

Study Area
The study area covered the Drawskie Lakeland in the northwestern part of Poland. The area features deep postglacial channels containing numerous lakes. The Drawskie Lakeland is enormously rich in terms of occurrence of lakes and other small water bodies. It includes approximately 500 of them (with a surface area of more than 1 ha). For 67 lakes, it was possible to determine detailed morphometric parameters, including: volume, maximum and mean depth, length, effective length, width, shoreline length and development, and exposure index ( Figure 1; Table S1) [26]. These lakes served as a further subject of the analysis of occurrence of ice cover. The largest lakes in the study area are Drawsko (18.1 km 2 ), Wielimie (16.9 km 2 ), Lubie (14.4 km 2 ), Siecino (7.2 km 2 ), Wierzchowo (7.0 km 2 ), and Komorze (3.9 km 2 ). The largest rivers in the Drawskie Lakeland include Drawa, Gwda, Parsęta, and Rega. Together with their smaller tributaries, they constitute a well-developed hydrographic system. Komorze (3.9 km 2 ). The largest rivers in the Drawskie Lakeland include Drawa, Gwda, Parsęta, and Rega. Together with their smaller tributaries, they constitute a well-developed hydrographic system.

Materials
This study used the USGS Earth Resources Observation and Science (EROS) Archive-Landsat Archives dataset: Landsat 4-5 TM, Landsat 7 ETM+, Landsat 8 OLI, and Landsat 9OLI + Collection 2 Level-1 Tier 1 and Tier 2. The compiled dataset refers to the cool half-year (November-April) and covers the period 1984-2022. The Drawskie Lakeland is located on path 191 and 192 on the interface of lines 22 and 23. The preliminary analysis covered scenes with cloud covers of up to 70%. Then, each scene was analysed in detail in terms of occurrence of cloud cover within the studied lakes. For the purpose of the identification of clouds and cloud shade, the CFMask algorithm [27] was applied. It was assumed that if no cloud cover occurred within the range of Lakes Lubie, Drawsko, Siecino, and Wielimie (the largest of the analysed lakes), the scene was qualified for further research procedures. As a result, a total of 181 Landsat scenes were selected for the analysis, including 152 scenes for Tier 1 and 29 for Tier 2 products. Each Landsat scene involves processing to one of the three levels: L1TP (Precision and Terrain Correction), L1GT (Systematic Terrain Correction), or L1GS (Systematic Correction) (h ps://www.usgs.gov/landsat-missions/landsat-levels-processing, accessed on 10 March 2023). Furthermore, Landsat 7 scenes had data missing as of June 2003 due to a failure of the Scan Line Corrector (SLC). The strip lines were gap-filled using Landsat toolbox's 'fix Landsat 7 scan line error' of ArcGIS software [28]. This study used 160 scenes processed to L1TP, and 21 scenes processed to L1GS.
This study also used the morphometric parameters of the lakes presented in the Atlas of Polish Lakes, which are: volume, maximum depth, mean depth, length, effective length, width, shoreline length, shoreline development, and exposure index. The data for individual lakes are presented in Table S1.

Materials
This study used the USGS Earth Resources Observation and Science (EROS) Archive-Landsat Archives dataset: Landsat 4-5 TM, Landsat 7 ETM+, Landsat 8 OLI, and Landsat 9OLI + Collection 2 Level-1 Tier 1 and Tier 2. The compiled dataset refers to the cool halfyear (November-April) and covers the period 1984-2022. The Drawskie Lakeland is located on path 191 and 192 on the interface of lines 22 and 23. The preliminary analysis covered scenes with cloud covers of up to 70%. Then, each scene was analysed in detail in terms of occurrence of cloud cover within the studied lakes. For the purpose of the identification of clouds and cloud shade, the CFMask algorithm [27] was applied. It was assumed that if no cloud cover occurred within the range of Lakes Lubie, Drawsko, Siecino, and Wielimie (the largest of the analysed lakes), the scene was qualified for further research procedures. As a result, a total of 181 Landsat scenes were selected for the analysis, including 152 scenes for Tier 1 and 29 for Tier 2 products. Each Landsat scene involves processing to one of the three levels: L1TP (Precision and Terrain Correction), L1GT (Systematic Terrain Correction), or L1GS (Systematic Correction) (https://www.usgs.gov/landsat-missions/landsat-levelsprocessing, accessed on 10 March 2023). Furthermore, Landsat 7 scenes had data missing as of June 2003 due to a failure of the Scan Line Corrector (SLC). The strip lines were gap-filled using Landsat toolbox's 'fix Landsat 7 scan line error' of ArcGIS software [28]. This study used 160 scenes processed to L1TP, and 21 scenes processed to L1GS.
This study also used the morphometric parameters of the lakes presented in the Atlas of Polish Lakes, which are: volume, maximum depth, mean depth, length, effective length, width, shoreline length, shoreline development, and exposure index. The data for individual lakes are presented in Table S1.

Methods
At the first stage, for all available scenes, top of atmosphere (TOA) reflectance ρ λ was calculated using the following equations [29]: where: ρ λ -TOA planetary spectral reflectance, without correction for solar angle (-); M-reflectance multiplicative scaling factor for the band (-); Q cal -level 1 pixel value in DN (-); A ρ -reflectance additive scaling factor for the band.
The values of M and A ρ were obtained from the metadata (MTL.txt) file distributed with the Landsat Collection 2 Level-1 data products. The ρ λ ' value is not the real TOA reflectance, because it does not contain a correction for the solar elevation angle. Therefore, in the second step, the solar elevation angle was taken into account, and the real TOA reflectance ρ λ was calculated using the following equation: where: ρ λ -real TOA planetary reflectance (-); θ SE -local solar zenith angle.
The purpose of calculating the real TOA reflectance is to standardise the input data for the subsequent stages of the analysis. Hence, the albedo range is usually between 0 and 1, which is unitless and represents the ratio of solar radiation that is reflected [30].
This procedure completed the preparation of data of the L1TP level. Data from the L1GS level required additional assigning of georeference due to the fact that the available scenes were not shifted in reference to the actual location. For this purpose, 1st-order polynomial transformation available in ArcGIS Pro software [ESRI, Inc., Redlands, CA, USA] was used. Each of the 181 available scenes was subject to the analysis of ice cover on lakes. This permitted assigning one of three classes to each scene: 0-no ice cover on any of the lakes; 1-partial ice cover, 2-full ice cover on all lakes.
The starting point of identification of the lakes' spatial extent was the 2022 Hydrographic Division Map of Poland at a scale of 1:10,000, made available by the National Water Holding Polish Waters (Państwowe Gospodarstwo Wodne-Wody Polskie). The data were accessed through Hydroportal (https://wody.isok.gov.pl/imap_kzgw/?gpmap=gpSIGW, accessed on 6 February 2023). Because the analysis covers 38 years, the water surface area could have been changing throughout that period due to the development of water vegetation within the range of the littoral zone or short-or long-term water level fluctuations in the lakes due to climate change. In this regard, the spatial extent of the lakes was verified each time. This was conducted using the Normalised Difference Water Index (NDWI) calculated based on the following formula [31]: where: ρ λ(green) -green spectral band; ρ λ(N IR) -near-infrared spectral band.
The NDWI values are within the range from −1 to 1, and for water objects, it adopts values higher than 0.2. The next stage involved the detection of the occurrence of ice cover within the range of particular lakes in the Drawskie Lakeland. The analysis of the lake ice cover employed a method combining the Normalised Difference Snow Index (NDSI) and blue band proposed by Cáceres et al. [30]. The detection of ice cover was based on NDSI calculated based on the following formula [32]: where: ρ λ(green) -green spectral band; ρ λ(SW IR1) -short-wave infrared spectral band.
NDSI can adopt values in a range from −1 to 1, whereas the probability of snow occurrence is proportional to how close the value is to 1. Moreover, the detection of occurrence of ice cover involved the application of TOA reflectance values for the blue band ρ λ(blue) . The following assumption was adopted at the stage of ice cover delimitation: where: During the determination of the range of ice cover, an NDSI value of 0.4 (a) was adopted as the threshold value. The analysis of ice cover occurring within the analysed lakes showed that it was impossible to determine universal threshold values b for the entire analysed period, or even for a single full winter season. This is related to the occurrence of ice cover in different phases (1-phase of development of ice cover on lakes, 2-intermelting phase with periodically occurring water layer on top of ice cover, or phase of disintegration of ice cover) and co-occurrence of snow cover. Due to the above, the detailed analysis covered two test lakes with the most frequently occurring partial ice cover, namely, lakes Drawsko and Siecino. They were subject to manual mapping of water surface and ice cover. The resulting learning dataset permitted the determination of threshold values b-Lake Drawsko, and the validation dataset allowed for the verification of the previously adopted values-Lake Siecino. Based on Lake Drawsko, the values of ρ λ(blue) were determined for the ice cover zone. The percentile 1% ρ λ(blue) value was used as the threshold value. Subsequently, this value was used to determine the extent of ice cover in Lake Siecino. In subsequent steps, percentile values of 2 and 3% were tested. The analysis was completed when a consistency was obtained between the ice cover determined manually and that based on percentile thresholds of 98%. The threshold value of parameter b during subsequent analyses of ice cover was in a range from 0.033 to 0.120. The analysis resulted in a dataset with data gaps for smaller lakes due to cloud cover. The dataset was used to determine a relationship between ice cover in individual lakes (the analysis was performed for each lake). Finally, the missing data were supplemented with data from the lake for which the strongest relationship was obtained.
Grouping lakes by probability of occurrence of ice cover employed hierarchical cluster analysis (CA). The cluster analysis was conducted by means of the Ward method, adopting square Euclidean distance as the probability measure. The identification of factors determining the course of ice cover in selected groups of lakes involved the application of the principal component analysis (PCA). The PCA was conducted separately for each group of lakes, and for the entire dataset combined. The PCA analysis made it possible to investigate general patterns in the relationship between the presence of ice cover in lakes and the lakes' morphometric parameters. The following morphometric lake parameters were considered: surface area (Area), altitude (Alt), maximum depth (MaxD), average depth (MeanD), Remote Sens. 2023, 15, 3030 6 of 18 volume (Vol), shoreline development (ShorD), exposure index (I), width (MaxW), length (MaxL), and effective length (EffL). The PCA analysis permits the presentation of data in new (orthogonal) spaces defined by principal components (PC). The number of PC was selected based on a Kaiser criterion of eigenvalues higher than 1. The relationship between PCs and analysed parameters was classified according to the values >0.75, 0.75-0.50, and 0.50-0.30 proposed by Liu et al. [33] as strong, moderate, and weak, respectively.
Finally, the degree of ice cover of lakes was analysed in reference to accumulated freezing degree days (AFDD) based on daily air temperature measurements for the Piła station ( Figure 1). Freezing degree days (FDD) were calculated for each day from November to April as follows: where: Ta-average daily air temperature in degrees Celsius.
A negative FDD value represented a temperature warmer than freezing, while a positive FDD value represented a temperature below freezing. The FDD values for each day of winter were added together to determine the AFDD for each winter season. Values of the coefficient of correlation between ice cover on the lake on a specific date and the AFDD value was calculated. Due to the strong relationship between lake water temperatures and air temperature at the Pila meteorological station during the AFDD calculation, no lapse rate adjustment was used to adjust the temperature from the meteorological station elevation to the lakes' elevations.

Temporary Changes in Lake Ice Cover
This section presents the research material for the analysis of ice cover on lakes of the Drawskie Lakeland, obtained for the period from November to April (winter halfyear). In that period, due to the periodical occurrence of air temperatures below zero, ice cover occurred on the lakes, beginning with partial ice cover leading to full ice cover. The number of obtained images in particular periods and by subsequent satellites of the Landsat programme is presented in Figure 2a.
The highest number of scenes (10) was obtained in the hydrological winter half-year 2020, and no cloudless scene was obtained for the hydrological winter half-year 2003. The highest number of scenes for analysis was provided by satellite Landsat 5-87 scenes-and satellites Landsat 4 and Landsat 9 provided 1 and 2 scenes, respectively ( Figure 2b). Only scenes from Landsat 5 documented full ice cover-defined as the occurrence of full ice cover on all the analysed lakes. Scenes from satellites Landsat 4 and Landsat 9 documented no occurrence of ice cover on any of the analysed lakes. The analysis of the collected data from particular months revealed that the highest number of scenes for analysis was available for April (70), and the lowest for December (8) (Figure 2c). In November, December, and April, full ice cover on all analysed lakes did not occur, and in January, February, and March, full ice cover on all lakes was documented in one, nine, and four images, respectively. Scenes from November, December, and April show no full ice cover on all analysed lakes. Scenes with partial ice cover on lakes are the most valuable from the point of view of the study. A total of 39 such scenes were obtained, 1 for November, 4 for December, 9 for January, 10 for February, 13 for March, and 2 for April. Results obtained from the analysis of these images were used for the assessment of the effect of lake morphometric parameters on the occurrence of ice cover. The occurrence of ice cover on particular lakes was determined as the ratio of ice cover surface area and water surface area on the lake. The range of ice cover on lakes in particular terms was highly variable (Figure 3). Six scenes showed a situation where parts of the lakes were fully covered with ice (value 1), and the remaining lakes remained devoid of ice (value 0). The lowest variability was observed in scenes from 9 March 2006, where the value was in a range from 0.97 to 1.00. In 22 scenes, at least one lake showed full ice cover; in 16 scenes, at least one lake was devoid of ice. The

Spatial Changes in Lake Ice Cover
Considering the variability in ice cover on particular lakes, it was revealed that it changed in a range from 0 (no ice cover) to 1 (full ice cover). The mean value of ice cover on lakes was at a level from 0.37 (Lake Drawsko) to 0.84 (Lake Wierzchnówko), and the median value was from 0.27 (Lake Lubie) to 0.98 (Lake Kiełpino) (Figure 4).

Spatial Changes in Lake Ice Cover
Considering the variability in ice cover on particular lakes, it was revealed that it changed in a range from 0 (no ice cover) to 1 (full ice cover). The mean value of ice cover on lakes was at a level from 0.37 (Lake Drawsko) to 0.84 (Lake Wierzchnówko), and the median value was from 0.27 (Lake Lubie) to 0.98 (Lake Kiełpino) (Figure 4).  The spatial distribution of mean ice cover on lakes and the median values are presented in Figure 5. West of Lake Drawsko, there are generally lakes with mean ice values ranging from 0.37 to 0.64, while east of it, there are mostly lakes with mean ice values ranging from 0.72 to 0.84 (Figure 5a). A similar tendency occurs with respect to median values. West of Lake Drawsko, there are lakes with median values of ice cover from 0.90 to 0.98, and east of it, there are lakes with median values of ice cover from 0.27 to 0.80 (Figure 5b).
Data on the occurrence of ice cover for particular lakes were used in the cluster analysis, aiming to group lakes depending on the degree of ice cover observed in subsequent scenes from Landsat satellites. The obtained results showed that the analysed lakes could be divided into two groups, A and B. Group A included 13 lakes, and group B 54 lakes ( Figure 6). In group B, three subgroups of lakes characterised by a different range of occurrence of ice cover were designated. Group B-1 included 14 lakes, group B-2-27 lakes, and group B-3-13 lakes. Data on the occurrence of ice cover for particular lakes were used in the cluster analysis, aiming to group lakes depending on the degree of ice cover observed in subsequent scenes from Landsat satellites. The obtained results showed that the analysed lakes could be divided into two groups, A and B. Group A included 13 lakes, and group B 54 lakes ( Figure 6). In group B, three subgroups of lakes characterised by a different range of occurrence of ice cover were designated. Group B-1 included 14 lakes, group B-2-27 lakes, and group B-3-13 lakes.  Data on the occurrence of ice cover for particular lakes were used in the cluster anal ysis, aiming to group lakes depending on the degree of ice cover observed in subsequen scenes from Landsat satellites. The obtained results showed that the analysed lakes could be divided into two groups, A and B. Group A included 13 lakes, and group B 54 lake ( Figure 6). In group B, three subgroups of lakes characterised by a different range of oc currence of ice cover were designated. Group B-1 included 14 lakes, group B-2-27 lakes and group B-3-13 lakes.  For lakes classified into group A, average ice cover varied from 0.37 to 0.61 (mean 0.52), and the median value was from 0.27 to 0.90 (mean 0.64). For lakes in group B, the ice cover values ranged from 0.56 to 0.84, and the median values from 0.73 to 0.98. Within particular subgroups, degrees of ice cover averaged from 0.63 to 0.81 for B-1, from 0.56 to 0.74 for B-2, and from 0.69 to 0.84 for B-3. The occurrence of ice cover within the designated groups and subgroups is presented in Figure 7. The statistical analysis of these results by means of a nonparametric Mann-Whitney U test showed that the average ice cover values in particular groups and subgroups differed at a significance level of 0.05. Considering the median values, no differences were found between groups B-1 and B-2, or between B-1 and B-3.

FOR PEER REVIEW 10 of 18
For lakes classified into group A, average ice cover varied from 0.37 to 0.61 (mean 0.52), and the median value was from 0.27 to 0.90 (mean 0.64). For lakes in group B, the ice cover values ranged from 0.56 to 0.84, and the median values from 0.73 to 0.98. Within particular subgroups, degrees of ice cover averaged from 0.63 to 0.81 for B-1, from 0.56 to 0.74 for B-2, and from 0.69 to 0.84 for B-3. The occurrence of ice cover within the designated groups and subgroups is presented in Figure 7. The statistical analysis of these results by means of a nonparametric Mann-Whitney U test showed that the average ice cover values in particular groups and subgroups differed at a significance level of 0.05. Considering the median values, no differences were found between groups B-1 and B-2, or between B-1 and B-3.

Factors Affecting the Lake Ice Cover Diversity
The analysis of lake morphometric parameters in particular groups designated by means of the CA method showed that group A includes lakes with the largest surface area, mean and maximum depth, volume, length, and width. The lakes also show the lowest altitude. The differences are statistically significant at the significance level of 0.05 (Table  1). Significant differences in morphometric parameters occur between subgroups B-1 and B-2. Subgroup B-2 includes lakes with almost all morphometric parameters showing values higher than those in lakes categorised to subgroup B-1, with the exception of altitude, maximum width, and exposure index. Lakes included in subgroups B-2 and B-3 differ only in terms of height of location, maximum and mean depth, and shoreline development index. No significant differences in morphometric parameters occur between lakes classified into subgroups B-1 and B-3. The spatial projection of results of the CA analysis revealed that lakes included in subgroups B-1 and B-3 are diverse in terms of location (Figure 8). Lakes included in subgroup B-1 are generally located in the west of the Drawskie Lakeland, and those in subgroup B-3 in the east.

Factors Affecting the Lake Ice Cover Diversity
The analysis of lake morphometric parameters in particular groups designated by means of the CA method showed that group A includes lakes with the largest surface area, mean and maximum depth, volume, length, and width. The lakes also show the lowest altitude. The differences are statistically significant at the significance level of 0.05 (Table 1). Significant differences in morphometric parameters occur between subgroups B-1 and B-2. Subgroup B-2 includes lakes with almost all morphometric parameters showing values higher than those in lakes categorised to subgroup B-1, with the exception of altitude, maximum width, and exposure index. Lakes included in subgroups B-2 and B-3 differ only in terms of height of location, maximum and mean depth, and shoreline development index. No significant differences in morphometric parameters occur between lakes classified into subgroups B-1 and B-3. The spatial projection of results of the CA analysis revealed that lakes included in subgroups B-1 and B-3 are diverse in terms of location (Figure 8). Lakes included in subgroup B-1 are generally located in the west of the Drawskie Lakeland, and those in subgroup B-3 in the east.
The principal component analysis showed that parameters that differentiate the degree of ice cover between the analysed lakes included surface area, lake volume, maximum and mean depth, shoreline length, and height of the lake's location above sea level (Figure 9). For lakes located at the lowest height, less ice cover was recorded, and lakes with high values of surface area, volume, maximum and mean depth, and shoreline length showed less ice cover. Table 1. Differences between morphometric parameters of lakes within groups designated by means of the CA method.

Parameters A vs. B-1 A vs. B-2 A vs. B-3 B-1 vs. B-2 B-1 vs. B-3 B-2 vs. B-3
Area (ha)  The principal component analysis showed that parameters that differentiate the degree of ice cover between the analysed lakes included surface area, lake volume, maximum and mean depth, shoreline length, and height of the lake's location above sea level (Figure 9). For lakes located at the lowest height, less ice cover was recorded, and lakes with high values of surface area, volume, maximum and mean depth, and shoreline length showed less ice cover.
The determination of the correlation between ice cover and climate conditions involved the analysis of the dependency of its surface area on accumulated freezing degree days (AFDD). AFDD constitutes the result of the accumulation of freezing temperatures for a given season (with consideration of snowmelt), and is a parameter commonly used in research on ice cover on lakes [34][35][36]. The obtained results show a significant correlation between ice cover surface area and AFDD for 54 lakes (Figure 10). The highest value was obtained for the largest lakes, namely, Drawsko, Lubie, Siecino, and Wielimie, where the correlation coefficient values were 0.71, 0.60, 0.56, and 0.54, respectively. The total surface area of ice on all the analysed lakes also correlated with AFDD at a level of 0.61. Some The determination of the correlation between ice cover and climate conditions involved the analysis of the dependency of its surface area on accumulated freezing degree days (AFDD). AFDD constitutes the result of the accumulation of freezing temperatures for a given season (with consideration of snowmelt), and is a parameter commonly used in research on ice cover on lakes [34][35][36]. The obtained results show a significant correlation between ice cover surface area and AFDD for 54 lakes (Figure 10). The highest value was obtained for the largest lakes, namely, Drawsko, Lubie, Siecino, and Wielimie, where the correlation coefficient values were 0.71, 0.60, 0.56, and 0.54, respectively. The total surface area of ice on all the analysed lakes also correlated with AFDD at a level of 0.61. Some limitations of the obtained results may be due to the distance of the Piła station outside of the Drawsko Lakeland area and the fact, that during the calculations, no lapse rate adjustment was used to adjust the temperature from the Piła station elevation to the lakes' elevations. On the other hand, air temperature showed high correlations with the temperature of the surface water layer (for monitored lakes: Drawsko, Lubie, Komorze), i.e., in the zone where ice cover develops.  The analysis of the spatial range of ice occurrence on particular lakes from subsequent terms permi ed the preparation of maps presenting the probability of occurrence of ice cover in reference to the largest lakes in the study area ( Figure 11).  The analysis of the spatial range of ice occurrence on particular lakes from subsequent terms permi ed the preparation of maps presenting the probability of occurrence of ice cover in reference to the largest lakes in the study area ( Figure 11). The analysis of the spatial range of ice occurrence on particular lakes from subsequent terms permitted the preparation of maps presenting the probability of occurrence of ice cover in reference to the largest lakes in the study area ( Figure 11). The results show a higher probability of occurrence of ice cover for lakes Komorze, Wielimie, and Wierzchowo than for lakes Drawsko, Lubie, and Siecino. Such a situation largely depends on the surface area and depth of individual lakes. In Lake Drawsko, in some areas of its central part, ice cover occurred only in 30-40% of the analysed Landsat satellite images. Ice cover was most frequently encountered in the narrowest northern and northwestern bays. On the other hand, Lake Wielimie, with a surface area approximate to that of Lake Drawsko, but with a smaller depth (mean depth = 2.2 m), showed a high probability of occurrence of ice cover.
Observations of ice cover in Polish lakes conducted to date have been based on in situ measurements, and cover a group of only several tens of lakes, narrowed down to objects with large surface areas. This context reveals a gap in reference to smaller lakes (the most abundant), important from the point of view of the functioning of the catchments in which they occur. So far, in the times of dynamic, or even explosive, development of teledetection, the tools have not been so far applied for the assessment of ice conditions in Polish lakes. It should also be emphasised that research on ice cover on lakes in many regions around the globe currently uses increasingly advanced remote sensing techniques [37] providing accurate and detailed information on the subject. The remote sensing of frozen rivers and lakes can, to a greater degree, contribute to the investigation of the role of water The results show a higher probability of occurrence of ice cover for lakes Komorze, Wielimie, and Wierzchowo than for lakes Drawsko, Lubie, and Siecino. Such a situation largely depends on the surface area and depth of individual lakes. In Lake Drawsko, in some areas of its central part, ice cover occurred only in 30-40% of the analysed Landsat satellite images. Ice cover was most frequently encountered in the narrowest northern and northwestern bays. On the other hand, Lake Wielimie, with a surface area approximate to that of Lake Drawsko, but with a smaller depth (mean depth = 2.2 m), showed a high probability of occurrence of ice cover.
Observations of ice cover in Polish lakes conducted to date have been based on in situ measurements, and cover a group of only several tens of lakes, narrowed down to objects with large surface areas. This context reveals a gap in reference to smaller lakes (the most abundant), important from the point of view of the functioning of the catchments in which they occur. So far, in the times of dynamic, or even explosive, development of teledetection, the tools have not been so far applied for the assessment of ice conditions in Polish lakes.
It should also be emphasised that research on ice cover on lakes in many regions around the globe currently uses increasingly advanced remote sensing techniques [37] providing accurate and detailed information on the subject. The remote sensing of frozen rivers and lakes can, to a greater degree, contribute to the investigation of the role of water bodies in hydrology and effects of environmental variability and change in the ice and hydrological system [38].
In this context, research was undertaken aiming to determine the usefulness of publicly available Landsat images for the analysis of ice cover on Polish lakes. This research subject corresponds with earlier studies employing such images [39][40][41]. Landsat images permit the analysis of ice cover on small lakes, although according to Zhang and Pavelsky [42], certain limitations also exist resulting from the revisit time (16 days for Landsat 4 onwards), as well as the occurrence of cloudiness, additionally limiting the availability of data from the images. On the other hand, there is a simple possibility of ice detection based on Landsat satellites with simultaneous application of the rich archive of the images [43]. This is particularly important for areas such as Poland, where this type of analysis has never been conducted before. The results showed that the designated study area is characterised by high variability in ice cover on its lakes. Therefore, the relatively small surface area of the Drawskie Lakeland (with generally uniform climate conditions) clearly points to the individual parameters of lakes as key for the occurrence of ice cover. This statement is reflected in the statistical analyses and significance of particular morphometric parameters. Lakes with the largest surface area, the middle parts of which were covered with ice the most, seldom are the most prone to wind action. According to Ke et al. [44] the delay in terms of the freezing of Lake Khanka is strongly correlated with an increase in wind speed. In the case of Lake Ulansuhai (China), wind was the main driver of development of ice cover [45]. The role of wind for the course of ice cover in five lakes in Poland was assessed by Bartosiewicz et al. [46]. According to the authors, the effect of wind delays ice development in autumn. It should be emphasised that particular parameters combined provide conditions less or more favourable for the occurrence of ice cover in lakes. Research by Williams et al. [47] evidenced that surface area and mean depth are positively correlated with the appearance of ice cover. The spatial pattern of freezing and defrosting of Lake Qinghai may, to a certain degree, reflect the difference in water depth in the lake within the designated areas; in the west, similar spatial patterns of freezing and ablation occurred, whereas in the eastern part, they were reversed [36]. The analysis of more than a hundred lakes in Norway showed a considerable effect of lake surface area because small lakes froze earlier due to their smaller volume [48]. The analysis of the frequency and range of ice cover in Lower Lake Constance with three basins with variable morphology revealed that full ice cover occurs almost each winter in the shallowest part, most isolated, with no substantial inflows [49]. Situations described in the said study, pointing to the variability in ice cover within particular lakes, owing to remote sensing techniques, were, for the first time, conclusively confirmed for the selected Polish lakes. In bays, the shallowest zones of the lakes protected from wind, ice cover was recorded the most frequently, as evident in the case of lakes Drawsko, Lubie, Wierzchowo, and Siecino. In the case of the latter, Ostrów Island plays a considerable role (in the southeastern part). It isolates the part of the lake with recorded high frequency of occurrence of ice cover. The image of spatial occurrence of ice cover in Lake Wielimie is interesting. Despite its large surface area, but simultaneously low mean depth, it maintains a relatively uniform distribution in terms of this characteristic. Two zones deserve particular attention. The first one is the northwestern part of the lake-the mouth zone of the Gwda River-where the probability of occurrence of ice cover is at a level of 50-60% (Figure 11c). The determination of the interaction between river flow and ice cover corresponds with earlier studies of the type, where, among others in the case of the Great Slave Lake, the inflow of the Slave River plays a substantial role in the occurrence of ice cover [50]. The second characteristic zone in Lake Wielimie is its southwestern part with the inflow of the urban Nizica River. Flowing through an urbanised area (Szczecinek), the river supplies waters with different physicochemical properties than those in the lake. In the winter period, the supply of salt used for snowploughing roads and pavements is important. With surface flow, it is supplied to the urban river and then reaches Lake Wielimie. The role of anthropogenic pollutants for the development of ice cover in water bodies was emphasised by Stolarski and Rzetała [51]. According to the authors, in such conditions, the ice layer does not form at all (even in the conditions of a substantial decrease in air temperature), or it is very thin.
The importance of ice cover on lakes is particularly great in the context of global warming, where numerous studies point to changes in the ice regime [52][53][54], and the situation is considered an evident indicator of climate change [55,56]. The changes will progress over the next several decades [46,57], considerably changing the current lake parameters shaped by the presence of ice cover. It is therefore important to use all the available techniques and possibilities of monitoring the distribution and course of ice cover on lakes. The issue is of importance not only in reference to cases subject to continuous observations, but also for ones to which such an approach is not applied. Appropriate decision making in the scope of water management and activities aiming to optimise of the functioning of lake ecosystems in the scope of both biotic and abiotic conditions requires detailed knowledge, and, as evidenced by the obtained results, remote sensing techniques can constitute an important source of information for Polish lakes.

Conclusions
Observations of ice cover in Polish lakes conducted to date provide several-decadeslong datasets, and are based on field measurements. They generally refer to lakes with relatively large surface areas covering a total of several tens of lakes, which is a small group in the context of the entire population. Remote sensing techniques have not been applied so far, whereas their importance and accuracy in data collection is dynamically increasing on a year-to-year basis in different regions of the globe. Satellite images from the Landsat mission permitted the description of ice cover in lakes never monitored before in those terms (mainly small ones), and the assessment of the spatial variability in the occurrence of ice cover within individual lakes. Despite the relatively small study area, considerable variability in the occurrence of ice cover was found, and the morphometric parameters of lakes proved key in this aspect. The implementation of the methodical assumptions referring to satellite images and the obtained results offers new possibilities in the context of expanding knowledge on the occurrence of ice in Polish lakes. In the times of climate change and the progressing transformation of the natural environment, this approach is complementary for activities in the scope of the management of lake ecosystems, in reference to both biotic and abiotic aspects. In the case of multitemporal analyses, attention should be paid to changes in water surface area caused by changes in water levels and overgrowth. Changes in water level in moraine lakes and coastal lakes with smooth shores can significantly affect the surface of the water table and the results of ice cover analysis based on satellite data. In tunnel lakes with steep shores, changes in water level affect the surface of the water table and, therefore, the analysis results, to a lower extent.