Lakes and their ice cover affect the local climate and weather by modulating the temperature, wind, humidity and precipitation (Bonan, 1995; Ellis and Johnson, 2004; Rouse et al., 2008; Thiery et al., 2015). Essentially, lake ice reduces heat and moisture transfer between the atmosphere and the lake water through its insulating effect from winter to spring. If the lake ice is covered with snow, this influence is further strengthened due to the lower net incoming solar radiation by the higher albedo and the lower thermal conductivity of the latter compared to the former. Therefore, the ice cover duration and the growth processes of lake ice are important for understanding the local/regional winter climate.
So far, seasonally ice-covered lakes have been studied mainly for North American and northern European lakes, and little attention has been paid to lakes at mid-latitudes in other regions. To understand them on a global scale, it is necessary to understand the properties and the formation processes of lake ice at mid-latitudes where growth conditions might be quite different from those at high latitudes. Although the lake ice properties at mid-latitudes were investigated for lakes in Hokkaido, Japan, by Kawamura and Ono (1978) and Muguruma and Kikuchi (1963), their studies were focused on crystallographic structure. The climate at mid-latitudes except arid regions, e. g. northern China, is characterized by relatively moderate air temperature, a considerable amount of snow and abundant solar radiation even in winter. Lake ice is commonly composed of congelation ice (CI) formed from the bottom freezing of lake water, and snow ice (SI) formed from water flooding into the snow on ice and subsequent freezing (Michel and Ramseier, 1971). It is well known that snow reduces the growth of CI thickness due to its low thermal conductivity, while it can also enhance the growth of ice thickness through SI formation (e.g. Leppäranta, 1983, 2015; Saloranta, 2000; Yang et al., 2012). Thus, to understand the lake ice growth processes at mid-latitudes, it is especially important to understand the role of snow on lake ice growth. This can apply to lake ice at high latitudes because an increase in the snowfall there is predicted in the future, associated with global warming (Dai, 2006; Simmons et al., 2010; Willett et al., 2012).
As a preliminary step, Ohata et al. (2016) investigated the lake ice properties at Lake Abashiri, Hokkaido, Japan, located at mid-latitudes (44.0 N 144.2 E) based on observations throughout the 2012/13 winter. The results showed that (1) the SI layer accounted for 29% to 73% of the total ice thickness, a much greater fraction than that for high latitude lakes (10–40%, Finland, Leppäranta, 2009), (2) a thermodynamic model including both SI and CI formation processes successfully reproduced the observed thickness. Based on observations and numerical experiments with this model, they presented a hypothesis on the role of snow in the lake ice thickening processes. The hypothesis is that snow plays a role in moderating the variability of the ice thickness caused by differences in meteorological conditions through the insulation effect and SI formation. Since these results were derived from only one winter observation, however, additional observations and a longer data period containing a wider variety of meteorological conditions are needed for model validation.
Several recent studies have reported that lake ice phenology (dates of freeze-up and break-up, ice cover duration) has changed significantly worldwide (IPCC 2013). Conventionally, the freeze-up date is defined as the first date when a lake is completely ice-covered, whereas the break-up date is the first date when it is completely ice-free (Kirillin et al., 2012). Magnuson et al. (2000) showed that on average the freeze-up date of 26 lakes and rivers in the Northern Hemisphere became later at a rate of 0.58 d/decade, and the break-up date earlier at a rate of 0.65 d/decade, over the period 1846–1995. According to more recent studies, ice cover duration reduction appears to be accelerating. For example, Benson et al. (2012) showed that the freeze-up date of 75 lakes in the Northern Hemisphere became later at a rate of 1.1 d/decade, and the break-up date earlier at a rate of 0.9 d/decade over the period 1855/56–2004/05, while the freeze-up date became later at a rate of 1.6 d/decade and the break-up date earlier at a rate of 1.9 d/decade over the period 1975/76–2005/06. For mid-latitude lake ice, there have been few phenology studies except at Lake Suwa, Japan (36 °N, 138 °E), where the freeze-up date was shown to have become later at a rate of 0.2 d/decade over the 550-year record (Magnuson et al., 2000). Due to a lack of ice thickness data and analysis of ice structure, however, factors involved in such trends remain unknown.
The purpose of this study is to verify the thermodynamic model developed for predicting the lake ice growth at Lake Abashiri and the hypothesis presented by Ohata et al. (2016), using the additional data obtained from the follow-up campaigns in two winters (2014/15, 2015/16), to clarify the role of snow in lake ice thickening processes, and then to show the expected properties of the ice phenology for comparison with other regions with the model applied for the long-term meteorological data records at Abashiri. As for the validation of the model, we paid special attention to the validity of the snow compression rate (β) and the timescale of the conversion of flooded snow (slush) layer to snow ice because these are critical assumptions in the model. For this purpose, during the observations in two winters (2014/15 and 2015/16), the freeze-up conditions were directly observed from a small boat in December, and the thickness measurements of snow, CI and SI were conducted at fixed sites once per month from January to March. The process of SI formation was monitored with a set of thermistors installed in the ice and snow.
To further examine the role of snow in the lake ice growth under various meteorological conditions, numerical model experiments were carried out for the three winters (2012/13, 2014/15 and 2015/16), and then the analysis period was extended using meteorological data records at Abashiri for the past 16 years (2000/01–2015/16) during which the freeze-up dates can be determined with minimal uncertainty from NASA Moderate Resolution Imaging Spectroradiometer (MODIS) satellite images. Finally, to examine the ice phenology for Lake Abashiri, the interannual variability of ice thickness and dates of break-up and freeze-up were calculated for the past 55 years (1961/62–2015/16) using our model, using nearby staffed weather station air temperature as proxy data. The results will be compared with phenology data from elsewhere in the Northern hemisphere (Benson et al., 2012), and meteorological factors responsible for the interannual trends will be discussed.
Although this study is limited to the properties of lake ice at mid-latitudes in Japan, the results are expected to apply to other lakes, such as the Great Lakes of North America, where meteorological conditions are relatively similar. In addition, since SI formation is regarded as an important sea ice growth process, especially in the Antarctic regions as well as in the Sea of Okhotsk, north of Hokkaido, Japan, the results are also expected to be applicable to this domain (Toyota et al., 2007, 2011).
Field observations were carried out at Lake Abashiri, located in the north-eastern part of Hokkaido, Japan (43.97° N, 144.17° E). The water depth is ~16 m at its maximum and 7 m on average, and the lake covers 32.3 km2, with dimensions of ~10 km long and ~3 km wide. Three observation sites were set along a line running from upstream (south-west) to downstream (north-east) on Lake Abashiri.
Observation details and results for the 2012/13 season were described in Ohata et al. (2016). To briefly summarize these, analysis of ice blocks collected at six sites on the lake showed that the ice was composed of two distinct layers: whitish SI and translucent CI, and SI occupied as much as 29–73 % of the total ice thickness. This is much thicker than for lake ice at high latitudes (e.g. 10–40%, in Finland; Leppäranta, 2009). To examine the interannual variability of these results, follow-up in-situ observations were carried out at sites 1, 3 and 6 once per month from December to March in the 2014/15 and 2015/16 ice seasons (Fig. 1). During these two ice seasons, we began by examining the freeze-up conditions with a small boat on the lake in December to confirm the freeze-up date. For the period January to March, the sampling of ice blocks, lake water and snow on ice was conducted at each site for structure analysis, and the vertical profiles of salinity and temperature of lake water from 0.15 m below the surface to the bottom of the lake were also observed at the same places with a conductivity/ temperature/depth (CTD) profiler (CastAway-CTD from Son Tek/YSI Inc.).
In addition, to understand the SI formation processes, the temperatures of air, snow, ice and water were monitored at site 3 from 29 January to 11 March in 2016, using the 12 thermo-recorders (RT-32S from ESPEC MIC Corporation), with a measurement accuracy of 0.1 °C. The evolution of the temperature profile is expected to reveal the conversion process of snow or the slush layer to SI because of their significant differences in thermal conductivity. Before the measurement, the sensors were calibrated in an ice water-filled bath in which temperature was set to 0 °C in equilibrium. Among the 12 thermistor sensors, 2 were installed in the air, 5 in snow, 1 at the snow/ice interface, 2 in ice and 2 in the water beneath the ice, as shown in Fig. 2. On the set-up day (29 January 2016), snow depth was 31 cm, containing a 22-cm thick slush layer composed of a flooded snow layer (15 cm) and the resultant wet snow layer (7 cm) just above it. The water level was at the top of the flooded snow. The sensors were set at an interval of 5 cm in the snow layer (4 in slush and 1 in dry snow).
Regarding model input parameters, we used meteorological data-sets recorded at the Abashiri Meteorological Observatory (AMO; 38 m a.s.l.) and Memanbetsu Airport (33 m a.s.l.), which are located about 10 km north-east and 9 km south, respectively, from the central site of Lake Abashiri. While the data-sets of Ta (air temperature), Va (wind velocity) and Hs (snow depth on land) are available at both stations on a daily basis, Qs (incoming solar radiation at the surface), C (cloudiness in tenths), r (relative humidity) and p (surface air pressure) were observed only at AMO. Therefore, to calculate ice thickness growth at site 3, Ta, Hs and Va calculated as the mean of the two meteorological station data, and Qs, C, r and p were given by the data observed at AMO.
To validate the ice thickness calculated by our model, daily lake ice thickness records observed at the near-shore area of nearby Yobito (Fig. 1b), recorded by the Abashiri Tourism Association, were also used in this study. This observation supports fishing activities in winter and the records cover the period from 2005/06 to 2015/16. Although the difference in measurement site might hamper exact validation of our model, the data are expected to help us check the yearly variation of ice thickness.
In our earlier work, a 1-D thermodynamic quasi-steady-state model was developed to calculate the ice thickness evolution at the centre of the lake (site 3), based on observations conducted during the 2012/13 winter season (Ohata et al., 2016). In this section, we describe the contents of this model briefly. The freeze-up date was determined from MODIS satellite images, photos taken from the top of the nearby mountain (Tento-san, Fig. 1b) and the surface heat budget. Although ideally it should be determined from the calculation of the heat contents at the surface layer, taking into account the effect of mixing with lower water layers, we took this method for a practical reason due to lack of data. Lake ice was assumed to be a horizontal uniform slab, composed of SI and CI, with a snow cover. To calculate the growth of CI, we followed the method of Maykut (1982), using meteorological data-sets obtained from the nearby meteorological observatories. The time step was set to one day. Given that there is still uncertainty in the freezing process and parameters of the slush layer, the following three assumptions were introduced to calculate the SI growth rate: (i) the slush layer freezes to form SI on the time scale of a few days (i.e. short enough to neglect details of the freezing process), then the freeboard recovers from negative to zero; (ii) snow is compressed in the transformation into SI by a compression rate β, which is defined as the flooded snow layer thickness divided by the produced SI layer thickness; (iii) β is kept constant throughout the season. Here, β is a parameter originally introduced by Leppäranta and Kosloff (2000). In this model, proper determination of the value of the tuning parameter β is essential to reproduce the ice thickness evolution. Judging from the fit of the calculated SI thickness to the observational data of the 2012/13 season, β was determined to be 2.0. The other parameters and constants used in the model are listed in Table 1. One of the purposes of this study is to examine the effectiveness of this model from additional observations.
Meteorological conditions were quite different between these three ice seasons. Time series of air temperature and snow depth from December to March are shown in Fig. 3. When we look at the mean air temperature and snow depth (given by the snow thickness increment from freeze-up date) from freeze-up date to 20 February to characterize the meteorological conditions during the ice growth period, it is found that the temperature was higher and the snow deeper in the 2014/15 winter than in the two other winters (Table 2). In Fig. 3, also note that air temperature was relatively high and the freeze-up was relatively late in December 2015/16. On the other hand, snow depth on land around the lake was relatively deep on the freeze-up date in the 2012/13 winter, while quite small in the 2014/15 and 2015/16 winters. Moreover, the snow depth in 2015/16 is characterized by relatively small accumulation rate for about three weeks after freeze-up date, while several heavy snowfall events took place after the freeze-up in 2014/15.
The weather conditions at the times of observation in 2014/15 and 2015/16 were as follows: for 2014/15, the freeze-up conditions observed from the boat were calm and cloudy with an air temperature of 0–1 °C on 11 December 2014 and calm and cloudy with 1–2 °C on 11 December 2015. The weather conditions on ice in early 2015 were cloudy with −5 to 0 °C on 30 January, cloudy sometimes fair with 0–4 °C on 24 February and fair with 0–5 °C on 12 March. For 2016, conditions were fair and sometimes cloudy with −5 to −4 °C on 29 January, fair with −8 to −4 °C on 19 February and fair with −4 to 3 °C on 11 March. On all days, winds were relatively weak and the observations were conducted between 08:00 and 15:00 local time.
First, we describe how we determined the freeze-up dates from our observations. On 11 December 2014, the lake had already frozen at the central part, the upstream part and along the shore region except for a few ice-free regions near site 6 and some offshore regions from shore (Fig. 1). The ice thickness along the shore and near the site 6 was 1 to 1.5 cm with translucent colour and the water temperature just beneath the ice (about 1–2 cm beneath the ice) was +0.1 °C. Based on this result, MODIS images and the negative surface heat budget estimated from meteorological data (−47 W m−2), we regard 10 December as the freeze-up date in the 2014/15 ice season. For the 2015/16 ice season, most of the lake surface was ice-free except for translucent thin ice plates (0.1–1.2 cm thick) floating near the shore and site 6 on 11 December 2015. The surface water temperature was still somewhat high (+0.4 to 1.3 °C). Based on this result and MODIS images, taking the meteorological conditions into account, we regarded 23 December as the freeze-up date.
Next, ice thicknesses and snow depths on ice in February at each site for three winters are listed in Table 3. To check the spatial representativeness of the obtained ice thickness, the spatial variability was examined by measuring thicknesses at nine points within a 10 m × 10 m grid on 19 February 2016. The average ice thickness was 40.2 ± 2.1 cm at site 1, 41.9 ± 3.2 cm at site 3 and 47.9 ± 1.5 cm at site 6, showing that the obtained ice thickness is considered to be representative of the surrounding area.
Ice structures were commonly composed of whitish (upper) and translucent (lower) layers, and it was confirmed from the δ18O measurements that they correspond to SI and CI, respectively. The contribution of SI to total ice thickness was significant: 29–73% for 2012/13, 40–60% for 2014/15 and 28–65% for 2015/16. In considerations of a considerable amount of snow and an abundant solar radiation at this lake, while there is a large variability of thickness, this result indicates that SI makes a stronger contribution to lake ice formation at Lake Abashiri, located at mid-latitudes, than at high latitudes. The densities of CI (ρci), SI (ρsi) and snow (ρs) at site 3 were 898 ± 14 kg m−3, 821 ± 81 kg m−3 and 380 ± 31 kg m−3 from the observation on 24 February 2015, and 872 ± 16 kg m−3, 838 ± 13 kg m−3 and 243 ± 37 kg m−3 on 19 February 2016, respectively, suggesting that SI has a somewhat lower density than CI.
From the CTD measurements, it was shown that lake water was composed of two layers with significantly different salinities during the winter: an almost fresh water layer (layer thickness 5–6 m, salinity 0.23 ± 0.06 psu for depth ≲0.6 m; 0.2–10 psu ≳0.6 m) for the surface water and a brackish water layer (layer thickness ~10 m, salinity ~20 psu for depth ≳10 m). Thus, we can assume that the lake ice here formed from fresh water.
To examine the validity of the model developed for the 2012/13 season (Ohata et al., 2016), the time series of hi, hsi, hci and hs was calculated in the same way with the same value of β for the 2014/15 and 2015/16 ice seasons. Here, we denote the thickness by h, using subscripts i, si, ci and s for total ice, snow ice, congelation ice and snow, respectively. The calculated results are shown in Fig. 4, along with the results for 2012/13. It is found that overall hci, hsi and break-up dates show a good agreement with the observed values except for hsi in January 2016 and hi in March in both years. The discrepancies between model and observational results were 4 cm for hi, 3 cm for hsi, 1 cm for hci and 12 cm for hs in February and 4 d for break-up date in 2015, and 0 cm, 4 cm, −4 cm, 6 cm and 1 d, respectively, in 2016. These results are almost the same as the values of −3 cm, −4 cm, 1 cm, 9 cm and 2 d in 2013 (Fig. 4). Thus, it is suggested that this model can reproduce the observed ice thickness evolution and break-up date well to some extent under different meteorological conditions, allowing us to examine the interannual variability of these parameters.
As for the strong discrepancy (11 cm) for hi in January 2016, in the light of the fact that the slush layer in snow amounted to as thick as 22 cm, the model assumption (i) in Section 3 did not apply in this case, i.e. the flooded snow layer was not converted to SI over a few days. On the other hand, the fact that no slush layer was observed on 19 February 2016 shows that slush layer was converted to snow ice completely during this period. Thus, this condition provided us with an opportunity to examine the formation process of SI in detail from the temporal evolution of temperature profile. The details will be discussed in Section 5.1.
As for the discrepancy (15 and 7 cm) for hi in March, in the light of the fact that it is attributed mainly to the underestimation of hsi (Fig. 4) and air temperature sometimes exceeded 0 °C in March (Fig. 3), the omission of the increase in hi due to the refreezing of melted snow in the model is presumably most responsible. This also explains the reason for somewhat earlier break-up prediction in the model than the observed dates. Despite this deficiency, it is noted that the modelled break-up date differs from the observed date by only less than four days. This result suggests that break-up date may be controlled more by meteorological conditions rather than by ice thickness. This speculation may be supported by the fact the net surface heat flux estimated, using the model with meteorological conditions, kept positive after ~15 March and it reached 100–200 W m−2 for 2014/15 and 2015/16 (the mean net surface heat flux from 15 March to the break-up date; 78 W m−2 for 2015, 71 W m−2 for 2016), due to the strong solar radiation, while the net surface heat flux turned positive on 24 March 2013 (the mean net surface heat flux after 15 March; 62 W m−2). The break-up dates estimated for 2015 and 2016 were the same dates (8 April) and that for 2014 was later (11 April).
Since the validity of our model was confirmed to some extent, here we investigate the role of snow in the ice-thickening processes quantitatively with this model by examining the sensitivity of the modelled ice thicknesses to snow. For this purpose, numerical experiments for calculating the thickness evolution were conducted under the meteorological conditions of these three ice seasons for the following three cases: where (1) snow accumulates accompanied with SI formation (control run), (2) snow accumulates without SI formation (i.e. CI only) and (3) no snow accumulates. Since the meteorological conditions were quite different among the three ice seasons, as described in Section 4.1, the comparison of the results is expected to serve to check the hypothesis that snow works to moderate the ice thickness variation caused by meteorological conditions, presented by Ohata et al. (2016) which was derived essentially from different air temperature conditions. Considering that the accuracy of the modelled ice thickness is highest in February, we focus on the results for February. The results are shown in Fig. 5. As shown in the previous subsection, it is confirmed that case 1 is in good agreement with observations for all the three winters. Therefore, we examine the role of snow by comparing the results individually for each winter.
Firstly, to investigate a case with similar snow depth but different air temperature, the results obtained for the 2012/13 and 2015/16 winters are compared. Here, the freeze-up date was 12 d later in 2015/16 than in 2012/13. For case 3 (no snow), the difference in hci is (68 – 56) 12 cm (Fig. 5). This is considered to be caused by higher air temperature and later freeze-up in 2015/16. Note that for case 2 (snow included, no SI growth), hci is by 10 cm thicker in 2015/16 than in 2012/13, contrasting to case 3. This contrast is attributed to the thermal insulation effect of snow because hs at the initial ice growing period is significantly less in the 2015–2016 winter (Fig. 4). For case 1 (control run), although the difference in hci still remains as case 2, the difference in hi is reduced to only 3 cm (= 45 − 42) by the increase in hsi and the difference in hi found for case 3 (no snow) was reduced by about 75% in the presence of snow. This indicates that the difference in hi which might have been produced in the absence of snow was significantly moderated through the thermal insulation effect of snow and the SI formation. This supports the hypothesis presented by Ohata et al., (2016) which was derived from observations under similar conditions.
Next, to investigate a case of both different snow depth and air temperature, two comparisons were made; (i) between the 2012/13 and 2014/15 winters, and (ii) between the 2014/15 and 2015/16 winters. As for (i), the result for case 3 shows that the difference in hi on 20 February is 13 (= 68 − 55) cm. For case 2, hi is significantly thinner than that for case 1 and is thicker in 2014/15 than in 2012/13 due to the thermal insulation effect of snow. In this case, however, hci for case 1 was calculated to be somewhat thinner in 2014/15 than in 2012/13. For case 1, this slight difference was compensated by higher SI growth in 2015 (32 cm) than in 2013 (23 cm). As a result, hi became rather thicker in 2014/2015 than in 2012/13, although the difference in hi was reduced about by half (45 − 52 = −7 cm) compared with that for case 3 (13 cm).
As for (ii), the result for case 3 shows that the difference in hi on 20 February is only 1 cm (= 56 − 55). This small difference in spite of the lower temperature in 2015/16 is attributed to a later freeze-up date. On the other hand, for case 2, hci is predicted to be thinner by 7 cm in 2014/15 because the thermal insulation effect of snow worked effectively due to thicker snow depth in this season. This can explain the difference in hci calculated for case 1, where hci is thinner by 7 cm in 2014/15 than in 2015/16. However, hi for case 1 is calculated to become thicker via larger SI conversion in 2014/15. This suggests the importance of SI formation over the effects of thermal insulation or air temperature to determine hi in February.
Thus, it is indicated from all these results that while they support the hypothesis by Ohata et al. (2016) that snow works to moderate the variability of hi caused by differences in air temperature, SI formation plays a more essential role in determining hi. In addition, comparison of the results between 2012/13 and 2015/16 suggests that although the timing of snowfall events affects hci due to the thermal insulation effect at the initial stage, hi in February depends more on the SI growth. Schematic pictures which explain how snow works to moderate the variability of hi through the thermal insulation effect and the SI formation under different meteorological conditions are shown in Fig. 6: the same snow amount and different air temperature (Fig. 6a) and different snow amount and the same air temperature (Fig. 6b).
In this section, the results obtained in the previous section will be examined further with our model under various meteorological conditions for the past 16 years (2000/01–2015/16) when MODIS satellite images are available to determine the freeze-up and break-up dates within a few days.
The evolutions of hi, hci, and hsi were calculated with our model in the same way by inputting meteorological data records from Abashiri. When determining the freeze-up date, although MODIS images serve to detect the freeze-up state efficiently in snow-covered conditions, they do not work as well in snow-free conditions due to the relatively low albedo of lake ice. In such cases, air temperature was used to judge freeze-up date. The predicted hi for each year was compared with the ice thickness records obtained in the nearshore region of Yobito from 2005/06 to 2015/16 for validation. The result for 15 January and 15 February is shown in Fig. 7. The figure shows that although the absolute values of calculated hi are smaller by 10–20 cm than the observed thicknesses, they are reasonably well correlated, with correlation coefficients (r) of 0.64 (p = 0.04) for 15 February and 0.52 (p = 0.10) for 15 January. When taking into account the difference in water depth between site 3 (16 m) and the nearshore region (~2 m), we consider that this result supports the validity of the calculated hi in our model to some extent. Therefore, we examine the role of snow in the thickening processes likewise for this period, based on the model results.
The calculated hi on 15 February ranged between 33 cm and 54 cm with an average of 42 cm, and showed a good (and statistically-significant) correlation with hsi, with r = 0.78 (p < 0.01). On the other hand, hi had almost no correlation with hci (r = −0.03, p = 0.92) and hsi showed a significant negative correlation with hci with r = −0.65 (p < 0.01). These results indicate that the SI formation certainly works to moderate the variation of hci caused by different meteorological conditions and essentially control the total thickness (hi). To elucidate the effect of snow on hi more, the correlations of hi on 15 February with mean air temperature and snow accumulation (both from the freeze-up date to 15 February) were examined (Fig. 8). This figure shows that hi is much better correlated with snow accumulation (r = 0.54, p = 0.03) than with mean air temperature (r = −0.37, p = 0.16). These results support those obtained in Section 4.1 and emphasize the importance of the role of snow in determining hi.
In this section, we examine aspects of the ice phenology and their long-term trends at Lake Abashiri with our model, and compare with other lakes in the northern hemisphere (Benson et al., 2012), in relevance to the climatological variability of air temperature and snow depth. To do so, we extended the calculation of hi for the period 1961/62 to 1999/2000 (39 years) during which time the meteorological data-sets at Abashiri are available. In this calculation, air temperature was used to determine the freeze-up date for each year as proxy data. The method to calculate hi, hci and hsi is the same as that described in the previous sections.
Here we describe how we determined the freeze-up date from the air temperature data. To develop the method, firstly we calculated ‘the mean air temperature on freeze-up date’ (Taf) determined with MODIS images for the 16-year period 2000–15 and observations. In estimating Taf, to avoid the temporary fluctuation of air temperature on a daily timescale, we obtained the regression of the temperature decrease from November to December as a function of time at Abashiri and then calculated the air temperature on the freeze-up date for individual years. By taking the average, Taf was estimated to be −1.44 °C. This method and criterion were applied to determine the freeze-up date for the extended 39 years. Yet, if the surface heat budget was calculated to be positive (surplus heat at the surface), the freeze-up date was delayed until the surface heat budget became negative continuously. This method is considered to be useful since the freeze-up dates determined in this method coincided with the results estimated from MODIS images within a difference of 4 days for 12 years (86%) among 14 years (for the remaining 2 years, judged from air temperature). By combining the results for these 39 years with those of the previous section, we could obtain the data-set of hi and break-up dates for the past 55 years (1961/62–2015/16).
The time series of the break-up dates calculated for the 55-year period were shown with the freeze-up dates in Fig. 9. In the figure, the break-up dates estimated from MODIS images for the most recent 16 years are also plotted with error bars. It is found that overall the dates predicted with the model are consistent with the MODIS observation, which allows us to examine the climatology and long-term trend in break-up dates. According to Fig. 9, the break-up date is estimated to be 17 April on average and shows a significant negative trend of −1.71 d/decade (p < 0.01) for 55 years, i.e. the break-up date has become earlier by 9 days for this half century. This trend is comparable with the mean trend of −1.9 d/decade for 75 lakes in the Northern Hemisphere for the 30-year period 1975/76–2004/05 reported by Benson et al. (2012). On the other hand, the freeze-up date is estimated to be 12 December on average and did not show a significant trend for the 55 years. The estimated slope of regression line, +0.74 d/decade (not significant), is smaller than the mean value of 1.6 d/decade reported by Benson et al. (2012). The ice cover duration is estimated to be 127 days and shows a significant negative trend of −2.4 d/decade (p < 0.01). This trend is smaller than the mean value of −4.3 d/decade reported by Benson et al. (2012).
As for calculated ice thickness on 15 February, the averages of hi, hsi and hci are estimated to be 43, 16 and 27 cm, respectively. It is notable that while mean hi is comparable with that observed in Finland (44 cm; Leppäranta, 2015), the fraction of hsi (37%) is larger than 30% in Finland. A significant decreasing trend (−1.4 cm/decade, p = 0.01) can be seen in hi for these 55 years (Fig. 10). The correlations of hi with hsi and hci are shown to be both positive unlike the results of Section 4.1 (Fig. 11). However, it is noticed that hi is much better correlated with hsi than hci, and hsi is weakly but negatively correlated with hci. This shows that hsi certainly worked to moderate the variation of hi caused by the difference in hci. Therefore, all these results basically support the concept of the hypothesis about the role of snow presented by Ohata et al. (2016).
As described in Section 3, no slush layer is included in this model, based on the assumption that the slush layer, caused by the flooding of lake water into snow on ice, freezes to form SI on a timescale of a few days. As elucidated by the January 2016 event, this assumption does not necessarily hold true when a big flooding event occurs. To reproduce hsi and hi well in the model, we will need to consider the freezing process of slush layer more carefully. However, to the authors’ knowledge, no direct measurements have yet been made on the freezing of the slush layer for lake ice. Therefore, in this section, we examine this process using our observational data obtained in the 2015/16 ice season, focusing on the timescale of SI formation.
According to in-situ observations, a 22-cm-thick slush layer just above the ice on the day of instrument deployment (29 January 2016) had been completely transformed to SI at the observation on 19 February 2016. Fortunately, we could monitor the evolution of the vertical temperature profile within the snow, slush and ice layers using the string of 12 thermistors during this period (Fig. 2). Due to the low thermal conductivity and large specific heat of slush, the amplitude of the diurnal variation, forced by the surface temperature, would be significantly reduced in the snow, slush and ice layers. Since SI conversion usually begins at the top of the slush layer due to the cooling by cold air (Leppäranta, 2015), it is expected that the thickness growth of SI can be detected by the decrease in temperature and the enhancement of temperature variation in the snow, slush and ice layers. In this way, we attempted to analyse the growth rate of SI as a case study.
Time series of temperatures monitored at 12 depths at site 3 are shown in Fig. 12. It is seen that the air temperatures showed a large diurnal variation, ranging from −25 °C to +10 °C or more. Except the temperature variation at T2 (Fig. 2) was significantly reduced due to accumulated drifting snow from 21 February to 7 March, the air temperatures at T2 and T1 varied almost in-phase during the period. Corresponding to this variation, the temperatures within dry snow at T3 (25 cm above the slush/ice boundary) show a relatively large variation, ranging from −5 °C to 0 °C (Fig. 12b). Note that the temperature at T4 (20 cm) in slush layer remained at nearly 0 °C at the beginning, started to decrease and vary with an amplitude of about 1 °C about one day after set-up. This indicates that SI formed up to the depth of T4. As for T5 (15 cm), the temperature was kept at nearly 0 °C until 1 February and then started to decrease and vary with an amplitude of less than 1 °C, indicating that SI formed down to the depth of T5 on 1 February. On the other hand, the temperatures at T6 (10 cm) and T7 (5 cm) were kept at 0 °C during the period except for several very small variations at T6 between 19 and 21 February. This is probably because the temperature variation attenuated significantly at such depths due to the insulating effects of snow and snow ice, and suggests the limitation of this method to detect SI formation. The temperatures at T8 (0 cm) at the interface between ice and slush, at T9 (−10 cm), T10 (−15 cm) in the ice and at T11 (−25 cm), T12 (−35 cm) in the water were kept at 0 °C during the observation, as expected.
To see it more clearly, the temperature gradients between the sensors in snow, slush and ice layers are shown in Fig. 13. It is found that about one day after set-up, the temperature gradient (T3–T4) started to decrease and instead the gradient (T4–T5) started to increase, indicating the formation of SI at T4 about one day after set-up. On the other hand, the gradient (T5–T6) was kept at 0 °C m−1 until 12:00 pm on 1 February, indicating that slush layer was maintained between T5 and T6 until this time. These results coincide with the above speculation based on Fig. 12. Besides, note that the temperature gradients between T4–T5 and T5–T6 varied in phase with equivalent amplitude from 12:00 pm on 7 February onward. This indicates that transformation from slush to SI between T6 and T5 was completed by this time.
As for the remaining 10-cm-thick slush layer between T6 and T8, in the light of the fact that the entire slush layer was transformed to SI by 19 February and SI formed up to the depth of T6 until 12:00 pm on 7 February, we consider that it was totally frozen for the preceding 12 days (7–19 February). On the other hand, the upward growth of SI should be negligible during the period, judging from the fact that hsi formed during this period (22 cm) almost coincided with the value (21 cm) estimated from the ice block collected in the vicinity on 19 February.
All these results are summarized in Fig. 14, which shows how the slush layer was transformed to SI during the observation period. It is found from this figure that the growth rate of SI in the slush layer is estimated to almost constantly 1.0 cm/d. Since the mean air temperature (−6.5 °C) and snow depth on ice (10–27 cm) during the observation period (29 Jan–19 Feb) were typical at the time of this season, here, the estimated SI growth rate might be assumed as representative at Abashiri. This value of 1.0 cm/d is also comparable with a value of 0.6 cm/d estimated similarly from the measurement of the vertical temperature profile in snow at the saturated basal snow layer on Antarctic sea ice with a surface air temperature of −11.6 °C and a snow depth of 14 cm (Toyota et al., 2011).
Finally, based on these results, we check the validity of the assumption used in the model, that slush is transformed to SI on the timescale of a few days. According to Fig. 4, the maximum accumulation rate of snow during 2015/16 is 20 cm/5 d (4 cm/d). Taking the compression rate of β = 2 during the formation of SI in the model into account, the SI growth rate is estimated to be mostly less than 10 cm/5 d. This means that in most cases the slush layer becomes frozen within 2 days. Therefore, we consider that the assumption used in our model basically holds true. The exact reason why such a thick slush layer existed on 29 January in 2016 is not clear, but we suppose that some ice crack might have occurred to induce a big flooding event before the observation.
To look at the relationship with ice phenology at Lake Abashiri, the interannual variability of air temperature and snow depth averaged for the period from the freeze-up date to 15 February for the 55 years are shown in Fig. 15. It is found that whereas the air temperature exhibited a significant increasing trend of +0.23 °C/decade, the snow depth did not have a significant trend for the period, although a slight difference in the variability of snow depth may be found before and after the mid-1990s, i.e. 26.3 ± 9.1 cm for 1961–2000 in contrast to 33.0 ± 14.3 cm for 2001–2016. The correlations of hi with mean air temperature and mean snow depth are shown in Fig. 16. This figure shows that hi is correlated negatively with air temperature (r = −0.49, p < 0.01) and positively with snow depth (r = 0.41, p < 0.01). This indicates that when we look at the long-term trend of ice phenology, while snow still plays an important role in determining hi in February, air temperature also affects it significantly. Keeping the interannual variations of meteorological data and the relationship with hi in mind, we conclude that the negative trend of the break-up date (becoming earlier) and hi (becoming thinner) may be attributed more to the warming trend of air temperature than the change in snow depth. This is possibly because air temperature tended to affect hi more due to thinner snow before 2000 compared with after 2001. Thus, it is found that the phenology of break-up date of lake ice at mid-latitude shows a negative trend for the past 55 years, similar to other regions at high latitudes, associated with the warming trend.
To examine the relationship between ice phenology and air temperature, time series for normalized anomalies of ice phenology and air temperatures are shown in Fig. 17. The shared variance between mean anomalies of freeze-up and the mean anomaly of the mean air temperature in November to December is 76%. Such a high value is not surprising because it comes mainly from our method to determine the freeze-up date with air temperature, as shown in Section 4.3. The shared variances for break-up date with the air temperature in March and ice cover duration with the mean air temperature in November to March also take high values of 35% and 47%, respectively. These high values explain how significantly ice phenology is affected by air temperature for the past 55 years. The high correlation between ice phenology parameters and air temperature becomes more prominent when taking eight-year running means. It is found in Fig. 17 that the interannual variability of the eight-year running mean of air temperature is well correlated positively with freeze-up date and negatively with break-up date and ice cover duration. These results are similar to those reported for lakes in Europe and North America, although the relationship seems to be stronger at lower latitudes (Weyhenmeyer et al., 2004; Livingstone et al., 2010). The reason why the trend of freeze-up date is less prominent than that of break-up date is considered to be firstly because the increasing trend of air temperature is twice as strong in March (0.43 °C/decade) as in November to December (0.21 °C/decade) and secondly because freeze-up date is affected not only by air temperature but also by lake-specific characteristics which determine the lake heat storage; for example, lake volume, as discussed in Vavrus et al. (1996), Williams et al. (2004), Duguay et al. (2006) and Benson et al., (2012).
This study aims to clarify the properties of lake ice at mid-latitudes, subject to relatively moderate air temperature, relatively heavy snow and abundant solar radiation even in winter, as contrasted with lakes at high latitudes. For this purpose, we conducted field observations at Lake Abashiri in Japan for three winters and developed a 1-D thermo-dynamical model for predicting ice growth at Lake Abashiri. Using this model forced by meteorological data-sets from Abashiri, we examined the role of snow in the ice thickening process and also estimated ice phenology parameters and associated interannual trends at this lake for the past 55-year period to compare with lakes at high latitudes.
Firstly, to confirm the observational results and verify the model, obtained from observations during the 2012/13 winter season (Ohata et al. 2016), follow-up campaigns were conducted throughout two additional winters (2014/15 and 2015/16). Based on the results obtained, the role of snow in ice thickening process was examined from numerical experiments with this model to verify the hypothesis that snow works to moderate the ice thickness variation caused by meteorological conditions, presented by Ohata et al. (2016). The results are summarized as follows:
Next, to further examine the validity of the model and the role of snow in the ice growth under various meteorological conditions, we extended the calculation of hi, hci and hsi for the past 16 years (2000/01–2015/16), during which MODIS images are available, to determine the freeze-up dates. The results were found to be consistent with those obtained from the field observations in the three winters, as shown below:
Furthermore, to estimate ice phenology parameters and their interannual trends at this lake for comparison with other lakes in the northern hemisphere, we extended the calculation for the past 55-year period (1961/62–2015/16), using our model. In the calculation, the freeze-up dates were inferred from the air temperature data at Abashiri, and the break-up dates were calculated with our model. The results are summarized as follows:
In this study, the properties of lake ice at Lake Abashiri were investigated as a case study for the lakes at mid-latitudes. As a result, it was shown that the main role of snow in ice growth is to moderate the variability of ice thickness caused by the change in meteorological conditions, and ice thickness is found to be more strongly affected by snow depth than by air temperature. Although the study is limited to a specific regional lake, the obtained results are expected to apply to other lakes where meteorological conditions are similar, irrespective of latitude. Given that air temperature is predicted to become higher along with enhanced snowfall at high latitudes in the future, the results obtained here might serve to understand the changes in lake ice properties there. In addition, in the light of the fact that SI formation is regarded as one of the most important growth processes of sea ice, especially in the Antarctic, the results are also expected to broader the understanding in this domain. In the case of sea ice, the SI formation processes appear to be more complex because the freezing of water which has infiltrated through brine channels is another SI formation process, in addition to the freezing of snow flooded by sea water (Maksym and Jeffries, 2000; Toyota et al., 2011). Since the former process can be neglected in the lake ice formation, the SI formation process of lake ice may be used to study the latter process in isolation.
Although many lake ice properties at mid-latitudes were revealed through this work, some details of the processes of snow ice formation and melting still remain unresolved. Further studies including field observation are required.
No potential conflict of interest was reported by the authors.
This study was carried out by the bilateral cooperative agreement between the Institute of Low Temperature Science (ILTS), Hokkaido University, and Abashiri City, and was financially supported partly by the Grant for Joint Research Program of ILTS. We sincerely thank Toshifumi Kawajiri of Nishi-Abashiri Fishermen’s Cooperative Association for his enormous support throughout the field observations, and the Fisheries Science Center of Abashiri City for their kind logistic support. We are very grateful to Takayuki Shiraiwa of ILTS for introducing us the field campaign at Lake Abashiri, Sumito Matoba of ILTS for kind arrangements regarding the measurement of δ18O and suggestions for our analysis and to Jun Nishioka of ILTS for financial support. We thank Eisho Kudo of the Okhotsk Ryu-hyo (drift ice) Museum for providing pictures. We also thank the Abashiri and Memanbetsu tourism associations for providing the records of ice thickness for ice fishing. We appreciate the help of Ryuta Kato and Jun-no Ishiyama of Hokkaido University with field observations. Comments from the anonymous reviewer were very helpful and to improve the manuscript.
Benson , B. J. , Magnuson , J. J. , Jensen , O. P. , Card , V. M. , Hodgkins , G. and et al. 2012 . Extreme events, trends, and variability in Northern Hemisphere lake-ice phenology (1855–2005) . Clim. Change 112 , 299 – 323 .
IPCC . 2013 . Climate Change 2013. The physical science basis . In: Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) (eds. T. F. Stocker , D. Qin , G. K. Plattner , M. M. B. Tignor , S. K. Allen and et al. ) Cambridge, UK , Cambridge University Press .
Kirillin , G. , Leppäranta , M. , Terzhevik , A. , Granin , N. , Bernhardt , J. and et al. 2012 . Physics of seasonally ice-covered lakes: a review . Aquat. Sci. 74 , 659 – 682 . DOI: https://doi.org/10.1007/s00027-012-0279-y .
Leppäranta , M. 2009 . Modelling the formation and decay of lake ice . In: The Impact of Climate Change on European Lakes (ed. D.G. George ) Aquatic Ecology Series , Vol. 4 , Springer , Berlin Heidelberg , pp. 63 – 83 .
Leppäranta , M. 2015 . Structure and properties of lake ice . Freezing of Lakes and the Evolution of their Ice Cover . Springer , Berlin Heidelberg , pp. 51 – 90 . DOI: https://doi.org/10.1007/978-3-642-29081-7_3 .
Livingstone , D. M. , Adrian , R. , Blenckner , T. , Geoge , D. G. and Weyhenmeyer , G. A. 2010 . Lake ice phenology . In: The Impact of Climate Change on European Lakes (ed. D. G. George ) Aquatic Ecology Series , Vol. 4 , Springer , Berlin, Heidelberg , pp. 51 – 61 . DOI: https://doi.org/10.1007/978-90-481-2945-4_4 .
Magnuson , J. J. , Robertson , D. M. , Benson , B. J. , Wynne , R. H. , Livingstone , D. M. and et al. 2000 . Historical trends in lake and river ice cover in the Northern Hemisphere . Science 289 , 1743 – 1746 . DOI: https://doi.org/10.1126/science.289.5485.1743 .
Ohata , Y. , Toyota , T. and Shiraiwa , T. 2016 . Lake ice formation processes and thickness evolution at Lake Abashiri, Hokkaido, Japan . J. Glaciol. 62 ( 233 ), 563 – 578 . DOI: https://doi.org/10.1017/jog.2016.57 .
Rouse , W. R ., Binyamin , J . , Blanken , P. D ., Bussieres , N . , Duguay , C. R . and co-authors . 2008 . The Influence of Lakes on the Regional Energy and Water Balance of the Central Mackenzie River Basin . In: Cold Region Atmospheric and Hydrogic Studies. The Mackenzie GEWEX Experience (ed. M. K. Woo), Vol 1. Springer , Berlin, Heidelberg , pp. 309–325 . DOI: 10.1007/978-3-540-73936-4_18 .
Simmons , A. J. , Willett , K. M. , Jones , P. D. , Thorne , P. W. and Dee , D. P. 2010 . Low-frequency variations in surface atmospheric humidity, temperature, and precipitation: Inferences from reanalyses and monthly gridded observational data sets . J. Geophys. Res. Atmos. 115 , D01110 . DOI: https://doi.org/10.1029/2009JD012442 .
Toyota , T. , Massom , R. , Tateyama , K. , Tamura , T. and Fraser , A. D. 2011 . Properties of snow overlying the sea ice off East Antarctica in late winter 2007 . Deep Sea Res. II 58 ( 9–10 ), 1137 – 1148 .
Vavrus , S. J. , Wynne , R. H. and Foley , J. A. 1996 . Measuring the sensitivity of southern Wisconsin lake ice to climate variations and lake depth using a numerical model . Limnol. Oceanogr. 41 ( 5 ), 822 – 831 .
Weyhenmeyer , G. A. , Meili , M. and Livingstone , D. M. 2004 . Nonlinear temperature response of lake ice breakup . Geophys. Res. Lett. 31 , L07203 . DOI: https://doi.org/10.1029/2004GL019530 .
Willett , K. M. , Williams , C. N. Jr , Dunn , R. J. H. , Thorne , P. W. , Bell , S. and et al. 2012 . HadlSDH: an updated land surface specific humidity product for climate monitoring . Clim. Past D 8 , 5133 – 5180 .