The role of snow in the thickening processes of lake ice at Lake Abashiri, Hokkaido, Japan

Abstract

To clarify the properties of lake ice at mid-latitudes subject to moderate air temperature, heavy snow and abundant solar radiation even in winter, we conducted field observations at Lake Abashiri in Japan for three winters and developed a one-dimensional (1-D) thermo-dynamical ice growth model. Using this model with meteorological data-sets, we examined the role of snow in the ice thickening process, as well as the Lake Abashiri ice phenology (including the interannual trend) for the past 55 years to compare with high latitude lakes. The ice was composed of two distinct layers: a snow ice (SI) layer and a congelation ice layer. The SI layer occupied a much greater fraction of total ice thickness than that at high latitude lakes. In-situ observations served to demonstrate the validity of the model. Freeze-up and break-up dates supplied by satellite imagery enabled further model validation prior to the availability of field data (2000/01–2015/16). Based on both observations and numerical experiments, it was found that one important role of snow is to moderate the variability of ice thickness caused by changes in meteorological conditions. Furthermore, ice thickness is more sensitive to snow depth than air temperature. When applied to an extended 55-year period (1961/62–2015/16) for which local meteorological observations are available, the mean dates of freeze-up and break-up, ice cover duration and ice thickness in February were estimated to be 12 December (no significant trend), 17 April (−1.7 d/decade), 127 d (−2.4 d/decade) and 43 cm (−1.4 cm/decade). For this long-term period, while snow still played an important role in ice growth, the surface air temperature warming trend was found to be a strong factor influencing ice growth, as reported for the high latitude lakes.

Keywords:

• lake
• ice thickness
• snow ice
• slush
• thermodynamic model
• freeze-up and break-up dates
• mid-latitudes

1. Introduction

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

2. Observations and data

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 km², 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.).

Fig. 1. Location map showing Lake Abashiri on (a) a wide scale and (b) a regional scale, showing three observation sites on the lake, two Japan Meteorological Agency (JMA) operational sites: Abashiri Meteorological Observatory and Memanbetsu Airport, the photograph spot at the top of Tento-san mountain and the near shore site for the past ice thickness record (Yobito). Modified from Ohata et al. (2016)

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

Fig. 2. A schematic picture of the installation of the 12 sensor thermistor string for measuring air, dry snow, slush, ice and water temperature at site 3. The vertical line represents the sensor cable with thermistor sensor depths referenced to the slush/ice interface on 29 January 2016.

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 T_a (air temperature), V_a (wind velocity) and H_s (snow depth on land) are available at both stations on a daily basis, Q_s (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, T_a, H_s and V_a calculated as the mean of the two meteorological station data, and Q_s, 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.

3. Summary of the thermodynamic ice model

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.

Table 1. Model parameters and constants. Modified from Ohata et al. (2016).

Parameter	Symbol	Value	Unit	Source
Congelation ice density	ρci	910	kg m^−3
Snow ice density	ρsi	870	kg m^−3
Snow density	ρs	300	kg m^−3
Air density	ρa	1.3	kg m^−3
Thermal conductivity of congelation ice	kci	2.07	W m^−1 K^−1	Yen (1981)
Thermal conductivity of snow ice	ksi	2.07	W m^−1 K^−1	Yen (1981)
Thermal conductivity of snow	ks	0.23	W m^−1 K^−1	Yen (1981)
Specific heat of air	Ca	1004	J kg^−1 K^−1
Latent heat of freezing	Lf	3.34 × 10^5	J kg^−1	Yen (1981)
Latent heat of sublimation	Lv	2.84 × 10^6	J kg^−1	Yen (1981)
Surface emissivity	ε	0.97		Omstedt (1990)
Transfer coefficient for latent heat	CE	137 × 10^−3		Andreas and Makshtas (1985)
Transfer coefficient for sensible heat	CH	1.37 × 10^−3		Andreas and Makshtas (1985)
Albedo, open water	α	0.07		Toyota and Wakatsuchi (2001)
Albedo, ice	α	0.3		Perovich (1998)
Albedo, snow	α	0.75		Pirazzini et al. (2006)
Transmittance, ice	I0	0.18		Grenfell and Maykut (1977)
Transmittance, snow	I0	0.0		Maykut and Untersteiner (1971)

4. Results

4.1. Field observations in the 2012/13, 2014/15 and 2015/16 ice seasons

4.1.1. Meteorological conditions

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.

Fig. 3. Time series of daily (a) air temperatures and (b) snow depths at AMO in 2012/13, 2014/15 and 2015/16. The vertical dashed lines in (a) indicate the freeze-up dates and the horizontal dashed lines in (b) indicate the snow depths on the freeze-up dates.

Table 2. Meteorological conditions from the freeze-up date to 20 February.

Year	Freeze-up date	Mean air temperature (^oC)	Mean snow depth (cm)	Increase of snow depth just after freeze-up (cm d^−1)
2012/13	11 December	−7.5	23.5	4.9
2014/15	10 December	−4.6	41.0	0.0
2015/16	23 December	−6.3	24.1	0.7

Note: The mean air temperature and mean snow depth (given by the snow thickness increment from freeze-up date) are the averages at the two meteorological stations; Abashiri Meteorological Observatory and Memanbetsu Airport. Increase of snow depth just after freeze up is the mean increase of snow depth for 2 days after freeze-up.

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.

4.1.2. Observational results

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⁻²), 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.

Table 3. Ice thickness and snow depth at sites 1, 3 and 6.

Year	Date	Thickness (cm)	Site
			No.1	No.3	No.6
2013	18–19 February	Total ice thickness	43	42	33
		Whitish ice layer thickness	12.5 (29)	19 (45)	24 (73)
		Translucent ice layer thickness	30.5	23	9
		Snow depth	28	26	23
2015	24 February	Total ice thickness	48	55	53
		Whitish ice layer thickness	19 (40)	31 (56)	32 (60)
		Translucent ice layer thickness	29	24	21
		Snow depth	26	17	29
2016	19 February	Total ice thickness	40	42	48
		Whitish ice layer thickness	11 (28)	19 (45)	31 (65)
		Translucent ice layer thickness	29	23	17
		Snow depth	14	20	11

Note: The numbers in parentheses denote relative portion of whitish ice layer to total ice thickness in ％.

Ice structures were commonly composed of whitish (upper) and translucent (lower) layers, and it was confirmed from the δ¹⁸O 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⁻³, 821 ± 81 kg m⁻³ and 380 ± 31 kg m⁻³ from the observation on 24 February 2015, and 872 ± 16 kg m⁻³, 838 ± 13 kg m⁻³ and 243 ± 37 kg m⁻³ 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.

4.1.3. Model results

To examine the validity of the model developed for the 2012/13 season (Ohata et al., 2016), the time series of h_i, h_si, h_ci and h_s 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 h_ci, h_si and break-up dates show a good agreement with the observed values except for h_si in January 2016 and h_i in March in both years. The discrepancies between model and observational results were 4 cm for h_i, 3 cm for h_si, 1 cm for h_ci and 12 cm for h_s 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.

Fig. 4. Model results of thickness evolution of total ice, SI and snow with observations shown by triangles (snow depth), squares (SI) and solid circles (total ice thickness) in 2012/13 (modified from Ohata et al., 2016), 2014/15 and 2015/16. The light blue curve indicates the snow depth on land obtained from the daily meteorological snow depth data by taking a 5-d running mean (the day and 4 d ahead). The purple, red and blue curves show the snow depth on the lake, SI thickness and total ice thickness predicted by the model, respectively. Freeze-up dates used in the model and break-up dates calculated by the model are indicated, along with those estimated with MODIS images (in parentheses).

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 h_i in March, in the light of the fact that it is attributed mainly to the underestimation of h_si (Fig. 4) and air temperature sometimes exceeded 0 °C in March (Fig. 3), the omission of the increase in h_i 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⁻² for 2014/15 and 2015/16 (the mean net surface heat flux from 15 March to the break-up date; 78 W m⁻² for 2015, 71 W m⁻² 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⁻²). The break-up dates estimated for 2015 and 2016 were the same dates (8 April) and that for 2014 was later (11 April).

4.1.4. Numerical experiments using the model

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.

Fig. 5. Observations of ice thickness (at site 3 on 19 February 2013, 24 February 2015 and 19 February 2016) and the ice thicknesses on 20 February 2013, 2015 and 2016, predicted by the model for three cases; (1) snow accumulates with SI formation (control run), (2) snow accumulates without SI formation and (3) no snow accumulates.

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 h_ci 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), h_ci 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 h_s 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 h_ci still remains as case 2, the difference in h_i is reduced to only 3 cm (= 45 − 42) by the increase in h_si and the difference in h_i found for case 3 (no snow) was reduced by about 75% in the presence of snow. This indicates that the difference in h_i 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 h_i on 20 February is 13 (= 68 − 55) cm. For case 2, h_i 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, h_ci 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, h_i became rather thicker in 2014/2015 than in 2012/13, although the difference in h_i 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 h_i 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, h_ci 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 h_ci calculated for case 1, where h_ci is thinner by 7 cm in 2014/15 than in 2015/16. However, h_i 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 h_i 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 h_i caused by differences in air temperature, SI formation plays a more essential role in determining h_i. In addition, comparison of the results between 2012/13 and 2015/16 suggests that although the timing of snowfall events affects h_ci due to the thermal insulation effect at the initial stage, h_i in February depends more on the SI growth. Schematic pictures which explain how snow works to moderate the variability of h_i 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).

Fig. 6. Schematic pictures showing the effect of snow and snow ice formation on the lake ice thickness for (a) the same snow amount and different air temperatures (cited from Ohata et al., 2016), and (b) different snow amounts and the same air temperature.

4.2. Validation for the past 16 years with MODIS images

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 h_i, h_ci, and h_si 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 h_i 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 h_i 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 h_i 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.

Fig. 7. The correlation coefficient between modelled and observed lake ice thickness at Yobito on 15 January and 15 February for 2005/06–2015/16. Diamonds and squares show ice thicknesses on 15 January and February, respectively. Correlation coefficients and p-values are shown by r and p in the figure. Yobito data are from the Abashiri tourism association.

The calculated h_i 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 h_si, with r = 0.78 (p < 0.01). On the other hand, h_i had almost no correlation with h_ci (r = −0.03, p = 0.92) and h_si showed a significant negative correlation with h_ci with r = −0.65 (p < 0.01). These results indicate that the SI formation certainly works to moderate the variation of h_ci caused by different meteorological conditions and essentially control the total thickness (h_i). To elucidate the effect of snow on h_i more, the correlations of h_i 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 h_i 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 h_i.

Fig. 8. The relationship between total ice thickness on 15 February and mean air temperature or mean snow depth from freeze-up date to 15 February, recorded at the Abashiri Meteorological Observatory for 2000/01–2015/16. The solid straight lines are the regression lines of best fit. Correlation coefficients and p-values are shown by r and p in the figure.

4.3. Estimation of the ice phenology for the past 55 years

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 h_i 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 h_i, h_ci and h_si 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’ (T_af) determined with MODIS images for the 16-year period 2000–15 and observations. In estimating T_af, 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, T_af 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 h_i 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).

Fig. 9. Freeze-up and break-up dates from 1961/62 to 2015/16, determined by the model and MODIS images. Error bars show the time ranges estimated for freeze-up and break-up dates from MODIS images. Diamonds and squares show freeze-up dates used for the model and break-up dates calculated by the model, respectively.

As for calculated ice thickness on 15 February, the averages of h_i, h_si and h_ci are estimated to be 43, 16 and 27 cm, respectively. It is notable that while mean h_i is comparable with that observed in Finland (44 cm; Leppäranta, 2015), the fraction of h_si (37%) is larger than 30% in Finland. A significant decreasing trend (−1.4 cm/decade, p = 0.01) can be seen in h_i for these 55 years (Fig. 10). The correlations of h_i with h_si and h_ci are shown to be both positive unlike the results of Section 4.1 (Fig. 11). However, it is noticed that h_i is much better correlated with h_si than h_ci, and h_si is weakly but negatively correlated with h_ci. This shows that h_si certainly worked to moderate the variation of h_i caused by the difference in h_ci. Therefore, all these results basically support the concept of the hypothesis about the role of snow presented by Ohata et al. (2016).

Fig. 10. The thicknesses of total ice, SI and CI on 15 February, determined by the model. The black solid line, dashed-dotted blue line and dashed red line show total ice thickness, CI thickness and SI thickness predicted by the model, respectively. The black solid straight line is the regression line of best fit for total ice thickness. The trend, shared variance (r²), correlation coefficient (r) and p-value (p) are shown in the figure.

Fig. 11. The relationships between total ice thickness (h_i), SI thickness (h_si) and CI thickness (h_ci) on 15 February for 55 years: (a) h_i vs. h_si, (b) h_i vs. h_ci and (c) h_si vs. h_ci. Correlation coefficients and p-values are shown by r and p in the figure.

5. Discussion

5.1. Validation of the SI formation process in the model

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 h_si and h_i 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.

Fig. 12. Time series of temperatures of (a) air, (b) snow, ice and water at site 3 from 29 January to 11 March 2016.

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⁻¹ 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.

Fig. 13. Time series of temperature gradients within snow and ice at site 3 from 29 January to 11 March 2016.

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 h_si 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).

Fig. 14. A schematic picture of snow ice growth at site 3. The vertical axis represents the heights with sensor names referenced to the ice/slush interface on 29 January. The horizontal axis represents time. The red line represents the boundaries between SI and slush or dry snow. The green line represents the snow surface. The red dashed line represents the supposed boundary between SI and dry snow which was not observed. The black dashed vertical and horizontal lines represent the dates and the corresponding heights of slush, respectively.

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.

5.2. Relationship between ice phenology and meteorological conditions

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 h_i with mean air temperature and mean snow depth are shown in Fig. 16. This figure shows that h_i 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 h_i in February, air temperature also affects it significantly. Keeping the interannual variations of meteorological data and the relationship with h_i in mind, we conclude that the negative trend of the break-up date (becoming earlier) and h_i (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 h_i 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.

Fig. 15. Meteorological conditions after freeze-up date to 15 February for (a) the mean air temperature and (b) the mean snow depth. The black solid straight line is the regression line of best fit. The trends, shared variances (r²) and p-values are shown in the figure.

Fig. 16. The relationship between total ice thickness on 15 February estimated by the model, and the mean air temperature or the mean snow depth from freeze-up date to 15 February for 1961/62–2015/16, recorded at the Abashiri Meteorological Observatory. The solid straight lines are the regression lines of best fit. Correlation coefficients and p-values are shown by r and p in the figure.

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

Fig. 17. Comparisons of anomaly time series for ice phenology parameters and air temperature from 1961/62 to 2015/16. All valuables are expressed as normalized anomalies (standard deviations of the measures being compared): (a) freeze-up date and the average of November and December temperatures, (b) break-up date and March temperatures and (c) ice cover duration and the average of November to March temperatures. Annual values are indicated by thin lines and eight-year running means by heavy lines; ice phenology parameters are blue and air temperatures are red. r² is presented for the linear relationship between the ice phenology parameters and the corresponding air temperatures.

6. Conclusions

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:

(1)	Commonly, lake ice was composed of whitish snow ice (SI) and translucent congelation ice (CI), and 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.
(2)	The model, which was developed based on the observational results in 2012/13, successfully reproduced the observed thicknesses of both hci and hsi for 2014/15 and 2015/16 as well, showing the validity of the model under different meteorological conditions. It was also found that the model can reproduce the break-up date well, with an accuracy of a few days.
(3)	The role of snow in lake ice growth was examined from numerical experiments with this model by comparing hi, hci and hsi calculated for three cases: (i) snow accumulates with SI formation (control run); (ii) snow accumulates without SI formation; and (iii) no snow accumulates. The results obtained for these three winters support the hypothesis that snow works to moderate the ice thickness variation caused by meteorological conditions, presented by Ohata et al. (2016), and it was shown that ice thickness is controlled more by snow depth than by air temperature.

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 h_i, h_ci and h_si 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:

(4)	hi calculated in the model are consistent with the ice thickness records measured at the near shore region of the lake.
(5)	There is a significant negative correlation between hsi and hci, indicating that snow works to moderate the change in hci through the formation of SI. The importance of snow was also endorsed by the fact that hi has a higher correlation with snow depth than with air temperature at Abashiri.

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:

(6)

The average freeze-up date, break-up date and ice cover duration were estimated as 12 December, 17 April and 127 days, respectively. On average, hsi (16 cm) accounts for 37 % of hi (43 cm), higher than about 30 % in Finland (Leppäranta, 2015).

(7)

Statistical analysis showed that the break-up dates have experienced a significant negative trend of −1.71 d/decade for the past 55 years, i.e. the break-up date has become earlier by 9 days for the recent half century. This trend is comparable with the trend for other lakes in the Northern Hemisphere for the 30-year period 1975/76–2004/05 reported by Benson et al. (2012). This is presumed to be because air temperature affected hi more strongly due to less snow before 2000 compared with that after 2001. Time series for normalized anomalies of ice phenology parameters and air temperature showed that they are highly correlated, which is consistent with the result reported for the lakes in Europe and North America. On the other hand, the freeze-up dates did not show a significant trend for the 55-year period. As for hi on 15 February, a significant decreasing trend (−1.4 cm/decade) was found, along with a much higher correlation with hsi than with hci. The observed negative correlation between hsi and hci supports our hypothesis about the role of snow.

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.

Disclosure statement

No potential conflict of interest was reported by the authors.

Acknowledgement

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 δ¹⁸O 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.