Climate change and Northern Hemisphere lake and river ice phenology from 1931–2005

At high latitudes and altitudes one of the main controls on hydrological and biogeochemical processes is the breakup and freeze-up of lake and river ice. This study uses 3510 time series from across 678 Northern Hemisphere lakes and rivers to explore historical patterns in lake and river ice phenology across five overlapping time periods (1931– 1960, 1946–1975, 1961–1990, 1976–2005, and 1931–2005). These time series show that the number of annual open-water days increased by 0.63 d per decade from 1931–2005 across the Northern Hemisphere, with trends for breakup and, to a lesser extent, freeze-up closely correlating with regionally averaged temperature. Breakup and freeze-up trends display a spatiotemporally complex evolution and reveal considerable caveats with interpreting the implications of ice phenology changes at lake and river sites that may only have breakup or freeze-up data, rather than both. These results provide an important contribution by showing regional variation in ice phenology trends through time that can be hidden by longer-term trends. The overlapping 30-year time periods also show evidence for an acceleration in warming trends through time. Understanding the changes on both longand short-term timescales will be important for determining the causes of this change, the underlying biogeochemical processes associated with it, and the wider climatological significance as global temperatures rise.


Introduction
One of the main controls on hydrological and biogeochemical processes at high latitudes is the freeze-up and breakup of lake and river ice (Bengtsson, 2011;Rees et al., 2008;Stottlemyer and Toczydlowski, 1999). Ice phenology is gov-erned by the geographical setting (heat exchange, wind, precipitation, latitude, and altitude) and the morphometry and heat storage capacity of the water body (Jeffries and Morris, 2007;Korhonen, 2006;Leppäranta, 2015;Livingstone and Adrian, 2009;Weyhenmeyer et al., 2004;Williams, 1965;Williams and Stefan, 2006). Though preceding surface air temperatures provide a seasonal energy flux that is well correlated with breakup and freeze-up (Assel and Robertson, 1995;Brown and Duguay, 2010;Jeffries and Morris, 2007;Livingstone, 1997;Palecki and Barry, 1986), cycles of temperature linked to large-scale climatic indices have also occasionally been observed to impact ice phenology (Livingstone, 2000).
The majority of lakes and rivers that seasonally freeze are in the Northern Hemisphere, and most research has focused on breakup and freeze-up dates, ice season length, and ice thickness (Duguay et al., 2003;Prowse et al., 2011). As acknowledged by the IPCC (2007), an assessment of changes in broader ice phenology is complicated by, among several factors, the tendency to consider only local areas. Although trends vary, there is a proclivity for breakup and freeze-up records to lean towards shorter ice seasons that are correlated with temperature trends (Table 1). Changes in ice breakup and freeze-up dates, therefore, provide an additional data source for investigating climate patterns (Assel et al., 2003). Whilst the current literature supports observations of a warming climate, the full spatiotemporal variation seen in smaller case studies has not been transferred to a hemispheric scale. This is important because over the next century temperature rise is expected to continue across the Arctic, where lakes and rivers subjected to freeze and thaw cycles are predominantly located (Collins et al., 2013). Understanding historical patterns and changes in lake and river ice phenology is Published by Copernicus Publications on behalf of the European Geosciences Union. 2212 A. M. W. Newton and D. J. Mullan: Climate change and NH lake and river ice phenology required to confidently project future evolution and climate system feedbacks (Brown and Duguay, 2011;Emilson et al., 2018). In the last century the number of ice phenology observations has increased markedly due to their importance for energy and water balances Weyhenmeyer et al., 2011) and infrastructure such as ice roads (Mullan et al., 2017). This paper explores the hemispheric spatiotemporal trends in ice phenology by investigating an extensive database containing 3510 individual time series from 678 Northern Hemisphere study sites. The aim of this work is to use this database to explore how spatiotemporal trends in lake and river ice breakup and freeze-up dates, as well as the number of annual open-water days, have changed across several 30-year-long overlapping time periods from 1931-2005. Sites with data available for the full 1931-2005 time period are used to investigate how short-term trends observed from 30-year-long records compare to longer-term changes. Sites with data for the full 1931-2005 time period are also compared with regional climate drivers (e.g. temperature) to investigate how much of the variability in lake and river ice phenology can be attributed to longer-term regional climate changes.

Materials and methods
The Global Lake and River Ice Phenology Database from the National Snow and Ice Data Center (NSIDC) (available at https://nsidc.org/data/lake_river_ice/, last access: 13 July 2020 - Benson et al., 2013) provides breakup and freeze-up dates for 865 Northern Hemisphere sites. In this database the freeze-up date is defined as the first day in which the water is completely ice covered and the breakup is the date of the last ice breakup before the open-water season. Whilst the specific definitions for breakup/freeze-up may vary between different sites, the precise definition is thought to be consistent at each site. Thus, if climate signals are present in the ice phenology data, then they should still be observable and broadly comparable. This database is supplemented with data from the Swedish Meteorological and Hydrological Institute (SMHI) which contain 749 lakes and rivers using similar terminology. Data for 122 lakes and rivers were provided by the Finnish Meteorological Institute. Several sites were already in the NSIDC dataset but were updated where necessary. The three datasets were integrated to create the Ice Phenology Database (IPD) containing data across North America, Europe, and Russia (Fig. 1). It is important to note that in the later part of the 1980s and 1990s data for many Russian and Canadian sites are not recorded in the database.
Prior to 1931 data are sparse, and many of the longer time series have been explored by Magnuson et al. (2000) and Benson et al. (2012). To investigate the spatiotemporal patterns of ice phenology, five overlapping time periods were studied : 1931-1960, 1946-1975, 1961-1990, 1976-2005, and 1931-2005. These are investigated across three broad areas: North America, Europe, and Russia. All study sites in the database which fall within these time periods and have a maximum percentage of missing values of 10 % were included. These specific time periods were chosen as they offer the opportunity to include as many data from the IPD as possible. Initial analysis showed that of the 1736 lakes and rivers in the IPD, 678 sites had ≥ 90 % of annual data for either freeze-up or breakup for at least one of the time periods within one of the three regions. The number of sites contained within each time period and for each geographical area is shown in Table 2. The final dataset provides 3510 individual time series spread across the Northern Hemisphere ( Fig. 1a) but primarily concentrated in North America (Fig. 1b) and Europe (Fig. 1c). Data on breakup, freeze-up, and annual open-water days for the 1931-2005 time period were available for 87, 48, and 37 sites, respectively ( Table 2). The majority of these sites are clustered around the Laurentian Great Lakes in North America and Sweden and Finland in Europe. In Russia there is only one site in the southwest of Lake Baikal.
Breakup and freeze-up dates were first converted to ordinal days. For some sites, freeze-up or breakup in a specific year occasionally fell in a preceding or succeeding year, and the ordinal date reflects this by providing a relative datei.e. if freeze-up for the 1941 ice season occurred on 5 January 1942, then the ordinal day allocated was 370. Likewise, if breakup for the 1943 ice season occurred on 28 December 1942, then the ordinal date allocated was −3. These records were adjusted as necessary to calculate the number of annual open-water days. The ordinal-day records were tested using the Mann-Kendall test where the null hypothesis of no trend was tested against the alternative hypothesis that there is a monotonic trend in the time series. The Mann-Kendall test is a nonparametric test which detects trends without specifying if they are linear or nonlinear. It does not, however, calculate trend magnitude, so Sen's slope was also used (Yue et al., 2002). A full description of these combined methods can be found in Salmi et al. (2002). These two statistical techniques are commonly used in climate and environmental science as they can account for missing values. These methods were applied to all sites with at least 90 % data coverage (Table 2) for each individual time period to document the significance (α < 0.1), the magnitude of the slope, and decadal change derived from that magnitude. The 90 % allowance means that the maximum number of sites was used for each of the five time periods. The trend magnitudes and directions were converted into the number of days per decade change in the date of breakup and freeze-up or number of annual open-water days at each site during each time period. The magnitude of the decadal change is mapped for all sites, with those that are statistically significant clearly identified in the symbology. To investigate short-term variations over the 75-year time period, residuals were calculated for breakup, freeze-up, and open-water days. Similarly to in Sharma and    Magnuson (2014), a range of running means were applied, with an 11-year window shown to be the most useful for the 75-year time series. A range of climate variables and atmospheric/oceanic modes of variability were downloaded from the KNMI Climate Explorer (http://climexp.knmi.nl/, last access: 31 August 2018) to facilitate examination of potential regional drivers of ice phenology change. Monthly mean temperatures and precipitation were extracted from the Climatic Research Unit (CRU) Time-Series (TS) Version 4.01 (Harris et al., 2014). CRU TS4.01 applies angular-distance weighting (ADW) interpolation to monthly observational data derived from national meteorological services to produce monthly gridded mean temperatures and precipitation at a spatial resolution of 0.5 • latitude × 0.5 • longitude. Wind speed data were extracted from the International Comprehensive Ocean-Atmosphere Data Set (ICOADS), which provides simple gridded monthly wind speeds for 2 • latitude × 2 • longitude grid boxes (Freeman et al., 2017). All these data were downloaded as a spatially averaged regional time series for three geographical regions encompassing only study sites with data for the full 1931-2005 time period: Europe (EUR), 57.5-68.5 • N, 12-29 • E; North America (NAM), 42.5-47 • N, 73.5-95.5 • W; and Russia (RUS), 51.5-52 • N, 104.5-105 • E. Data were extracted for 1931-2005 to correspond with the length of the IPD. We elected for this regionalised strategy because (1) the computational and human resources needed to analyse climate records for each individual site are vast and (2) we were interested in establishing broader regional climate drivers of ice phenology rather than in developing correlations with local climate, which we would expect to be very strong. For 1931, monthly data on the Arctic Oscillation (AO) (Thompson and Wallace, 2000), the Atlantic Multidecadal Oscillation (AMO) (van Oldenborgh et al., 2009), the North Atlantic Oscillation (NAO) (Jones et al., 1997), and the Southern Oscillation index (SOI) (Ropelewski and Jones, 1987) were also extracted.
Ice breakup and freeze-up records from the IPD were spatially averaged into three regional composite records corresponding to the three geographical regions (EUR, NAM, and RUS) defined above. Statistical relationships were then examined between ice phenology dates and climate records (maximum temperatures and modes of variability) using Pearson product-moment correlation. These relationships were analysed on a monthly basis, first for each of the 12 calendar months and second for 12 sliding windows of 3-month means (e.g. mean of January, February, March, then mean of February, March, April).

Results -ice phenology change
A climate regime with increasing mean air temperatures would be expected to increase the number of annual openwater days for sites that seasonally freeze through earlier breakup and/or later freeze-up dates. The decadal trend for the number of annual open-water days allows for an integrated observation of breakup and/or freeze-up date changes relative to each other -i.e. the longevity of open water, rather than a specific shift in the precise breakup and/or freeze-up dates. In this section the results from the Mann-Kendall and Sen's slope analysis are presented for the three main study areas: North America, Europe, and Russia. In total, 678 study sites provide at least one time series with ≥ 90 % complete annual data across the four 30-year time periods and the one 75-year time period, with 3510 individual time series available (Table 2). A summary of the breakup/freeze-up dates available for each of the four 30-year time periods is presented in Fig. 2. These data are used to determine decadal trend directions that have been summarised in Fig. 3 Table 3 as mean changes in breakup and freeze-up dates, as well as in the number of annual open-water days. The general trends are first presented, before looking at the spatiotemporal trends across the three study regions.

General trends
The combined time series and spread of dates for breakup and freeze-up across each time period is summarised in Fig. 2. In North America across all time periods the majority of sites are in a band of latitude between 42-55 • . There is a moderate correlation between median breakup dates and latitude, with the R 2 values typically ≥ 0.50 showing that the breakup date becomes later with increasing latitude (Fig. 2). The one exception to this is for the 1976-2005 time period where the R 2 value is 0.27. However, one site in the northwest of the region has a latitude 16 • more northerly than any other site and appears to skew the correlation as when this outlier is removed the R 2 value increases to 0.48. An additional caveat is that this time period also marks a reduction in the latitudinal range of the sites included in the database. Median freeze-up dates in North America also show a mod-erate correlation (R 2 = 0.49-0.59) with latitude, with freezeup occurring earlier in the year with increasing latitude. Like breakup trends, the 1976-2005 time period shows the weakest correlation, but it is not associated with an anomalous high-latitude site. Unlike North America, where sites cover a wide range of longitudes, in Europe the data are generally restricted to a narrower range in Sweden and Finland (Fig. 1). In all four of the 30-year time periods there is a strong correlation (R 2 = 0.77-0.86) between median breakup dates and latitude (Fig. 2). Freeze-up dates appear to show some association with latitude, but trends are very weak in the first two time periods (R 2 = 0.19-0.21) and weak in the last two (R 2 = 0.39-0.42). The range of breakup and freeze-up dates recorded at European sites (grey points in Fig. 2) becomes more scattered through time, especially south of 60 • N. This shows greater variability in breakup and freeze-up dates at lower-latitude sites and that the time window in which ice breakup and freeze-up occurs appears to have become wider from 1961-2005. These date shifts also show that in the latter two time periods, compared with the first two time periods, there is an increased occurrence of breakup dates within the first 40 d of the year and freeze-up dates shifting to a later part of the winter season -i.e. freeze-up not occurring until January and February of the following year. The wide longitudinal and latitudinal spread of a comparatively small number of lakes in Russia for any time period (Table 2) precludes any confident correlations or associations. Although it is sporadic and not consistent in study areas or time periods, additional analysis of all the lake and river sites show that occasionally median dates were weakly or very weakly (R 2 = 0.05-0.25) correlated with other criteria such as lake area and elevation.
For each 30-year time period the proportions of trends displaying warming and cooling trends have been summarised in Fig. 3. This shows that through time the proportion of sites displaying warming trends has increased. Freeze-up and the number of annual open-waters days display a gradual increase in warming trends through time and an increase in the proportion of sites with statistically significant warming trends. Mean decadal values show a gradual reduction in cooling trends from 1931-1960 and an increased warming during 1976-2005, albeit with high standard deviations when averaged across all sites (Table 3). Despite this consistent pattern, when observed at the three regional scales (discussed below), the proportions of warming and cooling patterns tend to fluctuate between the different time periods. It is only freeze-up changes in Europe that show a similar pattern to that observed for all freeze-up sites when combined, likely reflecting that data in Europe provide a larger proportion of the total number of sites (Fig. 3). What is common amongst all sites is that the 1976-2005 time period displays the largest proportion of sites with warming trends, with the exception of Russia (which has only one site), for freeze-up and the number of open-water days. For breakup the warming pattern for all sites also shows a longer-term increase through time that is interrupted by an increased proportion of sites The data are colour coded by region in the key. The numbers that are adjacent to the recorded dates are R 2 values for each set of regional data. These are also coloured coded on the key -e.g. light blue data show the median breakup dates for North America and an R 2 value between the median date and latitude of 0.52 from 1931-1960. The underlying grey points show the total ranges of dates that were recorded for each site in each time period. Note that some European breakup observations demonstrate that breakup occurred in the December preceding the start of that year's open-water season -i.e. a very early winter cessation of the ice season. Likewise there are sites in all study areas where freeze-up dates were sufficiently late that they did not occur until late in the winter season -i.e. January or February of the following year. displaying cooling trends from 1946-1975). This appears to be largely driven by an increase in the proportion of sites in Europe during that time period displaying either cooling or significant cooling trends. A similar interruption is also observed in North America but is followed during 1961-1990 by a major increase in the number of warming trends. Like trends for freeze-up and the number of annual open-water days, the mean decadal change for all sites shows warming trends develop and increase in magnitude by 1976-2005, again with the caveat that the standard deviation is high enough to switch the trend direction (Table 3). The limited number of Russian sites with breakup data show a decrease through time in the proportion of cooling trends (Fig. 3).
For breakup, freeze-up, and annual open-water days there is general pattern towards warming through time and mean values increase in the magnitude of change. This increase in magnitude is sufficient so that during 1976-2005 breakup was 2.81 d per decade earlier (σ = 2.18) and the number of annual open-water days increased by 5.83 per decade (σ = 4.08) for all sites. The standard deviation from these sites is lower than the mean magnitude of change, meaning variation higher than 1 standard deviation is required to potentially move across a zero value and change the trend direction -i.e. whilst the standard deviation is larger than most other time periods, the higher magnitude means that more of this variability is in one trend direction (Table 3). A difference is also observed for the evolution of lakes and rivers, where rivers appear to show a more consistent warming pattern for breakup, freeze-up, and the number of annual openwater days through time (Table 3).

North America
In North America, the only sites with consistent data are clustered around the Great Lakes. During 1931During -1960, in the east, earlier breakup dates dominate, and in the west, later breakup ( Fig. 4a), with a number of sites being statistically significant (Fig. 3). This variation explains the large standard deviation (σ = 2.56) of the mean trend towards a 0.36 d per decade earlier breakup (Table 3). An east-west pattern is reversed in the 1946-1975 period, with later breakup more common in the east (Fig. 4d). Mean trends show breakup dates be-   (Table 3), many of which are statistically significant (Fig. 3). From 1961From -1990, most sites display earlier breakup trends, with a mean change of 2.98 d per decade (σ = 1.88) (Table 3). Nearly half of all sites display significant breakup trends (Figs. 3, 4g), many of which previously displayed significant later breakup trends (Fig. 4d). Four sites show later breakup trends, of which one is geographically isolated and the others are surrounded by lakes with earlier breakup trends, many of which are significant. This suggests local fac-tors, such as human modification of water courses  or lake circulation patterns (Bennington et al., 2010), might account for local-scale heterogeneity. From 1976-2005 sites are clustered around the Great Lakes and demonstrate partial changes compared to the preceding time period (Fig. 4j). Whilst 72 % of sites trend towards earlier breakup (Fig. 3), in the east several sites now display low-magnitude earlier and later breakup trends. Earlier breakup decadal change for lakes, at 1.16 (α = 1.39), is double that for rivers, at 0.56 (α = 0.81) ( Table 3). The standard deviation continues to show considerable variation around the mean. Fewer sites with freeze-up data are available compared to breakup (Table 2) and remain generally clustered around the Great Lakes (Fig. 4b). From 1931From -1960 no clear geographical pattern exists, with 25 % of sites displaying significant later freeze-up trends for rivers and lakes (Fig. 3). Mean decadal trends show freeze-up was 0.85 d per decade later, but this is associated with a high standard deviation (σ = 3.16) and a large difference in the mean trends for lakes and rivers (Table 3). During 1946During -1975, spatial patterns remain varied (Fig. 4e) and sites with significant later and earlier freeze-up trends each account for 10.5 % of all sites (Fig. 3). Significant sites are both rivers and lakes and unlike breakup do not appear to be clustered east of the Great Lakes. The mean trend for lakes remains low at 0.17 d per decade earlier (σ = 2.29), whilst rivers are comparably higher with freeze-up occurring 1.22 d per decade (σ = 4.42) later (Table 3). Freeze-up date changes during  show that sites in the west more commonly trend towards earlier freezeup and in the east towards later breakup (Fig. 4h). Compared with the breakup trends for the same period and freeze-up trends for the preceding period, the proportion of sites with significant earlier (4.9 %) and later (3.3 %) freeze-up dates is smaller (Fig. 3). The mean decadal trend of 0.18 d per decade (α = 1.97) earlier freeze-up dates for lakes and rivers combined is weaker than observed for earlier breakup during the same period (Table 3). From 1976-2005, freeze-up trends demonstrate a clear pattern, with no sites displaying earlier freeze-up trends (Fig. 4k) and 39 % of sites showing significant later freeze-up trends (Fig. 3). This is markedly different to all other time periods where spatial patterns were much more varied in the Great Lakes region (Fig. 4h). There are no river sites with freeze-up data for this time period (Table 2), and mean values for lake changes show that freeze-up was becoming later by 3.61 d per decade (α = 2.32) ( Table 3).
Trends for annual open-water days during 1931-1960 are broadly similar to those for freeze-up, with a comparable number of sites showing more or fewer open-water days (Fig. 4c). Of the 17 sites, 4 show significant trends (note that 2 sites overlap in Fig. 4c), and this variability reflects the low mean value of 0.29 fewer annual open-water days per decade (σ = 4.82) and is mostly associated with lakes (Table 3). From 1946-1975 the number of annual open-water days closely matches breakup trends, with 20.7 % of sites displaying significant trends towards a decrease (Fig. 3), all of which are east of the Great Lakes (Fig. 4f). Reduced annual openwater days are observed for lakes rather than for rivers, which display a mean increase (Table 3). Annual open-water days during 1961-1990 are similar to breakup patterns during the same period, including in western Canada where freeze-up dates were earlier (Fig. 4i). The low magnitude of freezeup trends compared to high-magnitude breakup trends in the same area has a larger impact on the number of annual open-water days. The majority of sites trend towards more open-water days, with 26.3 % being significant (Fig. 3) and spread across North America (Fig. 4i). The mean magnitude of change shows the number of annual open-water days increased by 2.77 d per decade (α = 3.12), with the changes for lakes being larger in magnitude than for rivers (Table 3). Most sites with data for the number of annual open-water days in the preceding time period show the trend direction changed or reduced in magnitude, even when the 1946-1975 trend was a significant reduction in open-water days ( Fig. 4f and i). Patterns from 1976-2005 reflect that most sites display earlier breakup and later freeze-up dates, extending the length of the open-water season by 4.15 more days per decade (α = 2.84) (Table 3, Fig. 4l). In total, 36.8 % of sites display significant trends towards more open-water days (Fig. 3), maintaining warming trends from the preceding time period but with less variability in the magnitude of that change.

Europe
In Europe, 1931-1960 breakup trends show a proclivity for sites to display non-significant earlier breakup or no trend at all (Fig. 3). Most sites trending towards earlier breakup dates are at higher latitudes compared to those displaying later breakup (Fig. 5a). The lack of observable trends is reflected by the low magnitude of the mean trend towards earlier breakup by 0.10 d per decade (α = 1.29) for lakes and rivers (Table 3). In 1946-1975 most sites show later breakup dates by 1.75 d per decade (α = 1.31) (Table 3), with the only observable spatial pattern being that of the 22.1 % of sites displaying significant later breakup trends (Fig. 3); most are located in areas where earlier breakup was common in the preceding time period (Fig. 5d). By 1961-1990 decadal breakup trends switched from predominantly later to earlier breakup. Of the 261 sites, 53.6 % display earlier breakup, with a further 8.8 % being significant (Fig. 3), with a change towards earlier breakup dates by 0.82 d per decade (α = 1.25), but the variability remains large enough that 1 standard deviation of change is enough to switch the trend direction (Table 3). Northern sites make up the majority with significant earlier breakup trends for both lakes and rivers (Fig. 5g). There remains spatial variability, with 12.6 % of sites showing later breakup trends. The magnitude of the trend towards earlier river breakup dates is almost 3 times that of lakes (Table 3). From 1976-2005 most sites display earlier breakup trends (Fig. 5j), of which 72.3 % are significant (Fig. 3). During this period the breakup date has become earlier by a mean of 3.70 d per decade (α = 2.00), with the magnitude of change experienced in lakes over double that for rivers (Table 3).
During 1931-1960, a total of 45.2 % of sites display earlier freeze-up, with a further 22.6 % being statistically significant (Figs. 3, 5b). Freeze-up decadal trends show lake freeze-up became earlier by 2.31 d per decade (α = 3.43) ( Table 3). The large standard deviation reflects highly variable trend magnitudes towards both later and earlier freeze- The trend directions and magnitudes were derived from the Mann-Kendall and Sen's slope tests. The triangles and circles indicate whether the trend was or was not statistically significant. Sites with a dot in the centre of the circle are river sites. Thus, a red triangle symbol with a dot in the middle indicates a river site that has a statistically significant warming trend over that time period. The blue and red tones on the scales are related to cooling and warming trends, respectively. Note that in some places the symbols overlap. up (Fig. 5b). The five river sites trend towards later freeze-up dates by 3.52 d per decade (α = 3.17). From 1946From -1975, spatial patterns in southern Finland (Fig. 5e), where many sites previously displayed significant earlier freeze-up dates, there is now considerable variability, more so than for breakup (Fig. 5d), with both earlier and later significant freeze-up trends. Compared to 1931-1960 there is a considerable drop in the number of sites displaying significant earlier freezeup trends to 8.4 % (Fig. 3). Mean lake decadal trends show earlier freeze-up reduced to 0.78 d per decade (α = 3.27) but with considerable variation (Table 3). Rivers continue to have opposing trends but also experienced a reduction in trend magnitude. During 1961During -1990 there is a clear increase in sites displaying later freeze-up trends and a reduction in trend magnitude for sites showing earlier freeze-up (Fig. 5h). Both freeze-up and breakup trends in Sweden display a warming pattern, whilst in Finland they are generally opposed (Fig. 5g, h). The decline in earlier freeze-up lake trends is now characterised by a later freeze-up of 0.34 d per decade (α = 2.17) ( Table 3). In the final time period the region is characterised by later freeze-up trends (Fig. 5k), which are similar to breakup trends (Fig. 5j). Later freeze-up trends account for 52.9 % of sites, with another 21.5 % displaying significant later freeze-up. A small number of sites display significant earlier freeze-up trends, but these are out of synchrony with the wider area (Fig. 5k). This time period is the culmination of a gradual reduction in earlier freeze-up trend magnitude for lakes during 1931-1960, before a switch to later freeze-up dates and then a magnitude increase in later freeze-up dates to 2.51 d per decade (α = 3.05) ( Table 3). Through all four time periods rivers have displayed trends towards later freeze-up dates (Table 3).
Spatial patterns in the number of annual open-water days from 1931-1960 (Fig. 5c) are similar to those observed for freeze-up, with most sites displaying decreases (Fig. 3). Across all sites a mean reduction of 2.09 d per decade (α = 4.06) is associated with considerable variation, whilst lakes and rivers show opposing trends (Table 3). During 1946During -1975 open-water days (Fig. 5f) remain broadly similar to freeze-up trend patterns for the same period (Fig. 5e), albeit with local-scale changes that appear to be associated with significant later breakup trends in southern Finland (Fig. 5d). The proportion of sites showing fewer open-water days remains broadly the same, as do mean trend values (Fig. 3, Table 3). The increased trend magnitude for river open-water days is halved compared to the previous time period, but this reflects the fact that only two river sites have data. Spatial patterns for open-water days during 1961-1990 (Fig. 5i) closely resemble breakup (Fig. 5g), except for in southern Finland where earlier freeze-up trends (Fig. 5h) cause several sites to display fewer open-water days. Most sites show an increase in open-water days (Fig. 3), with a mean increase of 1.81 d per decade (α = 2.84) ( Table 3). From 1976From -2005 in the number of annual open-water days are similar to breakup trends (Fig. 5l, j), with a near-uniform increase and 50.7 % of sites significant (Fig. 3) (Table 3), with the trend being considerably stronger for lakes than for rivers.

Russia
In Russia there are only a few sites across the four 30year time periods with breakup, freeze-up, or open-water-day data, with the 1976-2005 time period only having one site at Lake Baikal (Table 2). The majority of the data are clustered in northwest Russia, with a number of individual sites spread out across the Kazakhstan border region and around Lake Baikal in the east (Fig. 6). The lack of spatiotemporal consistency makes it difficult to determine any prevailing trends. Broadly there is a reduction in the number of sites displaying later breakup dates through time (Fig. 3), as is also reflected by the changes in mean breakup date from 0.83 d per decade (α = 0.79) later in 1931-1960 to 0.83 d per decade (α = 1.83) earlier from 1961-1990 (Table 3), albeit with the latter associated with more variability. For breakup trends, in the northwest there are two sites with continuous data across the first three time periods (Fig. 6a, d, g), and these show a gradual change from later to earlier breakup through time. The adjacent sites in this area also show a tendency towards earlier breakup dates during these time periods. The border region sites generally display earlier breakup dates, many of which are statistically significant during the 1946-1975 time period. Around Lake Baikal there is considerable variation between different sites, with no dominant trends, even for the one continuous site through all four 30-year time periods (Fig. 6j).
Between the four 30-year time periods, sites with freezeup data covering at least two time periods demonstrate considerably more variation than those with breakup data (Fig. 6b, e, h). Between different time periods the freezeup dates for the same sites can move in opposing directions, and in some cases, such as in the Kazakhstan border region, these freeze-up date changes have been significant. Longterm there is an apparent reduction in the number of sites displaying earlier freeze-up trends (Fig. 3), but this is caveated by the low number of sites with data and the much larger standard deviations associated with decadal trends (Table 3). Changes in the number of annual open-water days across Russia capture a slightly more consistent pattern compared to changes in breakup and freeze-up dates, but it remains spatially chaotic, with no dominant spatial patterns observable (Fig. 6c, g, i). In all three regions, northwest Russia, the Kazakhstan border region, and around Lake Baikal, there is a shift through time for most sites with continuous data to display more annual open-water days per decade, a number of which are statistically significant (Fig. 3). However, these values are again associated with considerable variation around what is generally a low-magnitude decadal mean (Table 3). The one site with continuous data through all four time periods, Lake Baikal, shows a gradual switch from fewer annual water days during the first time period to no observable trend, before demonstrating more open-water days in the final two time periods, suggesting a gradual warming signal.  (c, f, i, l) in Europe for the four individual time periods. The trend directions and magnitudes were derived from the Mann-Kendall and Sen's slope tests. The triangles and circles indicate whether the trend was or was not statistically significant. Sites with a dot in the centre of the circle are river sites. Thus, a red triangle symbol with a dot in the middle indicates a river site that has a statistically significant warming trend over that time period. The blue and red tones on the scales are related to cooling and warming trends, respectively. Note that in some places the symbols overlap.

Sites with continuous data -1931-2005
Data covering the full 1931-2005 time period in North America are clustered around the Great Lakes region. Over this period mean breakup dates became earlier by 0.66 d per decade (α = 0.50) (Table 3), with 66 % of sites displaying earlier breakup trends and 29.5 % showing significant earlier breakup dates. No dominant spatial patterns are observed, with earlier breakup dates observed across the entire Great Lakes region, except for two sites displaying no trend (Fig. 7a). The extent of sites with freeze-up data limits spatial analysis, but of the nine sites with data, 55.6 % show statistically significant later freeze-up dates (Fig. 7d), with freeze-up, on average, occurring 0.84 d per decade (α = Figure 6. Decadal trends for breakup (a, d, g, j), freeze-up (b, e, h, k), and the number of annual open-water days (c, f, i, l) in Russia for the four individual time periods. The trend directions and magnitudes were derived from the Mann-Kendall and Sen's slope tests. The triangles and circles indicate whether the trend was or was not statistically significant. Sites with a dot in the centre of the circle are river sites. Thus, a red triangle symbol with a dot in the middle indicates a river site that has a statistically significant warming trend over that time period. The blue and red tones on the scales are related to cooling and warming trends, respectively. Note that in some places the symbols overlap. 0.78) later through time (Table 3). Of sites with both breakup and freeze-up data (Fig. 7g) (Fig. 4) appear to be superim-posed onto a longer-term warming pattern, particularly the cooling trend towards later breakup dates from 1946-1975 (Fig. 4d, Table 3). Breakup dates, when viewed as a running 11-year annual mean (Fig. 7a), show a weak (R 2 = 0.25) trend towards earlier breakup, whilst freeze-up trends display a moderate trend (R 2 = 0.48) towards later freeze-up (Fig. 7a, d). Breakup and freeze-up trends combined show The triangles and circles indicate whether the trend was or was not statistically significant. Sites with a dot in the centre of the circle are river sites. Thus, a red triangle symbol with a dot in the middle indicates a river site that has a statistically significant warming trend from 1931-2005. The blue and red tones on the scales are related to cooling and warming trends, respectively. Note that in some places the symbols overlap. that once shorter-term variability is removed, there is a moderate trend towards more annual open-water days per year (R 2 = 0.50) (Fig. 7g). Sites for the 1931-2005 period in Europe cover much of the length of Sweden and southern Finland. Breakup date changes over this time period suggest that it was becoming 0.52 d per decade earlier (α = 0.40) ( Table 3). Most of the significant trends are located in southern Finland, with Sweden characterised by sites with low-magnitude earlier breakup or no observable trend (Fig. 7b). Lakes and rivers both trend towards earlier breakup, albeit with rivers displaying a lower magnitude. Freeze-up trends demonstrate greater variability than breakup ones, with freeze-up dates becoming earlier by 0.20 d per decade (α = 0.97) (Table 3). However, lakes and rivers show opposing trends, with rivers demonstrating later freeze-up by 0.70 d per decade (α = 0.70) (Table 3). Spatial patterns in freeze-up dates vary more than breakup dates, with significant trends towards earlier freezeup dates in Finland and later freeze-up in Sweden (Fig. 7e). These heterogeneous changes in freeze-up dates are also reflected in the low-magnitude mean trend of 0.39 more openwater days per decade (α = 0.90) ( Table 3). The spatial patterns remain varied but are more similar to those of freeze-up dates (Fig. 7h). For all three phenomena, the 11-year running means of the residuals display very weak correlations through time.
In Russia only Lake Baikal has data for 1931-2005, and these show there is no observable change in breakup dates (Fig. 7c), in contrast to the freeze-up dates which have become significantly later by 1.04 d per decade (Table 3, Fig. 7f). Unsurprisingly, when the two are combined there is a significant trend towards 1.53 more annual open-water days over the 75-year time period (Fig. 7i). The 11-year running means of the residuals show strong trends towards later freeze-up (R 2 = 0.60) and moderate trends towards more open-water days per year (R 2 = 0.38).

Results -causes of ice phenology change
Correlations are investigated between ice phenology dates and a series of regionally averaged climatic variables and indices for each of the three study regions (Fig. 8) -Europe (EUR), North America (NAM), and Russia (RUS) -on a monthly basis and for 3-monthly means over the time period 1931-2005. Unsurprisingly, rising temperatures appear to be the dominant control on the shift towards earlier breakup and later freeze-up in the ice phenology records. Late winter and spring temperatures negatively correlate most strongly with breakup, which is expected since rising temperatures lead to more rapid ice melt and thus earlier breakup dates. Autumn and early winter temperatures positively correlate most strongly with freeze-up, which is entirely as expected as increasing temperatures lead to delayed freeze-up dates. In Europe and North America, the month preceding breakup (April and March, respectively) exhibits the strongest correlation with temperatures, whereas for freeze-up the strongest correlation with temperatures occurs in the month of freezeup (November and December, respectively). This may relate to the gradual build-up of rising air temperatures required to break up ice to depth, as opposed to the more rapid onset of freeze-up with falling autumn and winter air temperatures. The 3-month temperature means exhibit even stronger correlations with breakup and freeze-up, with March-May temperatures and February-April temperatures correlating most strongly with breakup in Europe and North America, respectively, and October-December temperatures correlating most strongly with freeze-up in both Europe and North America. These correlations are physically sensible, with breakup and freeze-up occurring towards the end of the 3-month means. In Russia, the strongest correlations with breakup occur in February -3 months prior to the mean breakup date in early May, which may relate to an increased ice thickness and hence longer time period required to cause breakup. However, when considering the 3-month temperature means, the strongest correlations with breakup occur during February-May -which fits more closely with the mean breakup date. Temperatures during the month preceding freeze-up (December) and particularly the 3-month mean period October-December correlate most strongly with freeze-up dates in Russia. This delayed response to falling winter temperatures in Russia compared to in Europe and North America may relate to the influence of other climatic or site-specific factors, especially since the Russian record applies to only one lake.
Although temperature exhibits the strongest correlations with both breakup and freeze-up, precipitation also appears to play an important role in some instances. Increasing winter precipitation (January and particularly the January-March mean) is associated with earlier breakup in Europe, while increasing spring precipitation (March and particularly the March-May mean) appears to exert a stronger influence on earlier breakup in North America. The latter likely relates to increasing precipitation as rainfall, which aids in the melting of ice (Beltaos and Burrell, 2003). The rising winter precipitation in Europe, presumably as snowfall, may also be associated with earlier breakup since snowfall settling on ice can insulate the ice surface and prevent further thickening during the winter (Park et al., 2016) -therefore potentially promoting earlier breakup. Rising precipitation in November (and to a lesser extent the November-January mean) is associated with later freeze-up in Europe. This may relate to increased discharge into lakes or rivers, making it harder for surfaces to stabilise and freeze. The correlations between precipitation and freeze-up are weak in both North America and Russia, while Russia also exhibits weak correlations between precipitation and breakup. There are also some relatively close associations between wind speed and breakup and freezeup in Europe and North America (no wind speed data were available for Russia). Higher wind speeds in summer correlate most strongly with later breakup and earlier freeze-up in North America. The latter seems counter-intuitive since high wind speeds are generally thought to disrupt the water surface and delay freeze-up, while the former does not have any particularly relevant temporal connection. These correlations are not particularly strong compared to those of temperature with breakup and freeze-up and to a lesser extent precipitation; thus, they could be a product of chance that relates to false positives.
In terms of the atmospheric and oceanic modes of variability, some strong correlations exist with breakup and to a lesser extent freeze-up in all regions. Most notably there are strong negative correlations between breakup and winter/early spring NAO and AO, i.e. when NAO and AO are in a positive phase, breakup occurs earlier. This is particularly true in Europe, where a strong positive phase of NAO and AO for the January-March mean and the February-April mean, respectively, is associated with earlier breakup. Correlations for Russia at a similar time of year are also apparent, while correlations in North America are much weaker. Positive correlations (albeit not as strong) between freezeup and NAO/AO occur in autumn in Europe and early winter in North America; i.e. when NAO and AO are in a positive phase, freeze-up occurs later. These findings are expected, since a stronger positive NAO/AO phase results in an increase in stronger westerly winds, drawing warmer air across northern Europe fed from the North Atlantic Drift and the Norwegian Current (Hurrell, 1995). A strong positive NAO/AO promotes later freeze-up in late autumn/early winter and earlier breakup in spring. Trends towards earlier breakup and later freeze-up throughout the latter third of the 20th century may relate to the positive trends of the NAO and the closely associated AO for much of the 1970s and 1980s, with historical highs in the early 1990s (Cohen and Barlow, 2005). Correlations with AMO and SOI for the full time period are generally not as strong, with the exception of negative correlations between late spring AMO and breakup in North America; i.e. when AMO experiences a warm phase, earlier break up occurs. During warm phases of the AMO, elevated sea surface temperatures in the North Atlantic bring about warmer and drier conditions across much of North America (Enfield et al., 2001) -hence the association between earlier breakup with the AMO in this region. There are also positive correlations between winter and spring SOI and freeze-up in Russia (i.e. when SOI experiences a positive phase, later breakup occurs).
All correlations were established between local ice phenology records and the broad regional climate for each region rather than the local climate corresponding to each site. Examining the latter, while more labour intensive, would likely reveal stronger correlations on a site-by-site basis -acknowledging the fact that synoptic and local climate forcings can greatly influence the timing of lake and river ice freeze-up and breakup. The broader geographical approach we have taken also has clear merit, however, as it demonstrates that wider regional climate exerts considerable influence over ice phenology. We also acknowledge the potential for "falsepositive" correlations when assessing so many correlations in a matrix as we do in Fig. 8. This provides reason to be cautious when interpreting these findings.

Discussion
The results presented for all three regions show that between different 30-year time periods there are fluctuations in the trend directions for breakup and freeze-up dates, as well as for the number of annual open-water days. The two most recent 30-year time periods in North America and Europe (Figs. 4,5) show that warming trends dominated. Warming trends for the number of annual open-water days were initially driven by earlier breakup dates before then being increased further by later freeze-up (see below). This is in line with other studies that capture long-term reductions in the ice season (Futter, 2003) and show that warming breakup trends are more common (Brammer et al., 2015;Jensen et al., 2007), and whilst freeze-up trends do move towards warming patterns, they are often more variable Hewitt et al., 2018). The 30-year time period analyses documented here also show that some short-term variations lead to variable spatial patterns through time. For example, in the Great Lakes region, the two most recent time periods show consistent trends towards earlier breakup dates (Fig. 4g, j), which is corroborated in more localised studies (e.g. Magnuson et al., 2005), but the trends in the first two 30-year time periods show variability in the trend magnitude and direction of phenology changes, with sites to the west and east displaying opposing trends (Fig. 4a, d). The trends in this region, as well as more broadly across North America, Europe, and Russia, are dominantly driven by regionally averaged temperature changes, with precipitation and teleconnections also helping to explain some of the variation (Fig. 8). Such a finding is not new (e.g. Blenckner et al., 2004;Duguay et al., 2006;Ghanbari et al., 2009;Hewitt et al., 2018;Livingstone, 1999;Sharma et al., 2013;Smith, 2000), but it does further confirm that the prevailing climate conditions can only partially account for some of the variability in ice phenology trends. What is interesting is the fact that strong correlations can be established despite correlating ice phenology records with spatially averaged climate data over large regions. This indicates that broad regional climate exerts a considerable degree of influence over changes in ice phenology. Such a finding is important because it means that site-specific climate data, either from in situ observations or from numerical downscaling of climate models, may not be explicitly required to explain a large amount of the phenology variation. Whilst there are clearly merits in looking at sites with local data to better understand the underlying processes taking place and how this relates to the observed climatological trends (e.g. the influence of wind on ice phenology), being able to regionalise and simplify the analysis to sites across broad areas that do not have local climate ob-servations is important for upscaling efforts to project largerscale climatic changes.
Across the longer 75-year time period the results broadly match those previously published (Table 1) and show general warming patterns for breakup across all regions (Fig. 7). Freeze-up patterns in Europe show less consistent patterns through time, with many sites showing earlier and later freeze-up trends (Fig. 7e). Whilst these freeze-up trends do evolve into warming trends in the latest time period, this is not fully captured in the 75-year time period studied here, but it is documented in other studies looking at longer records (e.g. Korhonen, 2006). It is also notable that the standard deviations of the trends derived from the Mann-Kendall and Sen's slope analyses for freeze-up tend to be higher than those for breakup (Table 3). Although temperature appears to be able to explain a large proportion of variations in freezeup dates, at least in Europe, it does not account for why the variability is larger than for breakup, which is also well correlated with temperature changes (Fig. 8). Whilst breakup is dominated by thermal characteristics of the climate, freezeup is a result of not just the thermal properties of the environment but also water kinetics -e.g. even if water temperatures are low enough to freeze, wind and water movement can mechanically prohibit freeze-up as the kinetic energy makes it harder for the water to stabilise or ice patches to agglomerate (Beltaos and Prowse, 2009). The complexity involved in water freeze-up likely acts as an important control on these fluctuating trends and would benefit from additional study to explore how this can be accounted for in models (e.g. Bruce et al., 2018). This likely explains why the breakup and freezeup patterns do not simply reflect observed increases in air temperatures.
Unlike for lakes, and with the exception of European river breakup trends from 1946-1975, the mean ice phenology trends for rivers show a more consistent warming pattern through most time periods for all regions (Table 3). Whilst acknowledging the caveats of a limited number of sites, the above evidence suggests that during the 20th century, rivers were responding to increased surface air temperatures faster than lakes. This may be explained, possibly, by the river flow gradient causing waves and ripples which instigates air turbulence and greater interaction of water and air causing a faster transfer of atmospheric heat. Whilst ripples and waves do form on still water bodies, this is likely limited compared to on actively flowing rivers, causing a slower response time in lake temperatures to air temperature increases. As the lakes gradually experience this warming the same reasons may also restrict heat exchange from the lake to the atmosphere. Though the physics require further study, it is possible this thermal legacy allowed lakes to gradually become a heat sink and might explain why over longer timescales the lakes begin to demonstrate larger-magnitude warming trends than rivers, particularly in the 1976-2005 time period (Table 3). Changes in the number of open-water days may relate to movements in breakup and/or freeze-up dates, allowing the relative influence of date changes to be compared. Figure 9 summarises sites with open-water data across all time periods in each region and separates breakup and freeze-up combinations into warming, cooling, and no trends -e.g. 35.2 % of North American sites during 1931-1960 had earlier breakup, later freeze-up, and more open-water days. In all three regions there is a gradual reduction through time in the proportion of sites displaying fewer open-water days caused by later breakup and earlier freeze-up dates. Most other sites are characterised by showing the same trend direction towards earlier or later dates, where either later breakup or earlier freeze-up trends (cooling trends) are stronger than later freeze-up or earlier breakup trends (warming trends), thus reducing the relative number of open-water days (Fig. 9). Through time there is a reduction in the number of sites displaying significant trends towards fewer open-water days, which contrasts with an increasing proportion of sites displaying more openwater days, with the trends at many sites becoming signif-icant in the later time periods 9). Most changes appear to be dominated by sites with both earlier breakup and later freeze-up dates or where earlier breakup or later freeze-up trends are larger in magnitude than later breakup and earlier freeze-up ones (Fig. 9). Some anomalous sites with no warming breakup or freeze-up trends relate to lowmagnitude trends close to zero. Considering all sites combined, most that display trends towards more or fewer openwater days do so because both breakup and freeze-up date trends are moving in opposite directions -earlier breakup and later freeze-up in the case of more open-water days, and the opposite trends where there are fewer open-waters days -suggesting changes across different seasons. During 1931-2005 most sites display an increase in open-water days that is predominately driven by earlier breakup and later freeze-up dates in North America. This aligns well with other studies looking at a range of different sites across the region showing ice season length was driven by either earlier breakup (Brammer et al., 2015;Futter, 2003) or both earlier breakup and later freeze-up (Latifovic and Pouliot, 2007). In Europe the pattern is more mixed with a number of sites showing that earlier freeze-up trends are enough to reduce the number of open-water days for ∼ 25 % of sites, irrespective of a warming pattern in earlier breakup dates. In some circumstances the ice-free season shifts -e.g. 17.3 % of sites in Europe during 1931-2005 display earlier breakup and earlier freeze-up -without actually changing its length, potentially having consequences on biogeochemical cycles in areas that have lakes responding at different rates and in different trend directions. The majority of sites do, however, display trends towards more open-water days but from a range of different breakup and freeze-up trend combinations, with most related to earlier breakup and later freeze-up dates (Fig. 9). This is similar to observations from Finland looking at a longer time period and documenting reduced ice season lengths (Korhonen, 2006). The one Russian site shows more openwater days being driven by later freeze-up dates. Combined, these results match well with numerous other studies across the Northern Hemisphere showing earlier breakup and later freeze-up dates (e.g. Magnuson et al., 2000; Table 1), with the addition that the strength of these changes has increased in more recent times and that the changes in the number of open-water days are not always associated with both warming breakup and warming freeze-up date trends 9). This suggests that not just the length of time but also which phenomena are being investigated is important for understanding the context of the trends that are documented (e.g. Wynne, 2000), as there are numerous examples where warming trends in either the breakup or the freeze-up date are not matched with a correlative increase in the number of open-water days. Thus, environmental changes inferred from only breakup or freeze-up data are potentially misleading as they might not reflect wider changes in the ice season length, and they should be interpreted cautiously.

Conclusions
Utilising a number of different datasets, a series of analyses have been used to investigate how the number of annual open-water days per year and the timing of breakup and freeze-up dates have changed for water bodies that ephemerally freeze across the Northern Hemisphere. Five overlapping time periods (1931-1960, 1946-1975, 1961-1990, 1976-2005, and 1931-2005) have been investigated across 678 sites with data in at least one of the time periods to provide 3510 time series of lake and river ice phenology change to be statistically, spatially, and temporally analysed. A warming signal has been observed that shows the number of annual open-water days increased by 0.63 d per decade across the Northern Hemisphere from 1931-2005. The breakup trends display a strong correlation with temperature observations in the weeks preceding breakup and during winter ice growth, suggesting that temperature can be confidently used to explain a large proportion of variability. Freeze-up trends show the greatest variability and are less easily predicted from air temperature changes compared to breakup trends. This is likely because freeze-up is not guaranteed to occur simply because temperatures have moved below 0 • C as water kinetics can prevent freeze-up. When the time series are investigated on smaller timescales to explore temporal changes, trends for the number of open-waters days show variation, with the two most recent 30-year time periods displaying a consistent trajectory towards more open-water days that is nonlinear with respect to magnitude. In general, the number of open-water days closely resembles breakup patterns, suggesting that breakup trends are the main driver in open-waterday trends. Four key conclusions have been drawn from this research: (1) an accelerating warming signal is clearly observable in breakup dates at many sites and is reasonably well aligned to broad regional temperature trends, (2) freezeup trends are more spatiotemporally complex and display weaker temperature correlations, (3) the length of the openwater season has generally increased through time and was predominantly driven by earlier breakup dates, and (4) care needs to be taken when interpreting the implications of ice phenology changes at sites that only have breakup or freezeup data. These results highlight the need for a more detailed understanding of historical changes and their causes to fully unravel the potential implications of ice phenology when projecting future climate changes.
Data availability. All of the raw data are available through the National Snow and Ice Data Center (https://nsidc.org/, last access: 13 July 2020) or by contacting the relevant meteorological institutes.
Author contributions. AMWN led the project analysis, writing, figure preparation, and revisions with input on all from DJM.