Ice Phenology in Eurasian Lakes over Spatial Location and Altitude

: Eurasian freezing lakes cover an almost 180 ◦ wide longitude sector between the latitudes 30 ◦ and 75 ◦ N, and their altitudes range from below the sea surface level up to 5 km elevation. Ice phenology varies widely in this region. However, these variations and their inﬂuence factors have been little studied. Analytic models are applied here to examine these variations supported by historical ice and weather data. These models are forced by a linear air–lake heat exchange formula based on local empirical ﬁts. The weather brings latitude–longitude–altitude patterns to the large-scale lake ice phenology. Freezing and breakup dates are forced by the local air temperature and solar radiation, and their rates of change are also important. In addition, freezing depends on lake depth and breakup depends on accumulated ice thickness. Lake depth provides a lag and radiation balance provides a shift with respect to the air temperature in cooling of the lake, and breakup is dictated by spring warming conditions and ice thickness. Due to solar radiation forcing, the common degree-day approach is biased for modelling ice phenology, especially in low latitudes. Analytic models provide a ﬁrst-order tool for climate sensitivity of ice seasons. The freezing date and breakup date both change by around ﬁve days per one-degree shift in air temperature away from the climatological ice margin; however, at this margin, the sensitivity is higher.


Introduction
In the Eurasian continent, lakes freeze between the latitudes of 30-75 • N. The extent is at most south over the high Tibetan Plateau, while freezing reaches only down to 50 • N in western Europe due to the heating provided by the North Atlantic Current (Figure 1). The northern boundary of the continent is the coast of the Arctic Ocean at 70-75 • N. Freezing and ice breakup bring major changes to the physics and ecology of lakes and to the human living conditions in lake districts [1]. The primary factors that drive the ice season are the air-lake interaction, radiation balance, and the characteristics of individual lakes [2]. Climate variations have a major impact on ice season characteristics; in other words, ice season characteristics are sensitive indicators of climate.
An ice phenology time-series consist of the dates of freezing and ice breakup and the number of ice days. There is a large set of these time-series, and many of them go back 100-200 years in history [3]. The thermal memory of lakes is too short in open water periods to show significant dependence between consecutive ice seasons. Therefore, the freezing date is independent of the previous ice season; however, in winter, time slows down, and there is a weak statistical connection between the freezing date and the following breakup date [4]. A connection through the ice season can also be seen between deep-water temperature and dissolved oxygen content [1]. During windy autumn seasons, forced deep mixing may continue until the water temperature is well below the temperature of Figure 1. The Eurasian continent with contours for lake ice season of 100 and 180 days [6], together with the mean January 0 °C contour of air temperature. The study sites are marked by red circles.
Ice phenology time series predominantly show decrease in the length of the ice season in Eurasian lakes over the last 100 years [1,3,[7][8][9][10]. The shortening is seen nearly equally in the dates of freezing and breakup. The trend comes from general climate warming in the 20th century in Eurasia and North America, where most of the lake ice investigations were made. Aperiodic variations and noise are superposed on the trend toward a milder ice climate. In the Tibetan Plateau, shortening of ice seasons is seen in satellite remote sensing data; however, some saline lakes show opposite behavior, because their salinity has been decreasing due to the increasing water volume that has influenced the cooling process [8,[11][12][13].
In addition to the phenology, the thicknesses of ice and snow have been recorded in many lakes, and their annual maxima are a measure of the severity of ice seasons. In large lakes, the evolution of ice coverage has been monitored with its annual maximum as another characteristic of season severity [9]. A statistical approach is normally taken in analyses of lake ice time series [4,7,14,15].
The connection of the statistical characteristics of ice phenology to the physics of lakes needs to be understood, particularly to evaluate the impact of predicted future climate scenarios on the ice season in lakes. This can be approached by analytical modelling [9] that guides to understand the role of climatological forcing factors and lake properties, opens more physically-based views into the properties of ice time series, and reveals the climate sensitivity and future of winter conditions in Eurasian lakes.
Lake ice phenology is strongly related to weather conditions, as shown by field experiment outcomes [16,17] and numerical modelling [18,19]. Air temperature depends both on altitude and spatial location, while solar radiation depends primarily on latitude. The annual air temperature cycle can be similar in low-latitude high mountains and lowaltitude polar regions, but they receive quite different solar radiation forcing. There is also a longitudinal temperature variation over Eurasia due to the variations in westerlies from the North Atlantic (the North Atlantic Oscillation, or NAO) and the continental climate impact that causes major differences in the heat balance, ice conditions, and primary production in Eurasian freezing lakes [2,20]. Latitude and altitude are expected to have different roles in ice phenology via radiation balances. This paper examines the spatial characteristics of lake ice phenology in Eurasia. The study reference sites cover the Tibetan plateau and boreal and tundra zone in northern Figure 1. The Eurasian continent with contours for lake ice season of 100 and 180 days [6], together with the mean January 0 • C contour of air temperature. The study sites are marked by red circles.
Ice phenology time series predominantly show decrease in the length of the ice season in Eurasian lakes over the last 100 years [1,3,[7][8][9][10]. The shortening is seen nearly equally in the dates of freezing and breakup. The trend comes from general climate warming in the 20th century in Eurasia and North America, where most of the lake ice investigations were made. Aperiodic variations and noise are superposed on the trend toward a milder ice climate. In the Tibetan Plateau, shortening of ice seasons is seen in satellite remote sensing data; however, some saline lakes show opposite behavior, because their salinity has been decreasing due to the increasing water volume that has influenced the cooling process [8,[11][12][13].
In addition to the phenology, the thicknesses of ice and snow have been recorded in many lakes, and their annual maxima are a measure of the severity of ice seasons. In large lakes, the evolution of ice coverage has been monitored with its annual maximum as another characteristic of season severity [9]. A statistical approach is normally taken in analyses of lake ice time series [4,7,14,15].
The connection of the statistical characteristics of ice phenology to the physics of lakes needs to be understood, particularly to evaluate the impact of predicted future climate scenarios on the ice season in lakes. This can be approached by analytical modelling [9] that guides to understand the role of climatological forcing factors and lake properties, opens more physically-based views into the properties of ice time series, and reveals the climate sensitivity and future of winter conditions in Eurasian lakes.
Lake ice phenology is strongly related to weather conditions, as shown by field experiment outcomes [16,17] and numerical modelling [18,19]. Air temperature depends both on altitude and spatial location, while solar radiation depends primarily on latitude. The annual air temperature cycle can be similar in low-latitude high mountains and lowaltitude polar regions, but they receive quite different solar radiation forcing. There is also a longitudinal temperature variation over Eurasia due to the variations in westerlies from the North Atlantic (the North Atlantic Oscillation, or NAO) and the continental climate impact that causes major differences in the heat balance, ice conditions, and primary production in Eurasian freezing lakes [2,20]. Latitude and altitude are expected to have different roles in ice phenology via radiation balances. This paper examines the spatial characteristics of lake ice phenology in Eurasia. The study reference sites cover the Tibetan plateau and boreal and tundra zone in northern Europe. Analytic models are employed for the dates of freezing and breakup based on spatial location, altitude, and lake depth, and forced by solar radiation and air temperature. The models are calibrated using the climatology of the study sites. Parameterized analytic models are used to examine the ice climatology of Eurasian lakes and the climate sensitivity of the lake ice season over the continent.

Lakes
A set of lakes was taken for the calibration of the climatological lake ice model (Table 1). They include lakes from tundra and boreal zones, and the Central Asian arid, cold climate zone. The latitude range is 34 • , the longitude range is 94 • , and the altitude range is 4.3 km in the sites. The data present the atmospheric and ice climatology obtained from the home laboratories of the authors, open data bases, and publications. The date of freezing normally refers to the date when the lake froze or the date when the observation site was frozen. In the past, phenological observations were made from the shore and the observation area was limited, while current satellite data provide surface information for the whole lake. The breakup date is the last day when ice was observed. Sometimes, the beginning of the ice decay period is included in the analysis. It is defined as the day the thickness of ice starts to decrease after the annual maximum, or as the day the continuous ice cover of a lake experiences final breakage. The latter approach is useful when air-lake gas exchange questions are examined. However, field data have shown that in a fixed lake, the length of the decay period does not vary much but its starting time does [2,9,17]. Thus, the final breakup day and the start of the decay are highly correlated.
If a lake does not freeze every year, the statistical treatment of freezing and breakup dates needs care since standard statistics such as the mean and variance cannot be properly defined. The binary variable (1 for ice and 0 for no ice) is then taken into the ice phenology data to account for the statistics of ice occurrence, and fractiles of the phenological distributions are used for the statistics [6,21], e.g., median (50 % fractile) date t 0.5 can be used instead of average when the probability of freezing is at least 0.5; the median freezing date t F0.5 is obtained from P(t F ≤ t F0.5 ) = 0.5, where t F is the annual freezing date. One lake included in this study was from the climatological ice margin where the probability of annual ice occurrence is between zero and one. In large lakes, although near-shore areas freeze annually, a complete freeze-over does not always take place. The maximum annual coverage is then recorded as another ice season characteristic [9].

Data
The wind, air temperature, and humidity were obtained from the Finnish Meteorological Institute and the China Meteorological Forcing Dataset for a detailed preparatory study of the heat budget. The freezing date and air temperature of lakes Kilpisjärvi, Kallavesi, and Pääjärvi were provided by the Finnish Meteorological Institute and Finnish Environmental Institute. In the climatology anaysis, air temperature  in Eurasia was otained from NCEP/NCAR Reanalysis 1, 2.5 • × 2.5 • latitude-longitude global grid.

Modelling Approach
Freezing of a lake, ice growth, and breakup are forced by the heat exchange with atmosphere-ice-water body-lake bottom interactions, and solar radiation. These interaction processes act on the boundary surfaces, while solar forcing provides heating throughout the ice and into the water body, and a part of the solar heat absorbed by the water is conducted to the ice bottom [15,22]. The basic analytical models are the slab model for lake cooling [23], freezing-degree-day model for ice growth [24], and positive-degree-day models for ice melting [25]. The connection of ice melting to air temperature is indirect, and models based on energy balance are more physically based [26].
Freezing date: Analytic modelling of autumn cooling is based on the slab model. The water temperature is assumed to be vertically uniform, and we have: where T w is the mixed-layer temperature, H is the mixed-layer depth or lake depth (fixed), t is time, ρ w and c w are the density and specific heat of the water, respectively, and Q 0 and Q b are the heat fluxes at the lake surface and bottom, respectively. The freezing date t F is taken as the date when the mixed-layer temperature has reached the freezing point, T w (t F ) = T f . Ice melting: Modelling the decay of lake ice is in principle straightforward, since in the melting season, ice loss is due to the net gain of energy from the external fluxes, and heat conduction is not significant. However, the space-time variability of the albedo and light attenuation coefficient of melting ice gives rise to parametrization problems for the solar flux. Ice melting progresses as: where h is ice thickness, Q w is the heat flux from water to ice, ρ i is ice density, and L f is the latent heat of freezing. The breakup date t B is taken from h(t B ) = 0 after the thickness has reached its annual maximum h = h max . To predict the breakup date, this thickness needs to be known. It can be estimated as h max = where a ≈ 2-3 cm ( • C d) −1 and b ≈ 10 cm are the model parameters; S is the cumulative temperature of freezing-degreedays [24,27]. Melting of ice also depends on the ice structure, in particular on the fractions of clear congelation ice and opaque snow-ice. Light transfer properties largely determine how much melting takes place on the top and bottom surfaces and in the interior. However, the net heat flux can be taken as the first order forcing factor.
Forcing: These models are forced by solar radiation, atmospheric heat flux, and heat flux from below. The surface+internal forcing is expressed in linear form, which is convenient in analytic analysis (see Appendix A): where K 0 and K 1 are the time-dependent model parameters, and T a and T 0 are the air temperature and surface temperature, respectively. The parameter K 0 is governed by solar radiation, while K 1 builds up from turbulent heat exchange. These parameters were estimated for the control sites in Finland and Tibet from climate data ( Figure 2). The forcing from below the mixed layer (or lake bottom) can be formulated as where K b is the heat transfer coefficient and T b is the temperature below the mixed layer (or lake bottom). and amplitude lower in the Tibetan site, compared to the sites in Jokioinen, southern Finland. The annual averages are 120.5 W m -2 in Ngoring (range 150 W m -2 in monthly means) and 17.8 W m -2 (range 170 W m -2 in monthly means) in Jokioinen, which represents Lake Pääjärvi. The parameter is rather stable and the level is lower in Tibet than in Finland, and the annual cycle is mostly due to the temperature dependence of the saturation water vapour pressure. The annual averages are 11.7 W m -2 °C -1 in Ngoring and 18.6 W m -2 °C -1 in Jokioinen. for Lake Pääjärvi (based on climate data from Jokioinen, 60°49′ N 23°30′° E) and Ngoring Lake (34°54′ N 97°42′ E).

Geography of Ice Phenology
The study lakes are located within latitudes 34-69° N, longitudes 14-109° E, and altitudes 0.1-4.3 km (Table 1). In all cases, the mean air temperature is below the freezing point in December-March. The northern lakes freeze in November and the Tibetan lakes freeze in December. The air temperature is quite low in Tibet in December (Table 2), which shows that solar radiation places a long delay on the freezing, compared to the north. The very shallow Lake Wuliansuhai in Inner Mongolia freezes early, and the very deep Lake Baikal freezes only in January. Close to the climatological margin of freezing lakes, the freezing date and breakup date have a large variability [15,28].
Ice formation, growth, and melting depend on the latitude , longitude Λ, and altitude via air temperature and solar radiation. Air temperature decreases upward by the adiabatic lapse rate Γ ≈ 6 °C km , while, according to climatology, latitude-wise and longitude-wise cooling in Eurasian winter are about 1.2 °C deg lat and 0.3 °C deg lon , respectively. Thus, 1 km increase in altitude corresponds to 5 degrees increase in latitude and 20 degrees in longitude. The location difference between Ngoring and Kilpisjärvi is 34 degrees in latitude, 77° in longitude, and 3.8 km in altitude, and therefore their temperatures should be close. In autumn and winter, Ngoring is colder but in spring and summer, the graphs are close (Figure 3). A similar condition exists between Qinghai and Kallavesi.
Comparing Kilpisjärvi with the Tibetan lakes, it was seen that the December air temperature is lower than in Qinghai but higher than in the other sites. The freezing dates then show that the latitudinal, i.e., solar radiation, effect accounts for about one month between the Tibetan plateau and the European Arctic tundra. For similar air temperature evolution, freezing takes place about one month later in the Tibetan sites, as compared to  The parameter K 0 has a strong annual cycle from solar radiation. The level is higher and amplitude lower in the Tibetan site, compared to the sites in Jokioinen, southern Finland. The annual averages are 120.5 W m −2 in Ngoring (range 150 W m −2 in monthly means) and 17.8 W m −2 (range 170 W m −2 in monthly means) in Jokioinen, which represents Lake Pääjärvi. The parameter K 1 is rather stable and the level is lower in Tibet than in Finland, and the annual cycle is mostly due to the temperature dependence of the saturation water vapour pressure. The annual averages are 11.7 W m −2 • C −1 in Ngoring and 18.6 W m −2 • C −1 in Jokioinen.

Geography of Ice Phenology
The study lakes are located within latitudes 34-69 • N, longitudes 14-109 • E, and altitudes 0.1-4.3 km (Table 1). In all cases, the mean air temperature is below the freezing point in December-March. The northern lakes freeze in November and the Tibetan lakes freeze in December. The air temperature is quite low in Tibet in December (Table 2), which shows that solar radiation places a long delay on the freezing, compared to the north. The very shallow Lake Wuliansuhai in Inner Mongolia freezes early, and the very deep Lake Baikal freezes only in January. Close to the climatological margin of freezing lakes, the freezing date and breakup date have a large variability [15,28].
Ice formation, growth, and melting depend on the latitude Φ, longitude Λ, and altitude Z via air temperature and solar radiation. Air temperature decreases upward by the adiabatic lapse rate Γ ≈ 6 • C km −1 , while, according to climatology, latitude-wise and longitude-wise cooling in Eurasian winter are about 1.2 • C deg lat −1 and 0.3 • C deg lon −1 , respectively. Thus, 1 km increase in altitude corresponds to 5 degrees increase in latitude and 20 degrees in longitude. The location difference between Ngoring and Kilpisjärvi is 34 degrees in latitude, 77 • in longitude, and 3.8 km in altitude, and therefore their temperatures should be close. In autumn and winter, Ngoring is colder but in spring and summer, the graphs are close (Figure 3). A similar condition exists between Qinghai and Kallavesi.
Comparing Kilpisjärvi with the Tibetan lakes, it was seen that the December air temperature is lower than in Qinghai but higher than in the other sites. The freezing dates then show that the latitudinal, i.e., solar radiation, effect accounts for about one month between the Tibetan plateau and the European Arctic tundra. For similar air temperature evolution, freezing takes place about one month later in the Tibetan sites, as compared to the Arctic (Figure 3). These Tibetan lakes freeze when the air temperature is 10-15 • C below the freezing point. A corresponding difference appears in the breakup date, and the breakup takes place in the Tibetan plateau at air temperatures below the freezing point. Figure 3 clearly illustrates the limitations of the normal freezing-degree-days and positive-degree-days formulae for the freezing and breakup dates, suggesting that they need to be modified in respect to the radiation balance. The Szczecin lagoon is shown as the reference, where the annual probability of ice occurrence is lower than one (0.78).  [21,30] 15 Jan -0.5 °C +6 °C 1 Mar Wuliansuhai [31] 7 Dec -10 °C +4 °C 31 Mar Baikal [32] 10 Jan -18 °C +1 °C 1 May Gyaring [33] 30 Nov ± 21 d -18 °C -2 °C 21 April ± 17 d Ngoring [33] 15 Dec ±10 d -17 °C -2 °C 1 May ± 8 d Qinghai [33] 29 Dec ± 4 d -10 °C +4 °C 1 April ± 7 d Next, analytic solutions were derived for the first-order explanations of the characteristics of the field data on ice phenology. Forcing of the system is taken care of by the linearized atmospheric heat flux (Equation (3)). In November-December, solar radiation is nearly zero in Kilpisjärvi, but in Ngoring the daily mean is around 150 W m −2 according to observations; thus, Ngoring is able to resist freezing until the air temperature is down to about -10 • C. The daily solar flux is above 150 W m −2 all winter and its penetration into ice is strong due to the absence of snow in dry climates, e.g., [29].   [21,30] 15 Jan. −0.5 • C +6 • C 1 Mar. Wuliansuhai [31] 7 Dec. −10 • C +4 • C 31 Mar. Baikal [32] 10 Jan. −18 • C +1 • C 1 May. Gyaring [33] 30 Nov. ± 21 d −18 • C −2 • C 21 April ± 17 d Ngoring [33] 15 Dec. ± 10 d −17 • C −2 • C 1 May ± 8 d Qinghai [33] 29 Dec. ± 4 d −10 • C +4 • C 1 April ± 7 d Next, analytic solutions were derived for the first-order explanations of the characteristics of the field data on ice phenology. Forcing of the system is taken care of by the linearized atmospheric heat flux (Equation (3)).

Freezing Date
The slab model (1) with the surface heat flux (3) is first written as: Here, λ a and λ b are the relaxation times of the mixed layer for the surface and bottom forcing, and T R (Figure 4) is a temperature shift term originating from the radiation balance. This equation can be directly integrated: where λ = λ a + λ b , and τ L = λ −1 ∝ H is the response time of the lake with the proportionality coefficient of about 2 d m −1 .

Freezing Date
The slab model (1) with the surface heat flux (3) is first written as: Here, and are the relaxation times of the mixed layer for the surface and bottom forcing, and ( Figure 4) is a temperature shift term originating from the radiation balance. This equation can be directly integrated:  The fundamental question in ice phenology is whether a lake is ice-free or freezes over. Take a parabolic form of the winter temperature evaluation, expressed as = + , where is the minimum air temperature in winter and > 0 is a shape factor. Freezing takes place in late fall when the radiation balance is rather stable, and therefore the term is taken as a constant ( Figure 4). The lake bottom temperature is fixed. The condition of lake freezing is: where we introduced the temperature = − − − as the critical temperature, which provides a necessary condition for any lake to freeze as < . Then from = the shape factor is obtained as = 4 − , where is the length of the 'cold period' < . The condition (6) becomes simply > 2 , i.e., for a lake to freeze, the cold period must be at least twice the response time of the lake. The temperature condition is external, but the required length of the cold period is lake-specific since ∝ . The heat flux from the sediments is generally much smaller than the surface heat flux in the cooling season, but geothermal lakes may stay ice-free by condition (6). The fundamental question in ice phenology is whether a lake is ice-free or freezes over. Take a parabolic form of the winter temperature evaluation, expressed as T a (t) = ct 2 +T a , whereT a is the minimum air temperature in winter and c > 0 is a shape factor. Freezing takes place in late fall when the radiation balance is rather stable, and therefore the term T R is taken as a constant (Figure 4). The lake bottom temperature is fixed. The condition of lake freezing is: where we introduced the temperature T c = T f − T R − λ b λ a T b − T f as the critical temperature, which provides a necessary condition for any lake to freeze asT a < T c . Then from T a (t) = T c the shape factor c is obtained as cτ 2 f = 4 T c −T a , where τ f is the length of the 'cold period' T a < T c . The condition (6) becomes simply τ f > 2τ L , i.e., for a lake to freeze, the cold period must be at least twice the response time of the lake. The temperature condition is external, but the required length of the cold period is lake-specific since τ L ∝ H. The heat flux from the sediments is generally much smaller than the surface heat flux in the cooling season, but geothermal lakes may stay ice-free by condition (6).
Assume that the heat flux from the bottom is negligible. Equation (6) shows that a necessary condition for freezing of any lake isT a − T f < −T R ; e.g., T R ≈ −3 • C in Lake Pääjärvi and T R ≈ 12 • C in Ngoring Lake ( Figure 4). Thus, in the north, radiational and evaporation losses can freeze the surface even when the air temperature is above the freezing point, but in low latitudes, solar radiation can compensate for the losses down to quite low temperatures. For a lake of finite depth, the minimum air temperature must be still lower so that the length of the cold period is longer than twice the lake response time. e.g., if T R = 0 and λ b = 0, for a 5 m deep lake, the air temperature must be below the freezing point at least 10 days. In Szczecin Lagoon, Poland (53.8 • N, 14.1 • E, mean depth 4 m), the probability of ice occurrence is 0.78, and therefore in an average winter, the basin is just about to freeze. We haveT a ≈ 0 • C, T R ∼ −1 • C and τ L ≈ 8 d. This means that a period of 16 days colder than +1 • C is needed to create an ice cover that corresponds to the observed ice climatology.
Assuming freezing takes place, to examine the freezing date, a linear cooling law T a (t) = T f − αt is taken, where α is the rate of cooling and t = 0 refers to the time the decreasing air temperature passes the freezing point, T a (0) = T f . The water temperature is obtained from Equation (4), and thereafter the freezing time is obtained from T t f = T f : Thus, the solution contains delay terms due to the lake response time, heat flux from bottom, and radiation balance. The local climate impact is included in α, the time zero, and T R , while the response time and bottom heating are lake specific.
Assuming that the bottom heat flux is small, the freezing date is delayed behind t = 0 by the lake response time plus the lag ∆t = α −1 T R , t F = τ L + ∆t, from the air temperature. It is interesting to compare the cases of Lake Kilpisjärvi in Arctic tundra and Ngoring Lake in Tibetan tundra, since their temperature cycles and depths are close (Figure 3). In Kilpisjärvi, we have ∆t =−15 days; with τ L ≈ 2H d m −1 , the estimated response time is τ L = 39 days, and thus the predicted freezing time is t F =24 days, 6 days less than in the true climatology. In Ngoring Lake, ∆t = 50 days, the estimated response time is 35 days, and the predicted freezing date is 85 days, 5 days more than in the true climatology. Thus, the analytic model provides an acceptable estimator of the freezing date, and it is feasible in analysing lake ice climatology across latitude and altitude. For the other lakes, corresponding results are obtained. e.g., in Szczecin lagoon, ∆t = −5 d and τ L = 8 d, i.e., the basin freezes soon after the air temperature has gone to the cold level. The lake response time is a lake characteristic, while the lag due to forcing depends on the local climate.

Breakup Date
Analytic models of ice melting are based on the heat gain by the ice that is used to decrease the ice volume, and 'ice thickness' represents the volume of ice per unit area. This approach corresponds to the physical picture where significant melting takes place at both boundaries and in the ice interior.
Using the linear surface heat flux (Equation (3)), the thickness of ice during the melting period is: where h max = h(t M ) is the ice thickness in the beginning of melting. In the melting period, net solar radiation increases fast and therefore we take K 0 = K 0 (t) and thus T R = T R (t) (see Figures 2 and 4). At the start-up time t M , the surface heat flux turns positive after the maximum ice thickness has been reached. This time is approximated from the solution This is a delicate formula, since in spring both air temperature and solar radiation increase, and the temperature T R depends largely on the latitude and altitude (Figure 4). The date of breakup, t B is obtained from the implicit equation h(t B ) = 0. Since melting Water 2022, 14, 1037 9 of 13 adds on the incoming heat fluxes, the heat flux from water can be absorbed to the parameter K 0 , as is clear from Equation (8).
The K 1 term in Equation (8) provides a positive-degree-day formula. For K 1 = 15 W m −2 • C −1 , we have the degree-day coefficient K 1 ρ i L f = 0.42 cm ( • C d) −1 . When the melting period has begun, the air temperature and the K 0 -term are assumed to increase linearly by the rates β and γ, respectively. Ignoring Q w (or absorbing to K 0 ), we have an equation for the breakup date: The square of the period t B − t M illustrates the rather low sensitivity of the breakup date to ice thickness and the model parameters.
The parameters of Equation (10) for our study sites were obtained from the data shown in Figures 2 and 3. For Ngoring lake, γ ≈ 1 W m −2 d −1 , for Kallavesi γ ≈ 1.5 W m −2 d −1 , and for both K 1 ≈ 15 W m −2 • C −1 and β ≈ 0.17 • C d −1 . An approximation for the date of ice breakup is: where A ≈ 15-20 cm −1 Thus, in the first-order approximation, the length of the melting period depends on the initial ice thickness, and for half-meter ice, it is about 30 days. Figure 3 does show that the melt rates overlap, but the breakup depends on the start-up time and initial ice thickness. The maximum ice thickness is typically 70 cm in Ngoring Lake and 90 cm in Lake Kilpisjärvi. Furthermore, the estimated time t M is 15 March for Ngoring Lake, 31 March for Lake Kallavesi, and 10 May for Lake Kilpisjärvi, and the observed mean breakup dates are for these sites, respectively, 30 April, 8 May, and 18 June. In the observations, melting took 35-45 days from the start in these sites, while the start-up dates ranged over 65 days (Table 2, Figure 3). In the Tibetan lakes, at the time of ice breakup, the air temperature had just reached the freezing point. At the time of start-up, in Kallavesi, the temperature is at the freezing point, but in Kilpisjärvi, it is about 2 • C higher than the freezing point.

Discussion and Conclusions
The ice phenology time series from 10 Eurasian sites illustrated the roles of local air temperature (determined by altitude, latitude, longitude) and solar radiation (determined by latitude) to determine the statistics of ice occurrence, freezing date, and breakup date.
The role of solar radiation shows up strongly in Eurasia, where the latitudinal extent of freezing lakes goes from 35 • N to 70 • N. In the Eurasian continent, winters also become more severe toward the east due to continental vs. maritime climate ( Figure 5).
In autumn, surface temperature of a lake follows the air temperature with a lag depending in the lake depth and with a shift depending on the net solar radiation. Ice phenology depends on the lake morphology and weather. For freezing, the air temperature must go down enough and stay there long enough for freezing to occur. The critical temperature is close to the freezing point at 60 • N but less by 10 • C in the Tibetan plateau at 35 • N. After freeze-up, ice thickness follows the freezing-degree-days. The melting period begins when the surface heat balance turns positive in late winter. This date depends on the solar radiation and air temperature. Again, at 60 • N, the air temperature is then close to the freezing point, but in Tibetan plateau it is still around 10 degrees lower. In low latitudes, solar radiation has a dominant role in the ice season heat balance, and therefore statistical degree-day models do not work there. Rather, the level of solar radiation determines how low the air temperature must go before ice formation is possible. Under similar temperature conditions, ice season is longer in a high latitude than in high altitude. In autumn, surface temperature of a lake follows the air temperature with a lag depending in the lake depth and with a shift depending on the net solar radiation. Ice phenology depends on the lake morphology and weather. For freezing, the air temperature must go down enough and stay there long enough for freezing to occur. The critical temperature is close to the freezing point at 60° N but less by 10 °C in the Tibetan plateau at 35° N. After freeze-up, ice thickness follows the freezing-degree-days. The melting period begins when the surface heat balance turns positive in late winter. This date depends on the solar radiation and air temperature. Again, at 60° N, the air temperature is then close to the freezing point, but in Tibetan plateau it is still around 10 degrees lower. In low latitudes, solar radiation has a dominant role in the ice season heat balance, and therefore statistical degree-day models do not work there. Rather, the level of solar radiation determines how low the air temperature must go before ice formation is possible. Under similar temperature conditions, ice season is longer in a high latitude than in high altitude.
The dependence of the freezing and breakup dates on the latitude in East Europe was already presented [1]. Lake sites were at lower levels than 0.5 km above the sea level, and there were clear trends for both dates ( Figure 6). Here, we added our sites from the Asian side of the continent. In addition, Lake El'gygytgyn (67°30′ N 172° E) from Northeast Siberia has been added in the plot [34]. This lake is a deep crater lake with an area of 14 km 2 . In principle, at high altitudes, the trend lines are shifted, as suggested by our new data points, but the eastward winter cooling effects are also included. The Lake Baikal breakup date fits in the picture but freezing date is quite late, because the large mixed-layer depth in Baikal delays the freezing date but does not affect the breakup. The dependence of the freezing and breakup dates on the latitude in East Europe was already presented [1]. Lake sites were at lower levels than 0.5 km above the sea level, and there were clear trends for both dates ( Figure 6). Here, we added our sites from the Asian side of the continent. In addition, Lake El'gygytgyn (67 • 30 N 172 • E) from Northeast Siberia has been added in the plot [34]. This lake is a deep crater lake with an area of 14 km 2 . In principle, at high altitudes, the trend lines are shifted, as suggested by our new data points, but the eastward winter cooling effects are also included. The Lake Baikal breakup date fits in the picture but freezing date is quite late, because the large mixed-layer depth in Baikal delays the freezing date but does not affect the breakup.  Analytic models were presented for the freezing date and breakup date forced by a linearized surface heat budget, heat flux from below the mixed layer, the mixed layer depth, and freezing point of lake water. The linearized heat budget was expressed as = + − including two climate parameters = and ≈ . (see Appendix A). These form a temperature shift parameter = / . Together with , the local air temperature cycle defines the timing of the freeze-up and melting periods. The former depends also on the lake depth and the latter on the ice thickness.
The analytic model can be employed to examine the climate change influence. For a change ∆ in air temperature, the first-order approximation is that the parameters and do not change, but the climate impact comes from the shift of the annual tem- Analytic models were presented for the freezing date and breakup date forced by a linearized surface heat budget, heat flux from below the mixed layer, the mixed layer depth, and freezing point of lake water. The linearized heat budget was expressed as Q 0 = K 0 + K 1 (T a − T 0 ) including two climate parameters K 0 = K 0 (t) and K 1 ≈ const. (see Appendix A). These form a temperature shift parameter T R = K 0 /K 1 . Together with K 0 , the local air temperature cycle defines the timing of the freeze-up and melting periods. The former depends also on the lake depth and the latter on the ice thickness.
The analytic model can be employed to examine the climate change influence. For a change ∆T a in air temperature, the first-order approximation is that the parameters K 0 and K 1 do not change, but the climate impact comes from the shift of the annual temperature cycle. First, the condition of freeze-up, τ f > 2τ L or the cold period, must be at least twice the memory of the lake, and becomes tighter with ∆T a > 0. Following the model of Equation (7), the freezing date would be changed by α −1 ∆T a , where α is the atmospheric cooling rate in fall. e.g., in northeastern Europe, α ∼ 5 • C month −1 , and the freezing date would change by 6 days for ∆T a = 1 • C. The breakup date would be changed by β −1 ∆T a , where β is atmospheric warming date in spring, β ∼ α. A secondary effect to the breakup date comes from less ice to melt (Equation (10)), since temperature change implies change in freezing degree-days S, and consequently ice thickness changes as δh ∝ δS/h. The length of the melting period is proportional to √ h max , and therefore the sensitivity of the breakup date to ice thickness is large only for thin ice.
The quality of ice cover can be categorized into weak and stable, where the latter refers to landfast ice cover thick enough for traffic and various human activities. The transition lies between 20-30 cm thickness, apart from the soft and porous ice cover in the melting season. In addition to ice phenology, the climatology of ice quality and thickness are also key characteristics for the society, which will be examined in a future study based on modelling and remote sensing.
Many ice phenology time series papers have shown that the climate warming in recent decades has had similar cuts of some days in the ice season in its both ends, which can now be understood from the analytic model [9,15]. Even though the model looks crude, the climate change prediction is a small differential, and one would expect these predictions to be robust.  The incoming terrestrial radiation from the atmosphere is a purely external factor, while the sensible heat flux is proportional to the temperature difference T a − T 0 and thus the proportionality coefficient goes to K 1 .
The outgoing terrestrial radiation and the latent heat flux contribute to both coefficients. The former is split into two parts by using the binomial formula, and then the net terrestrial radiation is where ε 0 and ε a = ε a (T a , R, N) are the emissivities of the surface and atmosphere, respectively, and σ is the Stefan-Boltzmann constant. The truncated formula is accurate to about 1% when |T a − T 0 | < 10 • C and T a ∼ 0 • C. The latent heat exchange is split using the Clausius-Clapeyron equation for the saturation water vapour pressure e s and assuming that humidity is at the saturation level at the surface, e o = e s (T 0 ). The latent heat flux is: where ρ a is air density, L E is latent heat of evaporation (water surface) or sublimation (ice surface), C E is turbulent transfer coefficient of moisture, p a is air pressure, and R v is the gas constant of dry air. Combining all together, the coefficients K 0 and K 1 become These expressions look complicated, but the resulting climatological values are rather smooth with seasonal variation coming from solar radiation and latent heat flux.