Estimation of the Dependence of Ice Phenomena Trends on Air and Water Temperature in River

: The identiﬁcation of changes in the ice phenomena (IP) in rivers is a signiﬁcant element of analyses of hydrological regime features, of the risk of occurrence of ice jam ﬂoods, and of the ecological e ﬀ ects of river icing (RI). The research here conducted aimed to estimate the temporal and spatial changes in the IP in a lowland river in the temperate climate (the Note´c River, Poland, Central Europe), depending on air temperature (TA) and water temperature (TW) during the multi-annual period of 1987–2013. Analyses were performed of IP change trends in three RI phases: freezing, when there appears stranded ice (SI), frazil ice (FI), or stranded ice with frazil ice (SI–FI); the phase of stable ice cover (IC) and ﬂoating ice (FoI); and the phase of stranded ice with ﬂoating ice (SI–FoI), frazil ice with ﬂoating ice (FI–FoI), and ice jams (IJs). Estimation of changes in IP in connection with TA and TW made use of the regression model for count data with a negative binomial distribution and of the zero-inﬂated negative binomial model. The analysis of the multi-annual change tendency of TA and TW utilized a non-parametric Mann–Kendall test for detecting monotonic trends with Yue–Pilon correction (MK–YP). Between two and seven types of IP were registered at individual water gauges, while di ﬀ erences were simultaneously demonstrated in their change trends over the researched period. The use of the Vuong test conﬁrmed the greater e ﬀ ectiveness of estimates for the zero-inﬂated model than for the temporal trend model, thanks to which an increase in the probability of occurrence of the SI phenomenon in the immediate future was determined; this, together with FI, was found to be the most frequently occurring IP in rivers in the temperate climate. The models conﬁrmed that TA is the best estimator for the evaluation of trends of the occurrence of IC. It was shown that the predictive strength of models increases when thermal conditions are taken into consideration, but it is not always statistically signiﬁcant. In all probability, this points to the impact of local factors (changes in bed and valley morphology and anthropogenic pressure) that are active regardless of thermal conditions and modify the features of the thermal-ice regime of rivers at speciﬁc spatial locations. The results of research conﬁrm the e ﬀ ectiveness of compilating a few models for the estimation of the dependence of IP trends on air and water temperature in a river.


Introduction
In a cool and temperate climate, ice phenomena (IP) in rivers are one of the elements of their winter hydrological and thermal regime. Information concerning river icing (RI) has a practical dimension, among others in regard to the identification of impediments to the operation of hydrotechnical equipment regulating water flows and levels, as well as the minimization of threats connected with ice jams and ice jam floods. As such, we are gathering data concerning the impact of these processes on the ecological state of waters, as well as the course of biochemical processes in the aquatic environment and time series and spatial trends regarding IP phenology and their connections with climate pattern trends. In the case of large spatial patterns of atmospheric circulation, emphasis has been placed on determining and indicating those systems that control regional air temperature fluctuations [34][35][36][37][38][39], but the empirical dependencies between dates of occurrence of ice on rivers and TA indices may not be reliable under future climatic conditions. Certain RI effects are strongly controlled by processes occurring on the scale of the catchment area and the river bed (e.g., the dynamics and intensification of river ice breakage fronts), which must be taken into consideration in the analyses of local factors. In the future, the verification of the results of model studies and the validation of RI models will require obtaining information concerning the scope of change of hydroclimatic conditions and the impact of heat exchange on the atmosphere-bed or bed-subterranean water levels, which will enable the identification of step changes in the thermal and ice regime of rivers. According to the Arctic Monitoring and Assessment Programme (AMAP) [28], in order to obtain better results of ice regime forecasting, it is necessary to achieve progress in the integration of models that consider the compilation of future changes in landscape hydrology; the exchange of energy between water, ice, and air; flow hydraulics; and ice mechanics.
The objective of this research was to estimate temporal and spatial changes in RI in connection with thermal conditions. The research focused on IP change trends in the lowland Noteć River (Poland, Central Europe), these connected with three RI phases: freezing, the persistence of IC and floating ice (FoI) drift. The identification of time trends while taking TA and TW into consideration was intended to indicate significant estimators that correlate with both the number and the probability of occurrence of IP in a given year. Due to zero inflation in observation series, which means that IP may not appear during a year, it was assumed that the process governing the number of occurrences of IP is accompanied by a dichotomous process that determines whether a phenomenon has a chance of occurring in a specific year. Additionally, by taking monotonic TA and TW trends into consideration, it was planned to provide an answer to the question of whether changes in TA have contributed to an increase in TW in relation to reference periods, while at the same time leading to changes in river thermal and ice regimes. The identification of IP variability and trends in rivers is significant due to the economic utilization and ecological importance thereof.

Study Area
The analysis was conducted for the Noteć River (391.3 km long), a lowland European river (Central Europe) that is a tributary of the Warta River in the catchment area of the Oder River that flows into the Baltic Sea (Poland) (Figure 1). The Noteć River, which constitutes an important element of the shipping lane connecting the largest rivers in Poland (the Vistula and Oder Rivers), has considerable economic significance and significant natural value. The largest tributaries of the Noteć River are: The Gwda, Drawa, Noteć Mała, and Łobżonka rivers. The river catchment area is located in the following macroregions: The Toruń-Eberswald proglacial stream valley, through which the main bed of the river runs, and the South Pomeranian Lakelands, through which the right-bank tributaries of the Noteć River flow.
The landscape is dominated by a relatively continuous sandy-loamy cover of the Quaternary deposits, associated with the accumulation activity of the ice sheet and the fluvioglacial and river waters. It is represented by frontal moraine hills, which are dissected by a network of glacial channels and river valleys. The substrate of the catchment consists of ablation sands and boulder loams. In the near-surface zone and on the surface, there are mainly sand and gravel formations of the Baltic Glaciation and the Holocene with a thickness of up to 40 m. These are river sediments of the middle and high terraces, dune sands, and peats. Peats build floodplain terraces, especially the valleys of the Noteć River and depressions with no outlets. The characteristic feature of the Upper Noteć is the presence of numerous lakes along its reach and near the tributaries, which occupy as much as 4% of the catchment area.  [23,40], partially modified).
The catchment area is predominantly agricultural and forest. The arable land and partly grassland are typically agricultural lands in the eastern part and forests in the western part ( Figure 1). The study area is diverse in terms of natural conditions and is characterized by a mosaic of habitats with priority riparian forests and with a high diversity of hydrogenic habitats. The river valley, due to its significant natural values, is protected by its inclusion in the Natura 2000 areas, such as "Dolina Noteci" (the Noteć River Valley) and "DolinaŚrodkowej Noteci i Kanału Bydgoskiego" (the Middle Noteć River Valley and the Bydgoszcz Canal). In its middle reaches, the river is canalized and accessible by ships, while its water levels are regulated by a sluice system. A detailed physical and geographical description of the catchment area was presented in earlier studies authored by Graf [23,40].
The study area is the point of convergence of a number of climatic regions, i.e., the Western and Eastern Wielkopolska regions and the Western Pomeranian and Eastern Pomeranian regions, which is the result of the catchment area being located within reach of a number of physical and geographical regions. The average TA in the winter season ranges from −0.7 to −1.2 • C. On average, there are 40-72 days with a slight frost, starting, at the earliest, around [28][29] November and, at the latest towards, the end of February or in March. Days with very frosty weather sporadically occur and not on a yearly basis [23]. The studied catchment area is characterized by frequent low-intensity rains and unstable snow cover, as well as by the considerable infiltration capacity of soils and numerous depressions with no outlets [40].
The rivers in the region are characterized by a snow-rainfall recharge type. The Noteć River is classified among rivers with a medium developed nival regime, which is attributed to rivers with a mean spring month flow (March-April) constituting 130-180% of the mean annual flow [41]. Its supply regime is characterized by a single fundamental maximum of levels and flows in spring and a minimum in the summer-autumn season. The maximum span of water level fluctuations in the lower reaches is about 3.5 m. The Noteć and the Drawa rivers are characterized by the longest duration of low water periods, which are 75 and 97 days, respectively. The average annual unit runoff from the catchment area for the period of 1971-2010 (according to the Institute of Meteorology and Water Management-National Research Institute (IMGW-PIB) data) in the Pakość section (the upper course of the river) is 3.5 dm 3 /s/km 2 (the average flow rate Q = 5.43 m 3 /h and the minimum flow Q = 0.53 m 3 / h), whereas it is 4.6 dm 3 /s/km 2 (the average flow rate Q = 74.8 m 3 /h and the minimum flow rate Q = 23.2 m 3 /h) in the Nowe Drezdenko section (the lower course of the river). The mean annual unit runoff from the catchment area is 90-100 mm. The Noteć River, as a tributary of the Warta River, has the lowest variability of daily flows and a very low variability of minimal flows, also in the winter season in which IP appear at the end of November/beginning of December and usually persist until the end of March [42]. In the winter season, the river supply and outflow conditions are modified by the occurrence of ice phenomena. The average duration of ice phenomena on the Noteć River is from 31 to 60 days. Permanent ice cover occurs at the end of December and lasts, on average, from 16 to 30 days.
The irregularity of extreme flows for the Noteć River results from the regulation of the river bed, by, among others, the erection of numerous fall stages that considerably modify the natural run-off regime. The mean daily water temperature in the lower reaches of the Noteć River changes seasonally from 3.6 • C (November-March) to 15.7 • C (April-October), reaching 18.8 • C in July. The temperature of flowing waters is sometimes increased by effluent discharges (thermal pollution and municipal and industrial sewage), particularly in the upper and middle reaches [40].

Materials and Methods
The analysis of IP change tendencies was conducted on the basis of daily data concerning the number and nature of phenomena occurring in the river, as well as daily TW and TA values for the period of 1987-2013. Use was made of data from four water gauge stations on the Noteć River (water gauge station symbols: Pakość-NP; Ujście-NU; Krzyż Wielkopolski-NK; and Nowe Drezdenko-ND) and the synoptic station in Piła (symbol: MP) ( Figure 1) obtained from the IMGW-PIB (Warsaw, Poland).
Data are presented in relation to the hydrological year, which lasts from 1 November until 31 October in Poland. The complete RI cycle comprises 8 forms of ice occurrence: stranded ice (SI), frazil ice (FI), stranded ice and frazil ice (SI-FI), ice cover (IC), floating ice (FoI), stranded ice and floating ice (SI-FoI), frazil ice and floating ice (FI-FoI), and ice jam (IJ). TW measurements were performed daily at 7:00 a.m. (GMT + 1) using automatic station probes every 10 min with an accuracy of 0.1 • C (source IMGW-PIB); if these were unavailable, measurements were done by means of mercury thermometers (depth 40 cm) over a period of 5-10 min with an accuracy of 0.1 • C. Measurements of air temperature (the ground surface) were taken at the meteorological stations of the IMGW-PIB by means of electrical sensors or mercury thermometers. The temperature measurement accuracy was 0.1 • C. Mean daily air temperature resulted from eight measurements and was calculated according to the method: (T00 + T03 + T06 + . . . + T21)/8. The hours are given in Universal Time Coordinated (UTC) time (source IMGW-PIB) [38].
During the first stage of research, a determination was made of the IP subset observed at a given station that may be subjected to a reliable analysis of change tendencies over time. These phenomena were taken into consideration in the analysis if they occurred on at least 1% of cool days in hydrological half-year periods in the years 1987-2013 (Table 1). In order to determine the independence of observations, use was made of the Ljung-Box test [43], which makes it possible to indicate whether any group of time series autocorrelation is different than zero. The test statistic is: where n is the sample size,p k is the sample autocorrelation at lag k, and h is the number of lags being tested. Under H 0 the statistic, Q asymptotically follows a χ 2 (h) . For significance level α, the critical region for rejection of the hypothesis of randomness is: where χ 2 1−α,h is the 1 − α -quantile of the chi-squared distribution with h degrees of freedom. Statistic Q(m) of the Ljung-Box test, which is intended to ensure greater strength of the test for completed samples, has an asymptotic distribution of χ 2 for m degrees of freedom. In practice, the selection of maximum delay m has a considerable impact on the result of the test. It is usually conducted for a few different values of m.
During the initial data preparation stage, it was determined that both the time series of annual numbers of IP occurrences and the time series of mean annual TW and TA comprised independent observations ( Table 2). For this reason, it was possible to analyze IP trends and their connections with temperature values by means of techniques that utilized the assumption of observation independence. TA ( • C) X-squared = 6.286 p-value = 0.791 Lag: df = 10 (delay).
The change tendencies of the appearance of IP on the Noteć River and their association with TA and TW were modelled in a set of non-negative integers using models from the general linear model (GLM) group, which are based on sizes with excessive dispersion and zero inflation. Use was made of the logistic regression model for count data with a negative binomial distribution, because the data were characterized by considerable dispersion. The applied logistic regression model makes it possible to interpret results as odds, which are a function of probability. The general logistic model is described by the following formula: where y is expected value of the dependent variable, x is value of the independent variable, exp is exponential function with a basis equal to e that is with a natural logarithm basis (exponent), and b is regression coefficient. Regardless of the regression coefficients and the value of x, the model provides a result in the form of predicted y values in the range from 0 to 1. The models of regression thus defined required the analyzed observations to be independent of each other, which was confirmed-as stated above-through the Ljung-Box test [43].
At the same time, there was considerable zero inflation in the observation series, which would suggest that the process governing the number of occurrences of IP is accompanied by a dichotomous process that determines whether a phenomenon has a chance of occurring in a specific year. In this situation, use was made of the zero-inflated negative binomial model, which is usually applied if data have greater dispersion than indicated by the Poisson distribution, in which variance is equal to the mean, and if their distributions have an over-representation of zeroes. Zero-inflated models are a class of models that are able to cope with additional zeroes [44,45]. Probability distribution for the negative binomial distribution has the following form: where Γ(·) stands for the gamma function, α ≥ 0, j = 0,1,2, . . . ,.
The expected value of distribution is determined by the following formula: where Xi is vector of explaining variables with sizes (1 × k) and β is vector of parameters with sizes (k × 1). Whereas variance is equal to: If α = 0, we receive the Poisson regression model. The statistical significance of the entire model, (the so-called global model significance test) was assessed using the likelihood ratio test (logLik). In the global maximum of the likelihood function, two explanatory variables were used: the variance of the slope χ2 and the p-value.
The necessity of using the regression model for the zero-inflated negative binomial model was confirmed by the Vuong test [46,47], which demonstrated that this model was better suited to data than the negative binomial model alone ( Table 3). The test statistic (Vuong z-statistic) is the relation of logarithmic probabilities of data between two competing models. The analysis of the multi-annual tendency of TW and TA utilized a non-parametric test of the monotonic Mann-Kendall (MK) trend (MAKESENS Excel template) with Yue-Pilon correction (MK-YP) [48]. The MK test is calculated as: where n is the number of data points; xi and xj are the data values in time series i and j (j[i]), respectively; and sgn (xj − xi) is the sign function as: The YP correction was used to decorrelate and detrend the data series in order for the MK test to provide reliable and unbiased results. The modified MK-YP test made it possible to effect an approximate determination of annual, semi-annual, and decadal temperature changes, as well as their changes throughout the entire period. Linear trend strength was estimated by means of the resistant Theil-Sen regression line slope estimator [49]. Analyses were conducted for decades only when data were available for at least 5 years. Years from two decades-the 1980s and the 2010s-were omitted from the research period (1987-2013) because it was not possible to estimate the trend of average annual temperatures using such a limited set of data. The calculations and the elaboration of the results of analyses were performed using the R calculational environment, version 3.3.2 (31 October 2016). The R suite (GNU R) is a software suite for analyzing data, as well as programming platform known as R [50]. The significance results of tests were assessed with reference to the adopted level of α = 0.05. If a different significance level was adopted, relevant information is included next to the calculation results.

Frequency and Thermic Conditions of Occurrence of IP
The analysis of the frequency of days with IP in the Noteć River disclosed the considerable differentiation of occurrences from 490 to 909 days (Table 4). On average, occurrences of the phenomena for between 21 and 40 days were recorded for individual years, with IP occurring more frequently in the middle and lower course of the river. The icing cycle, which encompasses the greatest number of forms of occurrence of ice in the river, was observed only in the middle course of the river (NU water gauge), while at the remaining water gauges, this cycle was incomplete. The RI structure was dominated by phenomena from the first river freezing phase (SI and FI) and a constant IC, whereas phenomena from the final icing stage, i.e., the disintegration of the IC and the appearance of FoI, occurred least frequently. The greatest number of days with IC (251 and 385 days, respectively) was recorded in the upper and lower course of the Noteć River (NP and ND water gauges) ( Table 4). A particular intensification of IP in the river occurred in the decades of 1990-2000 and 2001-2010, while the middle course (NU water gauge) had a period of IP intensification after 2010 ( Table 5). The predominant phenomena were SI, FI, and IC.  An increase in the frequency of the occurrence of IP in the Noteć River was observed in the years 1996, 1997, and 2006, as well as in 1987 in certain instances ( Figure 2). In 1996, their number exceeded 50 days in a year, while in the middle course of the river (NU water gauge), the number exceeded 100. In the winter season of the highlighted years there, an increase in the frequency of occurrence of days with a steadily lower TA (below −10 • C) was observed, and this led to an increase in the number of days with IP ( Figure 3). In the years 1987 and 1996, there was observed the lowest annual TA in the researched period at the station in Piła, which totaled 6.3 and 6.   An analysis of average annual TA trends (station in Piła) and TW trends in the Noteć River (excluding the NK station, for which there were no data concerning TW) for the entire multi-annual period and individual decades disclosed no significant differences (Table 6). In the upper course of the river (NP water gauge) in the period of 1987-2013, a slight decreasing tendency of TW (−0.19 • C) when compared with the growth tendency of TW at the remaining stations (from 0.25 to 0.78 • C) was determined. Explanations: dT-annual trend effect (regression line slope coefficient); dTn-total trend effect for the entire period; T0-estimated average annual temperature at the beginning of the period (regression equation constant); tau-non-parametric Tau-b Kendall correlation coefficient between average annual temperature and the time step; and p-significance of non-parametric correlation (significance of the monotonic trend).
During the cool hydrological half-year period (from November to April) when IP were present in the Noteć River, the average TW in the upper and lower course of the river (at the NP and ND water gauges) was 4.4 • C and was only slightly higher than the TW measured along the middle section (NU profile), whereas the TA in the winter season remained, on average, at a level of 2.3 • C. In the cool half-year period (1987-2013), the TW change tendencies in comparison with the TA trend-which totaled −0.82 • C during the entire researched period -tended to vary (Table 7). Negative TW trends in the cool half-year period were demonstrated by measurement series in the upper course of the river (−0.84 • C at the NP profile) and in the lower course (−0.52 • C at the ND profile). In the middle course of the river (NU profile), the total effect of the trend for the entire period was slightly positive. In the cool half-year periods of the 1990s decade, the TA and TW trends were negative. The TW change tendencies along the Noteć River ranged from −1.37 • C in its lower course (ND water gauge) to −2.60 • C in its upper course (NP water gauge) ( Table 7). TW change tendencies in the Noteć River and TA change tendencies at the MP in the cool half-year periods of the analyzed term are presented in Figure 4.  (Table 7 and Figure 4). Similar changes, albeit of lesser intensity, were recorded in the lower course of the Noteć River (ND water gauge). Meanwhile, the negative TW trend in the winter half-year period persisted throughout all the decades in the middle course of the river (NU profile), with the negative TA trend being maintained throughout the entire researched period and in individual decades.

Linkage of Ice Phenomena Trends With Air and Water Temperature
The analysis of IP tendencies in the Noteć River covered four measurement stations at which different RI cycles were identified. The phenomena occurred irregularly, while in certain instances, they were relatively rare. For this reason, the analysis was conducted only for those IP that, in any given instance, occurred on at least 1% of days of the cool half-year hydrological periods of the researched period.
The assessment of trends of the occurrence of phenomena from first RI phase concerned the following: SI, FI, and SI-FI. The model estimating the number of occurrences of SI in the Noteć River, taking connections with TW and TA into consideration, concerned three water gauges with the exception of ND (lower course) and was, in each instance, globally statistically significant at levels of χ 2 (6) = 13.74 and p < 0.01 (NP water gauge), χ 2 (6) = 27.35 and p < 0.001 (NU water gauge), and χ 2 (4) = 20.62 and p < 0.001 (NK water gauge). This means that the number and annual probability of occurrence were governed by systematic trends. Even though the model was globally significant, in the upper course of the river (NP), none of the individual estimators attained a level of statistical significance. The temporal trend for the zero-inflated model was very close to significance (p = 0.055) and negative (Table 8), which would suggest that in successive years, the probability of occurrence of stranded ice may be steadily greater, though this effect was relatively weak. In the middle course of the river (NU water gauge), the occurrence of SI was connected with the average annual TA, and the connection was negative and statistically significant at the level of α = 0.001 (Table 8). At the same time, the occurrence of a negative temporal trend in the zero-inflated model was determined, which means that the probability of occurrence of SI in successive years may be steadily greater, similarly to in the upper course. In the lower course of the Noteć River (NK profile), due to the lack of data concerning TW, consideration was only given to average annual TA, which was a significant individual change estimator in the model for the number of days with SI (Table 8).  In all instances, the change in trends over time displayed a similar dynamic, i.e., a growth trend in the initial period of observations, that thereafter steadily disappeared ( Figure 5). This was the effect of opposing temporal trends cancelling each other out, i.e., the trend of the number of occurrences and that of the zero-inflated process (NP and NU water gauges). In the case of integrated models, no clear evidence was obtained for the occurrence of any systematic trend, decreasing or increasing.
The model estimating the number of occurrences of FI in the Noteć River, when taking connections with TW and TA into consideration, concerned three water gauges, with the exception of NP, where no such phenomena were recorded (Table 9). For the lower course of the river (ND profile), the model was globally statistically significant at a level of χ 2 (6) = 20.65 and p < 0.001, whereas for the following water gauges-NU (χ 2 (6) = 7.02 and p = 0.135) and NK (χ 2 (4) = 6.43 and p = 0.169)-the model turned out to be statistically insignificant, which means that neither the number nor the annual probability of occurrence of FI were dependent on average annual temperatures and were not governed by systematic trends.   Explanations as in Table 8, significance levels: 0.0001 ***, 0.001 **, and 0.05 *.
In the model for FI at the ND water gauge, important individual estimator included average annual temperature: TA was negatively connected with the occurrence of FI, while TW was connected positively (Table 9). However, no clear systematic trend was noticeable in the observed data and trend graphs ( Figure 6). Explanations: The model for the NU and NK water gauges turned out to be statistically insignificant, and, thus, the distribution of change tendencies has not been presented graphically, with only an interpretation being provided of the number of occurrences of FI in individual years. The model forecasting the number of occurrences of SI-FI, which takes connections with temperature into consideration, was determined solely for the middle course of the river (NU water gauge). Globally, the model was statistically significant for χ 2 (6) = 30.12 and p < 0.001, while its analysis demonstrated a significant and rather clear growth trend of the number of occurrences of the phenomenon over time (Table 10 and Figure 7). However, none of the estimators were statistically significant in the case of the zero-inflated model.  The model for the occurrence of IC in the Noteć River, connected with average TW and TA, concerned all four water gauges. Globally, it was statistically significant at a level of χ 2 (6) = 24.70 and p < 0.001 (NP water gauge), χ 2 (6) = 12.54 and p = 0.051 (NU profile, on the threshold of significance), χ 2 (4) = 27.06 and p < 0.001 (NK profile), and χ 2 (4) = 27.06 and p < 0.001 (ND water gauge). In the upper course of the Noteć River (NP profile), the number of occurrences of IC was associated with average annual TA, which correlated with both the number of occurrences and the probability of occurrence of IC in a given year (zero-inflated process) ( Table 11). The strongest estimator (on the threshold of significance) turned out to be the effect of the growth trend of the number of occurrences of IC, whereas the trend for the zero-inflated process was insignificant but concordant regarding direction with the trend of the number of occurrences. Regarding the IC model, in the middle course (NU water gauge), the low predictive strength was connected with the very low significance of a part of the model for the zero-inflated process. In this instance, occurrences of IC were connected with average annual TW (negative correlation) and time-a significant, albeit slight, growth trend was determined (Table 11).
Concerning the number of occurrences of days with IC in the lower course of the river (NK water gauge), the time trend (negative) turned out to be a significant estimator, while in the case of the zero-inflated process, there was an average annual TA (Table 11). In the lower course of the Noteć River (ND profile), occurrences of IC were associated with the average annual TA, which correlated with both the number of occurrences and the probability of occurrence of IC in a given year (zero-inflated process). The sole estimator that attained statistical significance was the effect of the average annual TW, which turned out to be negatively correlated with the number of occurrences of IC (Table 11). Over time, in the model for the upper course of the Noteć River (NP profile), a slight growth trend of IC became apparent (Figure 8), whereas in the middle course (NU water gauge), the temporal trend started to grow only after temperature fluctuations were taken into consideration. This would suggest that in this instance, the issue of occurrence of the trend was, on the whole, ambiguous and unclear, i.e., the observed slight decreasing trend could have been exclusively caused by temperature fluctuations, while if these fluctuations had not occurred, the number of occurrences of IC in the river could have even increased slightly. In the lower course of the river (NK profile), the disclosed negative trend was practically unnoticeable due to the fact that there was present an opposing trend that bordered on significance in the zero-inflated model. This means that the later the year, the more chance there was of the occurrence of IC, and thus, finally, the two tendencies cancelled each other out. Meanwhile, in the final profile of the lower course of the Noteć River (ND profile), the estimated trend for time alone was, practically, ideally flat, which would point to the lack of systematic changes in the occurrence of IC (Figure 8).
The assessment of trends connected with the disappearance of IP in the Noteć River only concerned occurrences of SI-FoI at the NU water gauge (middle course). At the remaining stations, no phenomena from the third and final phase (disappearance of RI) were recorded for at least 1% of days of the cool hydrological half-year periods in the researched period. The model for SI-FoI in the river was globally statistically significant for χ 2 (6) = 18.72 and p = 0.005. At the same time, none of the individual estimators were significant for either the number of occurrences or the annual probability of occurrence (Table 12). Thus, it is difficult to determine which of the variables were most important in this instance.
The dynamic of changes in the occurrence of SI-FoI displayed a lack of any systematic temporal trend ( Figure 9); however, taking temperature-in this case TW-into account allowed for a better estimation of the fluctuation of the number of occurrences of IP in the Noteć River.

Discussion
Under conditions of the temperate climate, the RI process progresses with varying intensity, more or less regularly in various phases [14,16,32], and this was confirmed by the results obtained for the Noteć River. In Poland, the average dates of appearance of ice phenomena on rivers, as well as the dates of their disappearance, vary considerably [22,24,25,31,42,[51][52][53].
In the present case, a small number of situations (middle course) occurred where RI appeared in three successive phases: the formation of IP, the presence of stable IC, and the disappearance of icing (Tables 4 and 5), which is a rather rare phenomenon in lowland rivers today. Years in which IP did not appear at all in the Noteć River have also been highlighted, and these are indicated in the presently analyzed observation series (1987-2013) by the large number of zeroes. The structure of the analyzed data series-that contained, on the one hand, information about the number of days with icing, and on the other hand, zeros that pointed to the absence of the phenomenon (a considerable dispersion of data with zero inflation)-led us to adopt both the negative binomial model and the zero-inflated binomial model. The correctness of this selection for analyses of IP change tendencies was confirmed by the Vuong test. It was assumed that the process governing the number of occurrences of a given phenomenon is accompanied by a dichotomous process that determines whether a phenomenon has a chance of occurring in a specific year. Already in the years 1878-1934, the number of days without icing in the Noteć River gradually increased, while the lack of stable IC became steadily more frequent [52]; today, this was also demonstrated by observations conducted on the river by the IMGW-PIB (Warsaw, Poland) ( Figure 2).
Numerous studies have stressed that changes in the RI regime are the result of warmer winters, a fact that corresponds well with the positive TA trend in the cool half-year period [14,[54][55][56]; this may, in consequence, translate into positive change tendencies of TW and the disappearance of periods of icing [4,6,16,17,51,57,58]. In the case of the Noteć River, no identical models of the tendencies of the occurrence of IP in the researched period were obtained. This is further confirmed by research into periods of icing in other European rivers, for which no clear direction of change or statistically significant trends have been determined [28,59]. In the temperate climate zone, the RI process tends to intensify in periods that are at least very cool. The formation of stable ice cover is facilitated by days that are very frosty (t < −10 • C) or with a strong frost (t < −20 • C) [23], the occurrence of which is connected with high-pressure weather that ensures the advection of dry and cool continental air masses from the eastern sector [56,60,61]. The impact of changes in global TA on the long-term fluctuations of the dates of formation and disappearance of IC is characterized by spatial heterogeneity [62,63], which means that in individual regions, change trends that are either positive or negative occur [26]. On the whole, the increase in TA and TW in the cool half-year period of the year impacts, among others, the later appearance and earlier disappearance of IP in rivers and, thus, the shortening of the duration of icing. According to Yang et al. [7], the duration of RI is shortened most greatly in the Rocky Mountains, in the north-eastern part of the United States, and in Eastern Europe.
For all water gauges on the Noteć River, the analysis resulted in the development of an identical model of change tendencies for SI; this is the first-next to FI-symptom of icing that forms near river banks as a uniform ice structure that reaches out towards the middle of the river as icing progresses. It is typical of the conditions of insufficient water cooling and calm flow and, thus, of the months of November and December. In all cases analyzed in the study, there was noticeable an initial increase in the number of days with SI, where after, in the second half of the period, the number of occurrences started to fall. TA turned out to be the best estimator for forecasting SI. Because of the negative trend obtained for zeroes in the forecasts, we should assume that it is probable that the occurrence of SI will increase in future. Over a longer period of observation, it was determined that phenomena from the first RI cycle (SI and FI) are dominant in the Noteć River, while stable IC occurs less frequently. Regarding FI, the statistical significance of specific change trends of its occurrence was not confirmed, which may have been due to the nature of its connection with TW (positive impact) and TA (negative impact). This form of icing requires considerable water overcooling [21,23] and the effective discharge of solidification heat. It forms particularly plentifully in the presence of strong, cold winds, even if the decrease in TA is only slight (even at temperatures only a few degrees below 0 • C), which may usually be determined during its first, slight decreases at the beginning of the hydrological year (in November and December). At present, this predominance of the duration of phenomena from the first RI phase in relation to the duration of stable IC is a typical feature of the ice regime of the majority of rivers in Poland [21,22,24,25,31,42,[51][52][53]64,65] and other countries of Europe [5,13,66,67].
In the case of lowland rivers, among which the Noteć River is classified, the phase of IC development is the longest, due to only slight flows and drifts [22,24,53,65], while its formation is facilitated by the long-term persistence of negative TA. The models utilized in the study have demonstrated that the best predictor for forecasting trends of occurrence of IC is TA. Research conducted by Graf and Tomczyk [23] confirmed that ice cover appears on the Noteć River at certain specific threshold values of accumulated negative TA (so-called degree-days), which are different for individual water gauges and equal to −16 • C for NP (upper course), −21 • C for NU (middle course) and NK, and −73.5 • C for ND (the last two stations are situated on the lower course of the river). The more rapid increase in cumulative series of negative TA contributes to an increase in the probability of occurrence and persistence of IC on the Noteć River, though it varies along its course. On average, an increase of one degree-day contributes to an increase in the chance of appearance of stable IC of approximately 1.5% (NU), 2-3% (NP, NK), and 6.0% (ND) [23].
The period of disappearance of stable IC due to the increase in TA occurs in the river in stages, in consequence of which FoI that moves downriver is created; this is known as the period of "ice movement" (drifting of ice) [32,42,52,53]. The drifting of FoI in the Noteć River is usually accompanied by FI and SI phenomena, for which a tendency analysis did not demonstrate a systematic time trend nor confirm the statistical significance of individual predictors. The FoI phenomenon very rarely occurs in the river, while its drift is usually accelerated by the break-up of IC caused by the augmentation of the water level in spring. The accumulation of a considerable quantity of FoI along a small section of the river may lead to the formation of IJ and thus to ice jam flooding [58,68], which, in the case of the Noteć River, is recorded along its middle and lower course (NU and NK water gauges). The phenomenon of IJ in the Noteć River, similarly to FoI, is extremely rare, with the last such occurrence taking place in 2010.
Far-reaching modifications of ice regime features have been indicated by research conducted into the nature and time of occurrence of IP in European rivers for the period of the last four decades [4,16,19,20]. The majority of authors have focused on an assessment of change tendencies at the start and end dates of RI and its duration, which, in light of the recorded increase in TA and TW, usually gives a result pointing to the shortening thereof. Apart from warmer winters, an increase in TA in the autumn and spring periods has also been noticeable, especially during the last twenty years of the 20th century but in recent decades as well [38,39,[69][70][71][72], which delays the date of appearance of IP in rivers and also accelerates their disappearance. Analyses conducted for the last five-to-six decades in Polish rivers [21,22,24,25,31,42,51,64,65,[73][74][75] have confirmed that the duration of IP has been shortened, though their respective change trends have differed along individual river sections. These studies indicated that an increase in TA and TW in the winter hydrological half-year period, which usually corresponds well with the negative IP duration trend, was connected with the variability of atmospheric circulation. In their research on the time of formation and disappearance of ice cover in the rivers Vistula and Niemen,Ćmielewski and Grześ [64] showed that the correlation between air temperature and data concerning the formation and disappearance of IC explains approximately 35-40% of the variability of these phenomena. Marszelewski and Pawłowski [51] confirmed a strong correlation (R 2 from 0.69 to 0.81) between the duration of ice phenomena in the Oder River and the mean air temperature in winter . During the indicated period, the duration of all ice phenomena decreased (by up to 0.58 days·year −1 ), as did the duration of ice cover (by up to 0.46 days per year −1 ). In turn, Wolski et al. [76] showed that during the multi-annual period of 2006-2017, the number of days with an SI occurrence greater than the average could be observed in the middle of the Odra River.
An analysis of circulatory determinants impacting the occurrence of IC in the Noteć River was performed by Graf and Tomczyk [23], while an assessment of the relationship between the TW of Polish rivers and large-scale atmospheric circulation, which was accompanied by the disclosure of TW change patterns in Poland, was carried out by Graf and Wrzesiński [38,39] and Marszelewski and Pius [71,72]. Research conducted by Yoo and D'Odorico [5] in Northern Europe showed that the ice regime of surface waters is connected with changes in the North Atlantic Oscillation (NAO) index, the features of which confirm long-term climate change [16]. The increase in cyclonic activity, caused by (among others) a percentage increase in the positive NAO phase [77], results in warmer winters that start at a later date and are followed by early springs. The increase in the positive phase of the NAO that occurred in 1987-1989 has led to a considerable decrease in the severity of winters in Poland after 1989 [78]. This may have far-reaching consequences for the thermal and ice regimes of rivers, including the Noteć River. Climate changes have resulted in an increase in the frequency of the occurrence of mild winters, low snowfall, and frequent mid-winter thaws in the moderate climate zones [4,57]. The research on the periods of the formation and disappearance of IC in the Vistula and Niemen rivers has shown that the correlation coefficient of the dates of the IC break-up with the NAO index allows for the explanation of about 30% of changes in RI [64]. In the case of the Odra River, it was found that correlation between the duration of ice phenomena and the NAO index (R 2 from 0.42 to 0.48) was weaker than between the duration of ice phenomena and the mean air temperature in winter [51]; the pollution of the Oder River and icebreaking operations on the river were identified as the possible causes. In Poland, the strongest impact on river water temperature in winter  was observed for the positive and negative phases of NAO and AO, when deviations of water temperature from average values were correspondingly higher (by 0.6-1.0 • C) and lower (by 0.4-1.5 • C) [38]. The results indicated that relationships between the temperature of river water and macroscale types of atmospheric circulation are less obvious than in the case of air temperature.
In the case of the Noteć River, which drains the western part of the Polish Lowlands, TW change tendencies in the cool half-year period (1987-2013) were varied in comparison with the negative TA trend-negative in the upper and lower course of the river (NP and ND water gauges) and slightly positive in the middle course (NU water gauge). The results of studies have pointed to the complexity of processes determining the thermal regime of rivers and confirm the impact of additional factors, regional and local, on their temporal and spatial variability. Spatial differences between thermal regime features are typical of rivers with a marked seasonality of supply in the temperate climate [39,59,71,72]. At the same time, no single trend attained statistical significance at the level of individual decades, which may have been due to either of two reasons: obviously, trends for decades were analyzed on the basis of a smaller quantity of data, which automatically translated into lower strength, or the observed growth trend may have been unclear over short periods of time in comparison with random fluctuations. The IP change tendencies for the Noteć River are, to varying degrees, controlled by physical factors, such as TA and TW, and that the effect of the trend for the number of their occurrences is an important predictor that is especially visible in the upper course of the river (NP water gauge). TW turned out to be a significant estimator in the middle and lower course of the river (NU and ND water gauges). This shows that in order to obtain a precise IP forecast, it is necessary to take temporal trend models and models of changes in TA and TW over time into consideration.
Of considerable importance for the persistence of IP in the Noteć River is the transfer of thermal energy, which may occur on two levels: air-water and river bed-water [32]. The process of icing of flowing water commences when TA falls below 0 • C, consequently leading to a decrease of TW. Research into small water courses has shown that TA is the main predictor of changes in their TW [79,80], while the strongest connection occurs between them at equilibrium temperature, which arises when the exchange of heat between air and water is equal to zero [81]. The greatest differences in the intensity of the strength of correlation between TW and TA appear at extremely low values of TA in winter, when TW does not fall below 0 • C and, thus, at the time when IP form and persist in rivers. However, the impact of TA on the thermal characteristics of waters in the majority of rivers is not spatially identical, and the effects frequently differ at various locations, which has also been confirmed by studies performed on the Noteć River. Furthermore, the strength of correlation between TW and TA fluctuates over the year, and relationships may change seasonally depending, for example, on the flow and quantity of subterranean supply; this has been documented in numerous stream temperature models [59,82,83]. Higher temperatures of river water in the winter season, in excess of 2.0 • C, usually point to the presence of pollutants in the river [84,85]. Pollutants have been channeled to the Noteć River for many years, including municipal wastewater from various towns and cities and water from coal mines (from the Konin region) [40]. The thermal and ice regime of the Noteć River has also been strongly modified due to the anthropogenic transformation of the bed and valley of the water-course. In the case of the Noteć River, these are most frequently regulatory activities that lead to changes in bed morphology and the artificial control of water flow [41]. This may impact variations in measurement series of river water temperatures [40], which brings about their instability and, in consequence, changes in ice regime [24].

The performed analysis has demonstrated clear fluctuations in the occurrence of IP in the Noteć
River for the multi-annual period of 1987-2013, which is typical of rivers located in the temperate climate zone. In recent decades, the total lack of icing or periodic changes in its intensity has been observed with a greater frequency, and this holds true for the majority of rivers on the European Lowland, also in Poland. It has been determined that IP in Noteć River occur irregularly and periodically, on average from 21 to 40 days per year, while the RI structure along the river course is very diverse and is dominated by the ice phenomena from the first river freezing phase (SI and FI) and a constant IC. Meanwhile, the ice phenomena from the final icing stage, i.e., the disintegration of the IC and the appearance of FoI, occur least frequently. 2. The analysis performed in order to estimate the IP trends change in the Noteć River has demonstrated a diversified share of the positive and negative trends in the studied time series. An analysis of the average annual TA and TW trends in the Noteć River for the whole multi-annual period and its individual decades have disclosed no significant differences. In the upper course of the river (NP water gauge), a slight decreasing tendency of TW (−0.19 • C) was determined when compared with the growth tendency of TW at the remaining stations (from 0.25 to 0.78 • C). In the cold half-year period (1987-2013), the TW change tendencies in comparison with the TA trend-that totaled −0.82 • C during the whole analyzed period-tended to vary from 0.84 • C (in the upper course) and a slightly positive trend in the middle course of the river to −0.52 • C (in the lower course).
3. The more rapid increase in the cumulative series of negative TA in the Noteć River catchment area contributes to an increase in the probability of occurrence and persistence of IC in that river. However, it is varied along its course (one degree-day contributes to an increase in the chance of appearance by about 1-6%). The probability of occurrence of ice cover depends on the length of time with the advection of cool air masses and, thus, on the duration of the period in which a given type of atmospheric circulation occurs in winter. The negative temporal trend obtained in the zero-inflated model was statistically significant or on the threshold of significance, and this would suggest that in the successive years the probability of IP occurrence may be steadily greater; however, this effect is weak. The result concerned, in particular, phenomena from the first phase, i.e., river freezing (SI). This confirms that in the near and distant future, these may be the only (or predominant) forms of icing in the Noteć River. The analysis has demonstrated the identical model of change tendencies for SI for all water gauges on the Noteć River-there was a noticeable initial increase in the number of days with SI; in the second half of the period the number of occurrences started to fall, and TA turned out to be the best estimator of the observed change. Regarding FI, the statistical significance of specific change trends of its occurrence has not been confirmed, which may be due to the nature of its connection with TW (positive impact) and TA (negative impact). 4. In some of the water gauges, the number of IP occurrences and the annual probability of occurrence were not dependent on thermal conditions, and they were not governed by any systematic trends. This situation may be associated with different nature of relationships existing between TW and TA in autumn (beginning of IP) and spring (the disappearance of icing). Furthermore, in winter, due to the occurrence of RI, the heat transfer relationships between air and water are changing. The inclusion of annual fluctuations and trends of TA or TW makes it possible to improve the estimation process of the level of occurrence and change in the IP trends, which has been demonstrated by the analysis carried out for the Noteć River. The results of research confirm the effectiveness of compilating a few models that are based on sizes with excessive dispersion and zero inflation for the estimation of IP trends in a river. It should be noted that a time trend is a significant estimator for the number of occurrences of days with specific IP in a river.
5. The observed tendencies of change in components of the icing regime of the Noteć River have only been partially confirmed in the recorded winter TA and TW trends. These are frequently observed in rivers that have a distinct seasonality of water supply and are subjected to anthropogenic pressure. In its middle reaches, the Noteć River is canalized and accessible by ships, while its water levels are regulated by a sluice system that probably affects the irregularity of the IP occurrence. The anthropogenic factors have a varying intensity of influence in the catchment area of the Noteć River. The influx of pollutants into the river has been documented in earlier studies, but it was not confirmed by direct analysis of correlations of water temperature, IP, and wastewater [23,40]. This topic was not discussed in this paper due to the lack of adequate data to investigate such relationships.
The results obtained in this study indicate that despite common warming patterns, the individual features of rivers (hydrologic, thermal, morphological, and those associated with human impact) that are identified along its respective courses may determine the pace and future directions of change in the IP occurrence. They also constitute an important source of information for modelling studies and forecasts that provide valuable results that can be applied to water resource management in the periods of, for example, total RI and its consequences (e.g., the rapid drifting and accumulation of FoI and the formation of ice jams and occurrence of ice jam flooding), which are dangerous from the point of view of water management and undesirable for the biotic environment and ecological state of river water.
Funding: This research received funding from an internal grant from the Faculty of Geographical and Geological Sciences, Adam Mickiewicz University in Poznan, in Poland. This research received no external funding.