Impacts of Hydro-Climatically Varying Years on Ice Growth and Decay in a Subarctic River

Abstract

Predicting the future changes in river ice development and impacts on seasonal sediment transport requires more in-depth examination of present river ice cover growth processes. This paper therefore investigates: (1) the impacts of hydro-climatically varying years on river ice development in a Scandinavian subarctic meandering river and (2) the accuracy of existing analytical models for predicting ice thickness growth and ice decay. Stefan’s ice growth equation (version by Michel et al.) and Bilello’s ice decay equation are applied to varying hydro-climatic conditions experienced in the years 2013–2019. Estimates from these equations are compared with observed field conditions such as ice thicknesses, ice clearance dates and freeze-thaw days. Overall, the equations were most accurate in the winter of 2016–2017 when the maximum mid-winter snow thickness value was high, the number of freeze-thaw days was the closest to the long-term average of northern Scandinavia, and the rate of thermal snow-melt in the subsequent spring was slow. The equations would need to be adjusted to take into account expected future changes to conditions such as shorter winters, less snow formation and increased frequency of air temperatures crossing 0 °C.

1. Introduction

Estimating or predicting the thickness of river ice and timing of ice breakups in polar regions is made difficult due to the multitude of processes in ice formation, growth and breakup that are highly dependent on local river morphology and hydrology. River freezing tends to start with skim and border ice formation along the river banks and the edges of bars as soon as the water has cooled down [1]. The main parameters affecting the lateral growth of border ice are heat loss and flow velocity, with field studies showing a positive relationship between the cumulative number of days below 0 °C (Accumulated Freezing Degree-Days: AFDD) and the lateral extent of border ice. However, flow velocities also affect border ice development; a study done on the Nelson River in northern Manitoba showed flow velocities exceeding 1.2 m/s to decrease noticeably the lateral growth rate of border ice [1]. Frazil ice, slush or floe formation may also take place prior to stable ice cover formation in turbulent rivers [2]. Once a stable ice cover has formed, it thickens thermally downward due to heat transfer processes based on meteorological conditions.

Downward thermal thickening of stable ice can be relatively well described by the degree-day method using the Stefan equation [3,4]. However, this method is not suitable for estimating very thin ice i.e., <10 cm thicknesses [5]. Michel [6] defined coefficients to be used in the thermal ice growth equation by Stefan when the river is snow covered. Ashton [7] concluded that the ice thickening equation provides accurate predictions assuming that there is no formation of snow-ice. When a snow cover exists, snow-ice may form if the snow becomes saturated with water, either by river water overflowing the ice surface through cracks [8], or from rain during a freeze-thaw cycle. Thus, the freezing of snow and water mixture can result in upward thickening of ice, which is not described with the basic Stefan equation.

Dornan [9] has developed site specific statistical ice growth models for several hydrometric river stations and found that accumulated freezing degree-days (AFDD) had a significant influence on observed river ice thicknesses, whereas the cumulative snow thickness and density had a relatively lower influence. Significant literature exists where Stefan’s equation has been applied to match the simulated ice thicknesses with observations (e.g., [10,11,12]). Zhang et al. [13] applied the ice thickness calculation method by Ashton [14] and developed a simple model for snow wetness determination which was used to define the thermal conductivity of snow. Thus, the input parameters of their ice thickness calculation framework consisted of air temperature, snow thickness and the modeled thermal conductivity of snow in addition to constant values. The method used by Zhang et al. [13] resulted in better thickness estimates than the original Stefan [3] equation, which ignores the thermal resistance of snow on heat transfer between water and air. However, the equations of Zhang et al. [13] and Ashton [14] require a greater number of parameters for their calculation compared to the Stefan equation [6], which may hinder their application.

Two sources of the break-up of river ice are thermal break-up and mechanical break-up. The main driver of thermal break-up is meteorological conditions, with this break-up generally occurring where there is no strong spring runoff event [15]. A spring runoff event occurs due to rapid snow melt caused by rising temperatures and/or significant rainfall, with the resultant river flows resulting in a mechanical break-up. Temperatures are influenced by the solar radiation absorption of ice; thus, the albedo of snowpack and ice play a role in the ice decay process rate. While different ice decay equations have been developed, there lacks a universal method for estimating the ice cover break-up. However, hinge resistance at river banks could be applied to determine the onset of ice cover break-up in relatively straight rivers [16].

In addition, Bilello [17] used river ice decay data from Canada and Alaska and found reasonable site-specific relationships between incremental ice thinning and thawing degree-day parameters. Furthermore, Shen and Yapa [18] developed a modified method to estimate river ice growth and thinning by applying degree-day method. The capability of the degree-day method to estimate the timing of ice breakups has been analyzed, e.g., by [19,20], showing that the modeling of ice breakup dates and changes in ice thickness is possible in large rivers. Despite that previous studies analyzing for the changes in the river ice break-up dates within Europe, or fine-tuning of equations based on the data from North-American rivers, the previous studies have not fully analyzed how accurate the simple ice growth and decay equations (which are based on air temperature) work during hydrologically variable years in rivers of northern Europe. These rivers can be considered relatively small (e.g., few tens of meters wide) and are under the hydro-climatic influence of the Gulf stream [21].

The polar region is warming faster than any other region on Earth [22]. However, northern Europe is experiencing a greater sensitivity to warming when compared to areas of similar latitudes in North-America [23]. Klaviņš et al. [24] have already shown a clear decreasing trend in the duration of ice cover during the last 150 years in the River Daugava, which has become more intense during the last 30 years. Altogether, Klaviņš et al. [24] analyzed river ice break-up dates from 17 stations in the Baltic States and Belarus. The climatic variables (i.e., North-Atlantic Oscillation (NAO) and Arctic Oscillation (AO)) have been stated as the reason for variation in freeze-up and break-up dates in Nemunas River, Lithuania, during the last 150 years [25]. NAO and AO, which impact on thermal anomalies, have shown positive correlations with freeze-up dates during mid-November to early January period [25]. However, statistically significant negative correlations during the period of early January to late March have indicated earlier than normal break-up dates in positive AO/NAO phases [25]. Šarauskienė and Jurgelėinė [26] have studied data from 13 water measurement stations (within eight rivers) in Lithuania and found that there has been a negative trend in ice duration and positive trend in freeze-up dates. Records from the Northern Hemisphere from 1846 to 1995 indicate a reduction in the ice season by 12.1 days per 100 years; however, records extending to 2005 indicate a further reduction of the ice season to 19.5 day per 100 years [23].

The expected future changes in discharges/flooding, sediment transport and channel morphology [27,28,29,30] in areas within cold climate region, such as northern Europe, are ultimately controlled by these river-ice characteristics, which will be affected by future shortening of the frozen period, earlier ice break-up, diminishing seasonally frozen ground, increasing freeze-thaw cycles and increased frazil ice occurrences [12,22,31,32]. To be able to predict these future changes in river ice characteristics, and further their impacts on river hydro-morphodynamics and related flood and erosion hazards, analysis of thermal ice growth and decay equations based on observations of hydro-climatically varying years are needed from rivers of the polar region.

Therefore, the aims are to (1) investigate the impacts of hydro-climatically varying years on the river ice development in a subarctic meandering river and (2) parameterize and analyze the accuracy of realistic ice thickness growth and decay models to examine climate-induced changes in river ice. This includes analyses of whether the existing ice thickness growth models, and decay equations can sufficiently describe the reality at present. Ice growth and decay equations are applied for estimating the river ice development during 2013–2019 years of varying hydro-climatic conditions (i.e., from mild to severe winters). The established methods, namely, the Stefan’s ice growth equation (version by Michel et al. [6]) and the Bilello [17] decay equation, are now for the first time applied for detecting multi-annual ice growth and decay processes of a subarctic river, within the polar region of Northern Europe. The novelty of the approach is in the climatic aspect, in particular in showing how the freeze-thaw cycle impacts on the ice growth and ice decay processes. The performance of these equations has not been tested previously by using multi-annual freeze-thaw data from hydro-climatically varying years. Thus, these modeled results are for the first time compared against multi-annual field observations (e.g., ice thickness and snow depth) measured during the mid-winters, time-lapse camera pictures throughout autumn-winter-spring periods, and hydro-climatic observation data from the Finnish Meteorological Institute.

2. Study Area

The study area is located at the meandering Upper Pulmanki River (Northern Finland), which is a tributary of the Tana River (Figure 1). The Pulmanki River is within the polar region (circa 70° N latitude) and belongs to the cold climate category without a dry season and with a cold summer (Dfc category of [33]) (Figure 2 and Figure 3). The hydro-climatic conditions, particularly the generation of low pressure systems, are affected by the Gulf stream heating the North Atlantic Ocean [21]. The study site has large seasonal hydro-climatic variations, below 0 °C winter air temperatures, and annual maximum discharges during spring snowmelt events. Furthermore, the region is predicted to experience warmer winter seasons in the future [22,30]. The study site has the advantage of no upstream lakes, which has been previously found to be an influential parameter for the malfunction of the ice growth analyses [12]. The temperatures often reach below 0 °C in early October (Figure 2 and Figure 3), and begin rising back above zero in early April. Thus, the freezing period lasts approximately seven months a year. The mid-winter ice bottom roughness was classified visually as smooth–rough, following the definition of [34]. Similar ice roughness occurred each winter from 2013–2019, as seen from photos taken from drill holes. Summer discharges are typically around 4 m³/s, but during spring the discharges can be around 50 m³/s [35]. However, less than 2 m³/s discharges have been observed in autumn and winter [36]. Thus, this river represents a typical subarctic river with its hydro-climatological characteristics and ice cover processes.

3. Materials and Methods

3.1. Field Measurements

The ice cover thickness data was gathered on 27 February 2014, 4 February 2015, 16 February 2016, 15 February 2017, 8 February 2018 and 8 February 2019 from the symmetrical meander bend locating 0.5 km upstream (straight distance) from the river outlet (Figure 1). During the sampling dates, snow depths were measured manually and their locations were measured with a RTK-GNSS (Real-Time Kinematic—Global Navigation Satellite System). In February 2014, we made snow thickness estimations based on pictures taken during the field survey, as no snow measurements were taken since the river ice was mainly snow-free and drift snow was only along some shorelines, as indicated in [11]. Most of the river ice was without snow cover in 2014, but the depths next to the banks were approximated against measurement equipment seen in the pictures. We installed two time-lapse cameras by the river to film the river ice development. These cameras located approximately 1.5 km upstream (straight distance) from the studied bend to minimize the human disturbance (e.g., by hikers) to the cameras (Figure 1, bank location). Continuous data (every two hours) were collected since August 2015, as either one or both cameras worked well from that point onwards. Due to battery problems and malfunction of the cameras, we did not have data from 2013 autumn–2015 spring. In addition, we installed on 7 February 2019 one camera close to, i.e., approximately 400 m upstream from the ice thickness measurement site (Figure 1, the HOBO (i.e., Honest Observer By Onset) sensor location). This was done to compare how well the cameras at the bank location represent other areas of the river. All cameras took pictures every two hours. From the cameras, the timing of the freezing and thawing of the ice cover were identified and used to compare against the ice decay calculations. The cameras also had air temperature sensors, however, without radiation shields.

In addition, the study used the continuous air temperature, rain intensity and snow depth data from the automated weather station of Finnish Meteorological Institute (FMI). This station is located in Nuorgam village (Utsjoki municipality), Tana River valley, approximately 20 m above sea level and 20 km north from the study site (Figure 1). The altitude is similar to the Pulmanki River study site which is approximately 15 m above sea level. A separate air temperature sensor (HOBO microstation) was installed at the same location as one of the cameras on 8 February 2019 (Figure 1, the HOBO location). The HOBO temperature data covered the whole spring period until the end of May 2019. The temperature data sets were compared to check the applicability of the FMI’s continuous data for the decay calculations.

The correlation between the air temperatures of the FMI data and the Pulmanki River data (HOBO) was 0.964 (Figure 4). Therefore, we can state that the FMI data could be applied for estimating the ice growth and decay processes at the upper Pulmanki River. However, when making detailed comparisons of temperature differences, it shows that the diurnal temperature variation is greater in the Pulmanki River valley than at the FMI’s station, especially in spring, i.e., after the end of March. The HOBO station at the Pulmanki River had a radiation shield, and was not located straight at the sun. Due to these diurnal differences, we have conducted the 2019 ice decay calculations based on both FMI’s and the Pulmanki River HOBO station’s data. When the time-lapse cameras’ temperature values were compared to the FMI’s data (see Figures S1 and S2 from the Supplementary Materials), it was noticeable that the temperatures were higher than the FMI’s temperatures, for example on some dates during spring as the sun shined directly at the camera. Note that the accuracy of air temperatures measured with the time-lapse camera without a radiation shield is less than the accuracy of the FMI’s data. However, in all data sets, the temperature magnitudes and whether temperatures are below or above 0 °C are suitably aligned.

3.2. Applied Ice Growth and Decay Equations

The AFDD values and the accumulated thawing degree-days (ATDD) were calculated based on the data from FMI for the whole period of 2013–2019, and on the HOBO air temperature sensor for the spring 2019 period. AFDD is the sum of the daily average air temperatures on days less than 0 °C. The summing starts at the first day below freezing and continues to the last day of the freezing in spring [10,37,38,39]. We tested the calculations of the AFDD for a full season (until the last day of freezing) and until March, as it is not recommended to use it to estimate ice thickness after the peak AFDD [40].

Modeling of freshwater ice during the last decades have its origin in the link between air temperature and ice growth, which was first established by Stefan [3] on sea-ice growth. Michel [6] then elaborated on Stefan’s analysis based on field experiments in rivers and published recommended correction coefficients to be used in ice thickness estimation of snow covered rivers. The coefficients reduce the ice thickness growth rate and mimic the insulating effect of snow. Hence, ice thickness (m) was calculated based on the Stefan formula [10,38,41]:

(1)$\eta ={\alpha }_h\sqrt{S}$

where S is the same as AFDD and ${\alpha }_h$ is the empirical coefficient. ${\alpha }_h$ was estimated as 0.015 m/°C^−1/2 day^−1/2, based on the knowledge of the area and after Michel [6]. The river can be considered as “Average river with snow”, which according to Michel [6] can have coefficient values of 0.014–0.017 m/°C^−1/2 day^−1/2. The mid-winter maximum ice thickness is gained as a result of this equation. The empirical coefficient (always less than 1.0) is site specific. Most importantly, this coefficient accounts for local factors like river morphology and curvature, the temperature difference between air and surface and insulating effect of snow, which actually have an effect on ice growth [42].

The ice decay was calculated based on the empirical Bilello [17] equation:

(2)$h=h_{max}-A{S}_T$

where $h$ (cm) is the thickness of the decaying ice sheet, $h_{max}$ (cm) is the maximum ice thickness at the start of the decay, A is an empirical constant (cm/°C) (average slope value) after [17]. Bilello [17] suggested the use of accumulated thawing degree-days (ATDD = $S_T$ in Equation (2)) to describe the decay of ice cover from their maximum thickness, and for rivers the melting degree-day is calculated base of −5 °C as the majority of river ice covers exhibit 10% or more thinning before the temperature exceeds 0 °C. The Equation (2) indicates the timing of thermal river ice break-up, i.e., when there is no more ice in the river. We also applied the Equation (2) so that we used average ice thicknesses from the same sampling occasion instead of the maximum. That way, the equation revealed the beginning time of ice melt and cracking each spring. In any case, the calculations with the Equation (2) initiated from the field campaign date, when we measured the maximum ice thicknesses of the area.

The results based on different equations were compared to the field measurements and FMI’s observation data. Freeze-thaw days (n) were calculated based on the days when maximum temperature was above 0 °C and minimum temperature was below 0 °C (T_max > 0 °C, T_min < 0 °C) [43]. This was applied for the detection of the differences between years, and as an explanatory variable for detecting the goodness of fit of estimated ice thicknesses and decay calculations.

4. Results

4.1. The Annual Variation of the Ice Cover Development

The time-lapse cameras have captured the river ice development since August 2015 every two hours (Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 and Supplementary Materials: Videos S1-9). Each year (2015–2018) frazil ice was present prior to stable ice cover formation (Table 1). Frazil ice events were the most common in autumn 2018 with five events were observed. Annually, the first ice rafts (i.e., loose sheets of surface ice floating downstream) formed on 8 October 2015, 27 October 2016, 23 October 2017 and on 7 October 2018. The ice rafts started to pile on top of the mid-channel bars first in all autumns, except in 2018, and later border ice started to form at the shoreline (on 30 October 2015, 29 October 2016, 25 October 2017, 27 October 2018).

The river froze totally (i.e., full ice cover) for the first time by 14 December 2015, 28 December 2016, 23 November 2017 and 25 December 2018. However, each autumn, the first complete ice-cover melted and then reformed, with the outer bank side of the channel experiencing melting. In the winter of 2015–2016, there was only one melting period before a stable ice cover formed, which lasted from 27 December 2016 onwards. In the 2016–2017 and 2017–2018 winters, a complete ice cover was formed five times, and just after 31 January 2017 and 21 January 2018, the mid-winter conditions prevailed and the long-lasting stable ice cover formed. The 2018–2019 winter was relatively colder (Table 2), with only one melting period and the second complete ice cover appearing on 8 January 2019, lasting the whole mid-winter period.

The snow depths (Table 2, Figure 10) at the Pulmanki River varied between years, from almost zero (2014) to up to 62 cm (2016). Most of the river ice surface was visually estimated to be snow-free during the field survey in February 2014. The second thinnest snow layer was in February 2018, when the snow depths were measured to be less than 10–20 cm on the Pulmanki River ice. The FMI’s snow depth data (note that these were not measured at a river) showed that the 2019 winter had the thinnest snow-cover. The ice thicknesses were the greatest in 2014 winter when there was least snow on top of the ice (Table 1 and Figure 11). Whereas, the thinnest ice cover was measured in February 2019. There has been a significant decrease in ice thickness (Max ice thickness: R² = 0.725, p = 0.0314, Average ice thickness: R² = 0.897, p = 0.004; Figure 12). The average and minimum temperatures of February decreased at the same time, but the maximum temperatures of February were also the highest in the end of the observation period (Table 2).

During each spring, small openings in the river ice cover appeared and re-froze before the final spring melting started on 9 April 2016, 1 May 2017, 24 April 2018 and 14 April 2019 at the bank location (Table 3). Only at the HOBO location, the camera did not show any melting and re-freezing of small openings before the beginning of the final melt on 16.4.2019, two days later than at the bank location. The snow cover decayed the fastest in 2016 spring and 2018 spring (Figure 2 and Figure 3). The ice decay was observed to be the earliest in 2015, taking place in early April (Table 3, Figure 6). The ice decay was classified thermal during springs 2015–2018. In spring 2019, for the first time during the observation period 2013–2019, a mechanical ice break-up took place and was over within 2 h (Figure 8 and Figure 9). Mechanical ice break-ups are typically closely related to water surface elevation difference at the time of the freeze-up and break-up as well as the ice cover thickness prior to the start of melt [44]. In 2019 spring, the ice melting started thermally, but the final break-up was mechanical, which is reflected in the short time difference between the beginning of melting and the decay of ice (Table 3). The river ice disappeared at the flood peak, and the trees within the impact zone were peeled off from the bark and knocked over.

4.2. The Performance of the Ice Growth and Decay Equations

Overall, the ice thicknesses calculated with the Stefan equation corresponded well with the observations, but annual differences occurred (Table 2, Figure 11 and Figure 12). The Stefan equation (when using the full season AFDD) underestimated the maximum ice thicknesses by 8–16 cm in 2013–2014 and 2014–2015 winters compared to the measured maximum ice thicknesses (Table 2, Figure 11 and Figure 12). These years had 15 freeze-thaw days, the greatest measured ice-thicknesses (70 and 60 cm, respectively) and the maximum snow depths on the river ice during the field surveys were the thinnest (5 and 26 cm, respectively) (Table 1 and Table 2, Figure 4, and Supplementary Materials: Figures S1–2, Table S1 and Videos S1–9). Overestimation (by 7–9 cm) took place during the years when there was only two freeze-thaw days (2017–2018) or when there was the greatest number of the freeze thaw days (20 in 2018–2019) (Table 2, Figure S2 and Figure 6). Thus, overestimation occurred the greater the long-term average number of freeze-thaw days was, which was 13 for the whole period (2013–2019). Thus, the best correspondence (difference only from −1 to +2 cm) was during 2015–2016 and 2016–2017, when there had been seven and 16 freeze-thaw days, respectively, and when maximum snow thicknesses at the Pulmanki River study site were its greatest (62 and 54 cm, respectively).

The ice thickness estimates were calculated with the Stefan equation also based on a shorter period of AFDD data, i.e., until 1st of March, and compared to the observed average ice thicknesses. The greatest underestimation (9 cm difference) was again during the 2013–2014 winter, when there was hardly any snow on the ice. The greatest overestimation of ice thicknesses was in 2017–2018 and 2018–2019, 18 and 17 cm, respectively, which had the lowest and highest number of freeze-thaw days, respectively. The best correspondence (2 cm overestimation only) was in 2016–2017, when the number of freeze-thaw days (16) was the closest to the average of the whole 2013–2019 period and had the second thickest maximum snow depth measured at the Pulmanki River. The year 2016–2017, which had performed well regarding the maximum ice thickness measurement and calculation comparisons, did not show as good correspondence. It overestimated ice thickness by 7 cm, which could be explained by the low number of freeze-thaw days.

The comparison of the results of the Stefan equation based on the AFDD data until March and observed maximum ice thicknesses generally showed a slightly different pattern in comparisons. All years underestimated the ice thicknesses except the years 2017–2018 and 2018–2019, which had a correspondence within 1 cm. Thus, these correspondences were best of all comparisons. Of note is the year 2015–2016, which had only 2 cm of underestimation when compared to the maximum ice thickness measurement (Table 2). For the year 2016–2017, measured maximum ice thicknesses showed good correspondence with the calculations applying full season AFDD, and when the average ice thickness values were compared to the calculations using AFDD data until March, this year had the second smallest underestimation (8 cm). The least correspondence was again during the years 2013–2014 and 2014–2015 with 22 and 11 cm underestimation, respectively, which had the thinnest snow-cover and thickest maximum and average ice cover measured at the Pulmanki River. Thus, when considering all the comparisons, the ice thickness estimates were the best during 2015–2016 and 2016–2017, when the freeze-thaw days were close to the long-term average, and the maximum snow depths were the thickest.

When the maximum ice thickness was applied in the ice decay calculations, the ice decay equation (Equation (2)) was able to define the ice clearance date at the bank location perfectly in spring 2017 (Table 3: “decay calculations, max ice thickness”). The beginning of the ice decay was also perfectly defined in spring 2017, when applying average ice thickness in the calculations. The melting period was longer in 2017 (temperature rise was more gradual) than in the other years (Figure 5). Thus, the decay equation with maximum ice thickness and with average ice thickness worked the best, when the air temperature increase was less intensive. The 2016–2017 also had the lowest AFDD value.

Overestimations of the ice clearance dates at the bank location were in years, which had a rapid snow-melt (Figure 2 and Figure 3). Those were 6, 5 and 8 days, in springs 2016, 2018 and 2019, respectively (Table 3). In addition, the beginning of the ice melt was overestimated with 19 days in 2016 spring and underestimated with 13 days in 2018 spring. The overestimation in the beginning date of the ice decay was less in 2019 spring, being seven days. The winter 2017–2018 had the highest AFDD, and 2018–2019 had the second highest AFDD. Third highest value was in the 2015–2016 spring. These three years also had the number freeze-thaw days the furthest away from the longer-term average value. It was noted that the ice decay calculations with average ice thickness actually resulted close to the ice clearance dates during spring 2016, 2018 and 2019 (with one, three and four days, respectively). Thus, during those years the application of average ice thickness seemed to predict the ice clearance date and not the date for the beginning of the ice melt. In conclusion, the ice decay equation worked the poorest during years which had a rapid temperature increase in the spring.

When the HOBO location’s air temperature data was applied in spring 2019 for calculating the ice decay date, there was only three days overestimation when using maximum ice thickness value, and five days overestimation occurred in the case when calculating the beginning date for the ice melt with average ice thickness values. Thus, it seems that during the spring when the mechanical ice break-up occurred, the observed ice decay at the straight reach (i.e., the HOBO location) corresponded better to the calculation results than the observed ice decay at the bank location, which represents a meander bend.

Note that when the ice decay calculations based on the FMI’s data and the HOBO data were compared in spring 2019, there was only one day difference betweeen the calculated ice-clerance dates when these different temperature data sets were applied. This very small difference in results gives confidence that the FMI station’s data was indeed applicable to Pulmanki River’s ice decay calculations, despite the 20 km distance to the Pulmanki River study site.

5. Discussion

When summing up the findings, the overall best performance of both equations (Stefan and Bilello) was in 2016–2017, when the maximum snow thickness values were high, the number of freeze-thaw days was closest to the long-term average of freeze-thaw days, and when there was both slow snow-melt and air temperature increases in the spring. Lind et al. [12] have also found that the stability of the ice cover depends mainly on the number of freeze-thaw days and the length of the cold period. Thus, our findings are in line with their study, as the most unpredictable years were those with the parameters furthest away from the long-term average number of freeze-thaw days. In addition, the greater the AFDD values were, the less well the decay processes were possible to be defined.

The least correspondence between the estimated (Stefan [6]) and observed ice thicknesses (underestimation) was during the years which had the thinnest snow-cover and the thickest maximum and average ice cover at the Pulmanki River. The ice thickness development was not possible to observe continuosly from 2013–2019, but only from the mid-winter conditions. However, visual inspection of the ice type was made in winter 2013–2014 through the bore holes and a 10–20 cm thick layer of snow ice was detected on top of the thermally formed black ice [11]. It is likely that river water had inundated on top of the ice cover at some point during the winter since the amount of snow was low on the river ice during field work in February–March 2014, despite snow having been accumulated in the nearby forested areas. As the Stefan equation [6] does not account for snow-ice formation, the underestimation of the ice thickness would be at least partly explained by the snow-ice in the 2013–2014 winter. Contrary to the present study, Ma et al. [20], found overestimation of ice thickness calculations during the years with snow on the ice. They had applied Stefan’s equation without any correction coefficient. Their four observation sites located along the large Lena River, in Siberia. However, Ma et al. [20] do not state where the measurement locations of the maximum ice-thickness values were. As the cross-sectional variation of the ice-thickness can be expected to be greater in the wide and large Lena River compared to the 20 m wide Pulmanki River, the possible difficulties in measuring the representative maximum ice thickness could also explain the overestimation of the model applied by [20] in a large river. The empirical coefficient applied in the Stefan equation (from [6]), which was used in the present study, took into account the temperature difference between air and surface and insulating effect of snow. The next step further would be to test with the Ashton [7] equation, which takes into account the thermal resistance of the surface to air boundary layer [7,42,45]. However, these parameters are not available from as many study sites as the air temperature parameters for Michel equation [6], and therefore, the lack in data sets hinders the spatially wide applicability of the Ashton [7] equation. In any case, further studies in areas with longer and continuous snow pack thickness time series data are suggested. The application of ice thickness estimation methods, which require contiuous information on snow depths, such as by [13], would then also be possible.

The thermal ice growth equation is not expected to work in the future in the areas of polar region, which are similar to the cold climate category area [33] of Northern Fenno-Scandinavia, as there is expected to be less snow and a higher number of freeze-thaw days in the future [22,30,43], in addition to the increased frazil ice appearance [46,47]. Persson et al. [43] has stated that the average number of freeze-thaw days would increase from the average 10 to 30 by year the end of the 21st century in northern Scandinavia. The 2018–2019 winter had 20 freeze-thaw days at the Pulmanki River, which would be in line of the predictions that Persson et al. [43] had mentioned more than a decade ago. Furthermore, the maximum ice thickness values had a linearly decreasing trend during the observation period 2013–2019. The changes in air temperature and increased number of freeze-thaw days could be a possible explanation for this. Thus, overestimations of ice thickness estimates are expected in years when the number of freeze-thaw days differs from the long-term averages, and underestimations during years with diminishing snow depths.

In addition, adjustments to the ice decay equation and the applied parameter values would be needed for predicting ice decay processes in future hydro-climatic conditions, because earlier snow-melt floods [30] are expected, which also indicates an earlier and more intensive rise in air temperature than at present. Turcotte et al. [48] have stated that future earlier spring snowmelt periods influence the break-up process as the ice cover is less exposed to short-wave radiation (the sun angle increases over time) when the discharge begins to rise. Based on long-term data sets, changes towards earlier break-up dates have also been observed in Northern (polar/subarctic region) and Southern Finland (boreal region) [49,50,51]. Thus, an increasing potential for mechanical break-ups and higher ice-jam-induced water levels are expected [48]. As the ice decay equation of Bilello [17] did not explain the mechanical break-up timing at present, it is not either expected to predict those in the future. During springs, with rapid and intensive temperature rise, the results of this study indicated better performance of the equations when applying the average mid-winter ice thickness value. Therefore, we suggest that these average mid-winter ice-thickness values may also work in the future with changed hydro-climatic conditions (i.e., early and rapid melting). Thus, the applied ice thickness value (i.e., whether applying average or maximum value) in ice decay equation would require adjustment based on the observed hydro-climatic characteristics of each year. This need for analyzing the hydro-climatic characteristics of each spring season first, i.e., before calculating the ice decay with Bilello [17] equation, makes it difficult to predict the ice decay processes of the coming years.

Beltaos [52] has observed in Peace River, Canada, that despite a rise in winter air temperatures since the 1960s, solid-ice thickness has increased. They concluded that porous accumulation covers enhanced winter ice growth via accelerated freezing into the porous accumulation. They stated that together with the reduction in winter snowfall, this effect (i.e., freezing into the porous accumulation) could reverse the effect of warmer winters on ice thickness, thus explaining present conditions. We are not able to state, based on the Pulmanki River data, if this would be also observable in Pulmanki River in the long-term. However, it is important to note that the equations at present may not work in future changed conditions due to these changes in process magnitudes. Note also that the climatic region of the Pulmanki River is presently under the warming influence of the Gulf stream [21], which makes it different to the rivers of more continental climate areas, such as Peace River in Canada.

Lind et al. [12] had previously observed that there is great difference between inland and coastal river reaches (in Sweden) in how well the AFDD can be applied to predict river ice thicknesses. The inland reaches had a stable ice cover, but the coastal reaches had many melting periods during the winter, and differences in AFDD values. Thus, the performance of the Pulmanki River can be considered more similar to the studied inland reaches of Sweden regarding the AFDD values, even though the impact of the Arctic Ocean warms the area more than the areas further inland in Finland. However, the Pulmanki River resembled the coastal rivers of the study of Lind et al. [12] in the freeze-thaw day amounts (average of 15.6 freeze-thaw days). When compared to the global river ice classification model by Turcotte and Morse [53], it was found that some ice characteristics of the Pulmanki River differed from it. According to Turcotte and Morse, the riffle-pool meandering rivers which are in cold climates, may get the complete ice cover that remains in place all winter after the first freeze-up cold spell. In the Pulmanki River, despite being in subarctic region and a small meandering channel, many freeze-thaw periods occurred each year.

Recently automated image processing algorithms have been developed and used for detecting the timing of ice formation and break-up processes [2]. This method would be worth testing in different sites. In addition, Park et al. [54] have simulated with a land process model ice thicknesses, as well as changes in the river freeze-up and ice break-up dates, consistent with in-situ data and satellite observations over the pan-Arctic. Their study indicated that less insulating snow might introduce thickening ice cover in the region. Thus, the next step would be to analyses the temperature changes in other subarctic and arctic areas, and analyzing how simple equations based on thermal properties fit with the observed thermal melting and ice break-up timing (such as based on satellite or aerial data). If these equations also work in other locations during similar hydro-climatic years as we observed at the Pulmanki River, then future predictions of changes in ice melting would be possible. Hence, the timing of the spring snow-melt flood and its changes could also be predicted, and further studies of the impacts of river ice changes on hydro- and morphodynamics could be performed for varying hydro-climatic conditions.

6. Conclusions

Based on the results of this study, it was possible to get a more detailed picture of the annual variation of ice-growth and decay processes within a sub-arctic river. It was also possible to define how well the simple ice growth and decay equations, which are based on air temperature, work during hydrologically variable years.

The mid-winter maximum ice thicknesses decreased during the observation period, along with the decreasing average and minimum mid-winter (February, i.e., the field measurement period) temperatures. However, the maximum mid-winter temperatures had also increased at the same time, indicating that extreme air temperatures had increased. The ice clearance date varied greatly between years, with the last one in 2017. Each year, the freezing process included frazil and border ice formation prior to the full ice-cover development. The ice clearance was thermal in all but one year, when it was mechanical.

Overall, the best performance of the ice growth and decay equations was in 2016–2017 winter, when the maximum mid-winter snow thickness value was high, the number of freeze-thaw days was the closest to the long-term average of northern Scandinavia, and the rate of thermal snow-melt was slow during spring. Overestimations of ice thickness estimates are expected in those years, when the number of freeze-thaw days differs from the long-term averages, and underestimations are expected during years with diminishing snow depths. Thus, the thermal ice growth equation is not expected to work in future changed conditions as well as at present as less snow and a higher number of freeze-thaw days are expected. Due to an expected earlier snow-melt flood and earlier rise in spring-time air temperatures in the future, adjustments to the ice decay equation for its future application would be needed based on seasonal hydro-climatic conditions. It is suggested to apply the average mid-winter ice thickness value for calculating the ice clearance date when winter is short and a rapid spring temperature increase is experienced. Thus, future possibly changed hydro-climatic conditions will complicate the application of both ice growth and decay equations.

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4441/11/10/2058/s1, Figure S1: The differences in the time-lapse camera’s (bank location) air temperatures and the FMI’S observation station (Nuorgam) data from the key time step dates of river ice freezing and thawing during 2015–2018. During these years these are the only temperature references from the Pulmanki River area, in comparison to the FMI’s data. The time-lapse camera did not have sun radiation shield. Therefore, the temperatures were higher than the FMI’s temperatures for example on some dates during spring, as the sun shined directly at the camera. Also the time-lapse camera’s accuracy is less than the accuracy of the FMI’s data; Figure S2: The differences in the time-lapse camera’s air temperatures (both from bank and HOBO locations) and the FMI’S observation station (Nuorgam) and the own HOBO station’s (Pulmanki River: since 7.2.2019) data from the key time step dates of river ice freezing and thawing period in 2018–2019. The time-lapse camera did not have sun radiation shield. Therefore, the temperatures were higher than the FMI’s and HOBO station’s temperatures for example on some dates during spring, as the sun shined directly at the camera. Also the time-lapse camera’s accuracy is less than the accuracy of the FMI’s or HOBO station’s data, Table S1: Air temp data HOBO Pulmanki River and FMI Nuorgam 2013–2019, Video S1: Ice growth 2015 autumn Pulmanki bank.pm4, Video S2: Ice decay 2016 spring Pulmanki bank.mp4, Video S3: Ice growth 2016 autumn Pulmanki bank.mp4, Video S4: Ice decay 2017 spring Pulmanki bank.mp4, Video S5: Ice growth 2017 autumn Pulmanki bank.mp4, Video S6: Ice decay 2018 spring Pulmanki bank.mp4, Video S7: Ice growth 2018 autumn Pulmanki bank.mp4, Video S8: Ice decay 2019 spring Pulmanki bank.mp4, Video S9: Ice decay 2019 spring Pulmanki HOBO location.mp4). The ice thickness values, which had been applied for the equations, can be found from the Tables presented in the methods section. The hydro-climatic data used in our paper is mainly based on the Finnish Meteorological Institute’s temperature data sets from Nuorgam station, in addition to the air temperature measurements at Pulmanki River (i.e., HOBO station).

Author Contributions

Data curation, E.L. and M.K.; Formal analyses, E.L. and L.L.; Funding acquisition, E.L.; Investigation, E.L.; Methodology, E.L., L.L. and M.K.; Project administration, E.L.; Visualization, E.L. and L.L.; Writing—original draft, E.L., L.L.; Writing—review & editing, M.K.

Funding

The study was funded by the Academy of Finland (‘Extreme and annual fluvial processes in river dynamics—ExRIVER’ [grant number 267345]), by the Maj and Tor Nessling Foundation (‘Äärimmäisten ja vuotuisten fluviaaliprosessien vaikutukset jokidynamiikkaan’ [grant numbers 2013067 and 2015046], ‘The Influence of River Ice and Fluvial Processes on River Environments Now and in the Future’ [grant number 2016042]), and the Department of Geographical and Historical Studies (University of Eastern Finland). Funding was also received from the Emil Aaltonen foundation (kohdeapuraha: “Tenojokeen, Pulmankijokeen ja Kokemäenjokeen kohdistuvien tutkimusten maastomittauksista aiheutuviin matkakuluihin ajanjaksolla syksy 2016-kevät 2017” [grant number: 160154 KO]). Most recent funding came from the British Society for Geomorphology (research project title “Defrosting sedimentary systems: the impacts on the evolution and material transport of high-latitude rivers” [registered charity number: 1054260]).

Acknowledgments

The authors thank Tiia Tarsa, Marko Kärkkäinen and Santeri Ronkainen from the University of Eastern Finland and Claude Flener, Leena Laamanen, Jenni-Mari Vesakoski, Elina Kasvi, Jouni Salmela, Linnea Blåfield and Petteri Alho from the University of Turku, in addition to William Speirs (visiting from Australia), for field assistance. William Speirs is also thanked for checking the English language of the paper. The authors also thank the Kevo subarctic research station for their help with ice thickness measurements.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. (A) The Pulmanki River study site and the measurement locations. The cameras have been filming at the bank location since autumn 2015, and at the HOBO (i.e., Honest Observer By Onset) sensor location since February 2019. (B) The Upper Pulmanki River locates within the Tana River watershed. The detailed study location is marked with “A”. The Finnish Meteorological Institute’s (FMI) weather station is at Nuorgam, approximately 20 km north from the Pulmanki study site. (C) The Pulmanki River locates in Northern Finland and is a tributary to the Tana River that forms a border between Finland and Norway.

Figure 2. The air temperature, rain intensity and snow depth measured at the FMI’s Nuorgam weather station 2013–2016. The black dots depict the day of snow and ice thickness measurements at the Pulmanki River.

Figure 3. The air temperature, rain intensity and snow depth measured at the FMI’s Nuorgam weather station 2016–2019. Also shown is the winter-spring 2019, which was measured at the Pulmanki River (HOBO station), 400 m upstream from the ice-thickness measurement locations. The black dots depict the day of snow and ice thickness measurements at the Pulmanki River.

Figure 4. The temperatures compared in more detail between the FMI’s station at Nuorgam and the Pulmanki River’s HOBO station measurements. The diurnal temperature variation is greater in the Pulmanki River valley than at the FMI’s station.

Figure 5. The ice cover development during 2015 autumn. (A) Ice appears for the first time on 8 December 2015 (−15 Celsius degrees). The (frazil) ice is floating as rafts. (B) The ice starts piling on top of the mid-channel bar (24 December 2015). (C) This 8 November 2015 time step is the first time when the border ice has started to grow around the mid-channel. Before this, there had been four times frazil ice, ice piling on top of the mid-bar and melting. After the fourth appearance (5 November 2015), the ice did not anymore melt totally. (D) The thin ice is covering most of the channel on 18 November 2015. Only an open strip exists next to the outer bank. (E) There is more ice and snow on 10 December 2015. Some openings appear in the middle of the channel and next to the bank. (F) The river got full ice cover on 14 December 2015, but wet areas appeared around 20–25 December 2015. After 27 December 2015, there is permanent full ice cover, which lasts more or less until the melting.

Figure 6. The freezing and thawing during spring 2016 as an example of a thermal ice break-up. (A) There are small openings in the ice on 16 February 2016, but those still freeze during the end of the February. (B) Open strip appears next to the outer bank 6 March 2016. Between 24 February–6 March, there has been static closed ice cover and below 0 temperatures (of camera). (C) During the night preceding 17 March 2016, the temperatures have gone down, and the river has frozen again totally. (D) Water is on top of the ice, at both sides of the mid-channel bar, on 28 March 2016. After this time step, there is diurnal freezing and thawing, but the ice cover stays rather similar until 9 April 2016, when water clearly starts flowing on top of the ice. (E) The ice cover has clearly collapsed at both banks of the river on 13 April 2016. (F) 27 April 2016 the river is ice-free.

Figure 7. The ice cover development during 2018 autumn. (A) 12 October 2018 (at 8:15): There is now for the first time clearly frazil ice/ice rafts/ice clumps floating in the river. At 10:25:28, the amount even increased. (B) 23 October 2018 (at 14:15): Couple of ice rafts appear in the river. During the following night, it is −11 degrees. (C) 27 October 2018 (at 8:15): border ice starts to appear for the first time next to point bar and the mid-channel bar. This is the first year (of the whole period 2015–2019) where border ice develops and not just frazil ice/clumps of ice pile on the bars. (D) 19 November 2018 (at 10:15): for the third time, there is frazil ice/ice clumps/ice rafts floating on the river. The previous night had been around −8 degrees. (E) 24 November 2018 (at 10:15): Frazil ice/ice clumps/ice rafts are floating on the river again. The ice has started to pile up on the mid-channel bar. (F) 25 December 2018 (at 12:15): Full ice cover appears for the first time this autumn. The ice will slightly melt at the outer bank side after this time step, before the full ice cover comes back on 8 January 2019 and stays until the late winter/spring melting periods.

Figure 8. The thawing during spring 2019 as an example of a mechanical ice break-up. The camera has been situated at the HOBO location. (A) 11 February 2019: this time step is an example of a solid ice cover situation and cold mid-winter and snowy conditions. (B) 24 February 2019: water is flowing on top of the ice. However, there is solid ice cover below. Open-channel conditions must have been somewhere upstream. (C) 8 March 2019: before this time step, more snow has fallen and there has been a freezing period. This picture is an example of solid ice conditions during the cold period. This cold and static period lasted until 14.4. (D) 14 April 2019: water appears on top of the ice at true-left and true-right banks. Open-channel conditions must have been somewhere upstream. (E) 22 April 2019 at 15:37: the ice has risen up along with the water level. The rising started on 20 April 2019 at 15:37 time step. This is the highest the ice has risen before its break-up. There is still similar ice cover two hours later at 17:37. (F) 22 April 2019 at 19:37. The mechanical ice-break takes place. The river shows no more ice cover, as only small rafts are floating. Note that the trees/bushes are also lining, i.e., they are no longer standing straight up. The ice and water had peeled of the bark from the trees from c. 1 m height.

Figure 9. The thawing during spring 2019 as an example of a combined thermal and mechanical ice break-up. The camera was at the bank location. There had been similar mid-winter melting occasions in the 2016 spring. However, in this figure, only the melting and break-up times are presented to show the difference of this location to the HOBO location (Figure 8), where the processes were more mechanical than in this upstream location. (A) 14 April 2019: This is the first time step when there is a clearly continuous flowing water area at the outer bank. (B) 15 April 2019: This is the time step when both the water level and the ice cover have started to rise (difference to the previous 10:33 time step’s water level). Thermal decay of the ice cover is taking place. (C) 17 April 2019: a mechanical ice break-up has happened between this 12:33 and previous 10:33 time step. The water level has kept rising gradually between 15–17 April 2019. (D) 18 April 2019: there is no more ice in the river (no rafts, nor rim ice). The water level is still rising. (E) 20 April 2019: big ice rafts appear in the picture, as they have been transported form the upstream. (F) 21 April 2019: ice rafts are all gone, and the ice melt has ended. The water level kept rising after the picture was taken.

Figure 10. The snow depth variation in the Pulmanki study site during the mid-winter measurements. Least snow was present on the river ice during winters 2014 and 2018.

Figure 11. The ice cover thickness variation in the Pulmanki study site during the measurement periods. The mid-winter ice cover was the thinnest during February 2018 and February 2019.

Figure 12. The results of the Stefan equation (ice thickness, cm), freeze-thaw days (number) and AFDD (degree Celsius). The Stefan equation was calculated based on the whole winter season and until March. Average snow depth was calculated from December to February.

Table 1

The freezing characteristics and ice growth calculations with the selected equations. Average snow depth of FMI’s station is calculated for the three winter months (December–February). The overall average freeze-thaw day number throughout the whole observation period 2013–2019 was 13 days, when rounded.

Winter	Frazil Ice Occurrence Days in Autumn (n)	Average Snow Depth (FMI, cm)	Freeze-Thaw Days (n)	Max Ice (Pulmanki, cm)	AFDD (Degree Celsius, Whole Winter)	Stefan Equation (cm, Whole Winter)	Difference (cm, Whole Winter Calculation versus Measured Max Thickness)	Average Ice (Pulmanki cm)	AFDD (until March)	Stefan Equation (until March)	Difference (cm, March Calculation Versus Measured Max Ice Thickness)	Difference (cm March Calculation Versus Measured Average Thickness)
2013–2014	no data	no data	15	70	1274	53.6	−16.4	57.4	1041	48.4	−21.6	−9
2014–2015	no data	42.0	15	60	1221	52.4	−7.6	46.6	1109	49.9	−10.1	3.3
2015–2016	3	40.5	7	53	1350	55.1	2.1	43.5	1166	51.2	−1.8	7.7
2016–2017	2	52.4	16	54	1237	52.8	−1.2	44	931	45.8	−8.2	1.8
2017–2018	2	43.2	2	55	1840	64.3	9.3	38	1390	55.9	0.9	17.9
2018–2019	5	37.0	20	50	1450	57.1	7.1	33.9	1144	50.7	0.7	16.8

Table 2

The mid-winter (February) snow depths at Pulmanki study site based on field observations. The 2014 values were estimated based on pictures, as we did not have measurements from that mid-winter condition. Additionally, the min, max and average air temperatures (FMI’s data) of February of each year are presented.

Winter	Minimum Snow Depth (cm)	Average Snow Depth (cm)	Maximum Snow Depth (cm)	Minimum Air Temperature (°C)	Average Air Temperature (°C)	Maximum Air Temperature (°C)
2013–2014	0	3	5	−30.3	−4.3	4.8
2014–2015	10	21	26	−30.1	−6.8	2.1
2015–2016	22	33	62	−31	−8.1	2.7
2016–2017	10	26	54	−28.5	−8.4	3.9
2017–2018	4	17	30	−32.6	−14.4	2.7
2018–2019	15	27	43	−35.6	−14	5.1

Table 3

The observed (cameras at the bank location) and calculated ice decay/clearance dates based on Equation (2) and the measured maximum and average ice thicknesses during field surveys. The ice clearance date means that there was no more ice in the river on that date. The bolded dates denote when calculations done with the max and average ice thicknesses corresponded the best to the observed dates of ice clearance and beginning of ice melt, respectively. The dates in italics denote when the ice decay calculations done with the average ice thickness corresponded to the observed ice clearance date. Calculated dates within parentheses are based on the HOBO temperature data.

Calculated Ice Clearance Date (Decay Calculations with Max Ice Thickness)	Observed Ice Clearance Date (Cameras)	Difference (Days)	Beginning of Ice Melt/Cracking (Decay Calculations with Average Ice Thickness)	Observed Date of Beginning of Ice Melt/Cracking (Cameras)	Difference (Days)
23 May 2014	no data		20 May 2014	no data
8 May 2015	no data		2 May 2015	no data
3 May 2016	27 April 2016	6 days overestimation	28 April 2016	9 April 2016	19 days overestimation
8 May 2017	8 May 2017	exact	2 May 2017	1 May 2017	1 day overestimation
11 May 2018	6 May 2018	5 days overestimation	9 May 2018	24 April 2018	15 days overestimation
25 April 2019 (24 April 2019)	17 April 2019 at the bank location; 22 April 2019 at the HOBO location	8 days overestimation; 3 days overestimation	21 April 2019 (20 April 2019)	14 April 2019 at the bank location; 16 April 2019 at the HOBO location	7 days overestimation; 5 days overestimation