Next Article in Journal
Novel Quantitative Evaluation of Biotreatment Suitability of Wastewater
Next Article in Special Issue
Spatial and Temporal Variation Characteristics of Glacier Resources in Xinjiang over the Past 50 Years
Previous Article in Journal
A Novel Radial Basis Function Approach for Infiltration-Induced Landslides in Unsaturated Soils
Previous Article in Special Issue
Prediction of Ice-Resistance Distribution for R/V Xuelong Using Measured Sea-Ice Parameters
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Ice Phenology in Eurasian Lakes over Spatial Location and Altitude

1
Institute of Atmospheric and Earth Sciences, University of Helsinki, 00014 Helsinki, Finland
2
Key Laboratory of Land Surface Process and Climate Change in Cold and Arid Regions, Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou 730000, China
*
Author to whom correspondence should be addressed.
Water 2022, 14(7), 1037; https://doi.org/10.3390/w14071037
Submission received: 1 March 2022 / Revised: 20 March 2022 / Accepted: 23 March 2022 / Published: 25 March 2022
(This article belongs to the Special Issue Sea, River, Lake Ice Properties and Their Applications in Practices)

Abstract

:
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 influence 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 fits. 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 first-order tool for climate sensitivity of ice seasons. The freezing date and breakup date both change by around five days per one-degree shift in air temperature away from the climatological ice margin; however, at this margin, the sensitivity is higher.

1. 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 maximum density. Then, after freeze-up, the oxygen storage in the cold lake is high enough to influence the oxygen content until spring turnover [5].
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.
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.
Water 14 01037 g001
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 low-altitude 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.

2. Materials and Methods

2.1. 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 F 0.5 is obtained from P ( t F t F 0.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].

2.2. 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 (1980–2020) in Eurasia was otained from NCEP/NCAR Reanalysis 1, 2.5° × 2.5° latitude–longitude global grid.

2.3. 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:
d ( T w H ) d t = Q 0 + Q b ρ w c w
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:
d h d t = ( Q 0 + Q w ) ρ i L f 0
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 m a x . To predict the breakup date, this thickness needs to be known. It can be estimated as h m a x = a 2 S + b 2 b , where a 2 3   cm   ( ° C   d ) 1 and b 10 cm are the model parameters; S is the cumulative temperature of freezing-degree-days [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):
Q 0 = K 0 + K 1 ( T a T 0 )
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 Q b = K b ( T b T w ) , where K b is the heat transfer coefficient and T b is the temperature below the mixed layer (or lake bottom).
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.

3. 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).
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].
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)).

4. Analytic Model

4.1. Freezing Date

The slab model (1) with the surface heat flux (3) is first written as:
d T w d t = λ a ( T a + T R T w ) + λ b ( T b T w ) ; λ a = K 1 ρ w c w H ,     λ b = K b ρ w c w H ,   T R = K 0 K 1
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:
T w ( t ) = t e λ ( t τ ) [ λ a ( T a + T R ) + λ b T b ] d τ
where λ = λ a + λ b , and τ L = λ 1 H is the response time of the lake with the proportionality coefficient of about 2 d m−1.
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 ) = c t 2 + T ^ a , where T ^ 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:
c τ L 2 < { T ^ a [ T f T R λ b λ a ( T b T f ) ] } = ( T ^ a T c )
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 as T ^ a < T c . Then from T a ( t ) = T c the shape factor c is obtained as c τ f 2 = 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 is T ^ 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 have T ^ 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 :
t F = τ L + 1 α [ T R + λ b λ a ( T b 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 2 H   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.

4.2. 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:
h m a x h ( t ) = 1 ρ i L f t M t [ K 0 + K 1 ( T a T f ) + Q w ] d τ
where h m a x = 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 Figure 2 and Figure 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 t = t M of the equation K 0 ( t ) + K 1 [ T a ( t ) T f ] = 0 :
T a ( t M ) + T R ( t M ) = T f
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 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:
1 2 ( γ + K 1 β ) ( t B t M ) 2 = ρ L f h m a x
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 Figure 2 and Figure 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:
t B = t M + 2 ρ L f γ + K 1 β h m a x t M + ( A Δ h m a x ) d
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.

5. 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.
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 km2. 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 Q 0 = K 0 + K 1 ( T a T 0 ) including two climate parameters K 0 = K 0 ( t )   and   K 1 c o n s t .   (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 m a x , 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.

Author Contributions

M.L. put together the data, performed the data analysis and wrote the original draft. L.W. contributed to the study conception and design, collected the Chinese data, and took part in the manuscript write-up. All authors have read and agreed to the published version of the manuscript.

Funding

The work was supported by the National Key Research and Development Program of China (2019YFE0197600), CAS “Light of West China” Program (E129030101, Y929641001), PIFI Fellowship of CAS 2021, the Academy of Finland project 333899, the National Natural Science Foundation of China (41930759, 41975081), and the Opening Research Foundation of the Key Laboratory of Land Surface Process and Climate Change in Cold and Arid Regions, Chinese Academy of Sciences (LPCC2019005).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The climate data from Jokioinen can be obtained by contacting M.L. The China meteorological forcing dataset can be downloaded from http://data.tpdc.ac.cn (accessed on 1 December 2021); NCEP/NCAR Reanalysis data can be downloaded from https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.pressure.html (accessed on 1 December 2021).

Acknowledgments

Two anonymous referees are thanked for their constructive comments.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Linearized Atmospheric Heat Flux

Linear solar and atmospheric heat flux is expressed as Q a = K 0 + K 1 ( T a T 0 ) , where T a and T 0 are the air temperature and surface temperature, respectively, and the coefficients K 0 and K 1 do not depend explicitly on T 0 . It is obtained as follows: The common formulations for the terrestrial radiation and turbulent fluxes are based on the standard weather variables: air temperature, relative humidity R , wind speed U a , and cloudiness N . The net solar radiation is ( 1 α ) Q s , where Q s is incident solar radiation and α is albedo. 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
Q n L = ε 0 σ ( ε a T a 4 T 0 4 ) ε 0 σ ( 1 ε a ) T a 4 + 4 ε 0 σ T a 3 ( T a T 0 )
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:
Q e = ρ a L E C E 0.622 e s ( T a ) p a [ ( 1 R ) L E R v T a 2 ( T a T 0 ) ] U a
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
K 0 = ( 1 α ) Q s ε 0 σ ( 1 ε a ) T a 4 ρ a L E C E 0.622 e s ( T a ) p a ( 1 R ) U a
K 1 = 4 ε 0 σ T a 3 + ρ a c a C H U a + ρ a L E C E 0.622 e s ( T a ) p a L E R v T a 2 U a
These expressions look complicated, but the resulting climatological values are rather smooth with seasonal variation coming from solar radiation and latent heat flux.

References

  1. Kirillin, G.; Leppäranta, M.; Terzhevik, A.; Granin, N.; Bernhardt, J.; Engelhardt, C.; Efremova, T.; Golosov, S.; Palshin, N.; Sherstyankin, P.; et al. Physics of seasonally ice-covered lakes: A review. Aquat. Sci. 2012, 74, 659–682. [Google Scholar] [CrossRef]
  2. Leppäranta, M. Freezing of Lakes and Evolution of Their Ice Cover, 2nd ed.; Springer: Berlin/Heidelberg, Germany, 2022. [Google Scholar]
  3. Korhonen, J. Long-term changes in lake ice cover in Finland. Nord. Hydrol. 2006, 37, 347–363. [Google Scholar] [CrossRef]
  4. Leppäranta, M. Interpretation of statistics of lake ice time series for climate variability. Nord. Hydrol. 2014, 45, 673–683. [Google Scholar] [CrossRef]
  5. Salonen, K.; Leppäranta, M.; Viljanen, M.; Gulati, R.D. Perspectives in winter limnology: Closing the annual cycle of freezing lakes. Aquat. Ecol. 2009, 43, 609–616. [Google Scholar] [CrossRef]
  6. Bates, R.E.; Bilello, M.A. Defining the Cold Regions of the Northern Hemisphere; Technical Report 178; US Army Cold Regions Research and Engineering Laboratory: Hanover, NH, USA, 1966. [Google Scholar]
  7. Magnuson, J.J.; Robertson, D.M.; Benson, B.J.; Wynne, R.H.; Livingstone, D.M.; Arai, T.; Assel, R.A.; Barry, R.G.; Card, V.; Kuusisto, E. Historical trends in lake and river ice cover in the Northern Hemisphere. Science 2000, 289, 1743–1746. [Google Scholar] [CrossRef] [Green Version]
  8. Yao, X.; Long, L.I.; Zhao, J.; Sun, M.; Jing, L.I.; Gong, P.; Lina, A.N. Spatial-temporal variations of lake ice phenology in the Hoh Xil region from 2000 to 2011. J. Geogr. Sci. 2015, 31, 70–82. [Google Scholar] [CrossRef]
  9. Karetnikov, S.; Leppäranta, M.; Montonen, A. A time series of over 100years of ice seasons on Lake Ladoga. J. Great Lakes Res. 2017, 43, 979–988. [Google Scholar] [CrossRef]
  10. Sharma, S.; Blagrave, K.; Magnuson, J.J.; O’Reilly, C.M.; Oliver, S.; Batt, R.D.; Magee, M.R.; Straile, D.; Weyhenmeyer, G.A.; Winslow, L. Widespread loss of lake ice around the Northern Hemisphere in a warming world. Nat. Clim. Chang. 2019, 9, 227–231. [Google Scholar] [CrossRef]
  11. Kropacek, J.; Maussion, F.; Chen, F.; Hoerz, S.; Hochschild, V. Analysis of ice phenology of lakes on the Tibetan Plateau from MODIS data. Cryosphere 2013, 7, 287–301. [Google Scholar] [CrossRef] [Green Version]
  12. Cai, Y.; Ke, C.Q.; Duan, Z. Monitoring ice variations in Qinghai Lake from 1979 to 2016 using passive microwave remote sensing data. Sci. Total Environ. 2017, 607–608, 120–131. [Google Scholar] [CrossRef]
  13. Cai, Y.; Ke, C.Q.; Li, X.; Zhang, G.; Duan, Z.; Lee, H. Variations of lake ice phenology on the Tibetan Plateau from 2001 to 2017 based on MODIS data. J. Geophys. Res. 2019, 124, 825–843. [Google Scholar] [CrossRef]
  14. Mishra, V.; Cherkauer, K.A.; Bowling, L.C.; Huber, M. Lake ice phenology of small lakes: Impacts of climate variability in the Great Lakes region. Glob. Planet. Chang. 2011, 76, 166–185. [Google Scholar] [CrossRef]
  15. Bernhardt, J.; Engelhardt, C.; Kirillin, G.; Matschullat, J. Lake ice phenology in Berlin-Brandenburg from 1947–2007: Observations and model hindcasts. Clim. Chang. 2012, 112, 791–817. [Google Scholar] [CrossRef]
  16. Wang, C.; Shirasawa, K.; Leppäranta, M.; Ishikawa, M.; Huttunen, O.; Takatsuka, T. Solar radiation and ice heat budget during winter 2002–2003 in Lake Pääjärvi, Finland. Verh. Int. Verein Limnol. 2005, 29, 414–417. [Google Scholar] [CrossRef]
  17. Leppäranta, M.; Lindgren, E.; Shirasawa, K. The heat budget of Lake Kilpisjrvi in the Arctic tundra. Hydrol. Res. 2016, 48, 969–980. [Google Scholar] [CrossRef]
  18. Mironov, D.; Ritter, B.; Schulz, J.-P.; Buchhold, M.; Lange, M.; Machulskaya, E. Parameterisation of sea and lake ice in numerical weather prediction models of the German Weather Service. Tellus A 2012, 64, 17330. [Google Scholar] [CrossRef] [Green Version]
  19. Huang, W.F.; Cheng, B.; Zhang, J.R.; Zhang, Z.; Vihma, T.; Li, Z.J.; Niu, F.J. Modeling experiments on seasonal lake ice mass and energy balance in the Qinghai-Tibet Plateau: A case study. Hydrol. Earth Syst. Sci. 2019, 23, 2173–2186. [Google Scholar] [CrossRef] [Green Version]
  20. Song, S.; Li, C.Y.; Shi, X.H.; Zhao, S.N.; Tian, W.D.; Li, Z.J.; Bai, Y.L.; Cao, X.W.; Wang, Q.K.; Huotari, J.; et al. Under-ice metabolism in a shallow lake in a cold and arid climate. Freshw. Biol. 2019, 64, 1710–1720. [Google Scholar] [CrossRef]
  21. Jevrejeva, S.; Drabkin, V.; Kostjukov, J.; Lebedev, A.; Leppäranta, M.; Mironov, Y.U.; Schmelzer, N.; Sztobryn, M. Baltic Sea ice seasons in the twentieth century. Clim. Res. 2004, 25, 217–227. [Google Scholar] [CrossRef] [Green Version]
  22. Leppäranta, M.; Lindgren, E.; Wen, L.; Kirillin, G. Ice cover decay and heat balance in Lake Kilpisjrvi in Arctic tundra. J. Limnol. 2019, 78, 163–175. [Google Scholar] [CrossRef] [Green Version]
  23. Rodhe, B. On the relation between air temperature and ice formation in the Baltic. Geogr. Ann. 1952, 34, 175–202. [Google Scholar]
  24. Barnes, H.T. Ice Engineering; Montreal Renouf Publishing Co.: Montreal, QC, Canada, 1928; pp. 1–364. [Google Scholar]
  25. Mueller, D.R.; Hove, P.V.; Antoniades, D.; Jeffries, M.O.; Vincent, W.F. High Arctic lakes as sentinel ecosystems: Cascading regime shifts in climate, ice cover, and mixing. Limnol. Oceanogr. 2009, 54, 2371–2385. [Google Scholar] [CrossRef]
  26. Zhang, T.; Jeffries, M.O. Modeling interdecadal variations of lake-ice thickness and sensitivity to climatic change in northernmost Alaska. Ann. Glaciol. 2000, 31, 339–347. [Google Scholar] [CrossRef] [Green Version]
  27. Ashton, G.D. Dynamics of Snow and Ice Masses; Academic Press: New York, NY, USA, 1980; pp. 261–304. [Google Scholar]
  28. Barańczuk, K.; Barańczuk, J. The ice regime of Lake Ostrzyckie (Kashubian Lakeland, northern Poland). Limnol. Rev. 2019, 19, 105–112. [Google Scholar] [CrossRef] [Green Version]
  29. Huang, W.; Zhang, Z.; Li, Z.; Leppäranta, M.; Arvola, L.; Song, S.; Lin, Z. Under-ice dissolved oxygen and metabolism dynamics in a shallow lake: The critical role of ice and snow. Water Resour. Res. 2021, 57, e2020WR027990. [Google Scholar] [CrossRef]
  30. Mietus, M. The Climate of the Baltic Sea Basin; WMO/TD-No. 933; MMROA Report-No. 41; World Meteorological Organization: Geneva, Switzerland, 1998. [Google Scholar]
  31. Yang, F.; Li, C.Y.; Shi, X.H.; Zhao, S.N.; Hao, Y.Z. Impact of seasonal ice structure characteristics on ice cover impurity distributions in Lake Ulansuhai. J. Lake Sci. 2016, 28, 455–462. [Google Scholar]
  32. Kouraev, A.V.; Semovski, S.V.; Shimaraev, M.N.; Mognard, N.M.; Legrésy, B.; Rémy, F. The ice regime of Lake Baikal from historical and satellite data: Relationship to air temperature, dynamical, and other factors. Limnol. Oceanogr. 2007, 52, 1268–1286. [Google Scholar] [CrossRef]
  33. Qiu, Y.; Guo, H.; Ruan, Y.; Fu, X.; Shi, L.; Tian, B. A Dataset of Microwave Brightness Temperature and Freeze-Thaw for Medium-to-Large Lakes over the High Asia Region 2002–2016; Science Data Bank: Beijing, China, 2017. [Google Scholar] [CrossRef]
  34. Nolan, M.; Liston, G.; Prokein, P.; Brigham-Grette, J.; Sharpton, V.L.; Huntzinger, R. Analysis of lake ice dynamics and morphology on Lake El’gygytgyn, NE Siberia, using synthetic aperture radar (SAR) and Landsat. J. Geophys. Res. 2002, 107, ALT 3-1–ALT 3-12. [Google Scholar] [CrossRef]
Figure 2. The parameters K 0   ( a )   and   K 1 (b) of the linear surface heat transfer formula Q 0 = K 0 + K 1 ( T a T 0 ) 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).
Figure 2. The parameters K 0   ( a )   and   K 1 (b) of the linear surface heat transfer formula Q 0 = K 0 + K 1 ( T a T 0 ) 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).
Water 14 01037 g002
Figure 3. Annual average temperature cycle in the lakes under study. The black dots show the freezing and breakup dates.
Figure 3. Annual average temperature cycle in the lakes under study. The black dots show the freezing and breakup dates.
Water 14 01037 g003
Figure 4. The temperature parameter T R in Lake Pääjärvi and Ngoring Lake. This temperature indicates how much solar radiation compensates for air temperature forcing in the annual cycle of the lakes. Data from Figure 2.
Figure 4. The temperature parameter T R in Lake Pääjärvi and Ngoring Lake. This temperature indicates how much solar radiation compensates for air temperature forcing in the annual cycle of the lakes. Data from Figure 2.
Water 14 01037 g004
Figure 5. Mean air temperature (1980–2020) in Eurasia in December (a) and April (b).
Figure 5. Mean air temperature (1980–2020) in Eurasia in December (a) and April (b).
Water 14 01037 g005
Figure 6. Long-term averages of the freezing date and breakup date in several Eurasian lakes. Original plot with white squares and black squares from Kirillin et al. [1]. Blue ovals have been added in this work for Asian lakes.
Figure 6. Long-term averages of the freezing date and breakup date in several Eurasian lakes. Original plot with white squares and black squares from Kirillin et al. [1]. Blue ovals have been added in this work for Asian lakes.
Water 14 01037 g006
Table 1. Study lakes in the calibration of analytical models. The data are from local sources.
Table 1. Study lakes in the calibration of analytical models. The data are from local sources.
LakeLatitude and LongitudeAltitude
m
Area
km2
Depth
m
Lake TypeClimate Zone
Kilpisjärvi69° 03′ N 20° 50′ E4733719.5Oligotrophic Arctic tundra
Kallavesi62° 46′ N 27° 47′ E824789.7OligotrophicBoreal
Pääjärvi61° 04′ N 25° 08′ E1031315.0Meso-oligotr.Boreal
Baikal53° 18′ N 108° 00′ E45631,772758OligotrophicBoreal
Szczecin lagoon53° 49′ N 14° 08′ E06874.6EutrophicBoreal
Wuliansuhai 40° 56′ N 108° 52′ E10192331.1EutrophicCold arid
Qinghai36° 53′ N 100° 11′ E3260423721.0OligotrophicTibetan
Gyaring34° 56′ N 97° 17′ E42945268.9OligotrophicTibetan
Ngoring34° 54′ N 97° 42′ E427261017.6OligotrophicTibetan
Table 2. Ice phenology in the selected Eurasian lakes. Notation ± nn refers to the standard deviation in days (d). The ice phenology data are based on the periods of 1960–2020 (Finland), 1991–2020 (China), 1896–2000 (Szcezin logoon), and 1960–2000 (Baikal), and the air temperature data are for the period of 1981–2010.
Table 2. Ice phenology in the selected Eurasian lakes. Notation ± nn refers to the standard deviation in days (d). The ice phenology data are based on the periods of 1960–2020 (Finland), 1991–2020 (China), 1896–2000 (Szcezin logoon), and 1960–2000 (Baikal), and the air temperature data are for the period of 1981–2010.
LakeFreezing DateDecember Mean Air TemperatureApril Mean Air TemperatureBreakup Date
Kilpisjärvi8 Nov. ± 8 d−12 °C−5 °C18 Jun. ± 7 d
Kallavesi30 Nov. ± 14 d−6 °C+1 °C8 May ± 8 d
Pääjärvi13 Dec. ± 16 d−5 °C+4 °C5 May ± 6 d
Szczecin lagoon [21,30]15 Jan.−0.5 °C+6 °C1 Mar.
Wuliansuhai [31]7 Dec.−10 °C+4 °C31 Mar.
Baikal [32]10 Jan.−18 °C+1 °C1 May.
Gyaring [33]30 Nov. ± 21 d−18 °C−2 °C21 April ± 17 d
Ngoring [33]15 Dec. ± 10 d−17 °C−2 °C1 May ± 8 d
Qinghai [33]29 Dec. ± 4 d−10 °C+4 °C1 April ± 7 d
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Leppäranta, M.; Wen, L. Ice Phenology in Eurasian Lakes over Spatial Location and Altitude. Water 2022, 14, 1037. https://doi.org/10.3390/w14071037

AMA Style

Leppäranta M, Wen L. Ice Phenology in Eurasian Lakes over Spatial Location and Altitude. Water. 2022; 14(7):1037. https://doi.org/10.3390/w14071037

Chicago/Turabian Style

Leppäranta, Matti, and Lijuan Wen. 2022. "Ice Phenology in Eurasian Lakes over Spatial Location and Altitude" Water 14, no. 7: 1037. https://doi.org/10.3390/w14071037

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop