1. Introduction
Changes in the amount of sea-ice coverage of the Baltic Sea in the last millennium were related to significant climate changes within that period. In the Baltic Sea Basin, the Medieval Warm Period (MWP) lasted from the beginning of 10th to the first half of 14th century. It was followed by a gradual cooling for ∼200 years (the Transitional Period), which initiated the Little Ice Age (LIA) that predominated until the mid-19th century (Brázdil and others, Reference Brázdil, Pfister, Wanner, Von Storch and Luterbacher2005; Niedźwiedź and others, Reference Niedźwiedź2015). The mean winter air temperature in the Baltic Sea oscillated ∼4°C, i.e. ∼2°C colder than present conditions (Eriksson and others, Reference Eriksson, Omsted, Overland, Percival and Mofjeld2007). Similar differences in winter temperature between the LIA and the modern period have been obtained for Poland by Przybylak and others (Reference Przybylak2005) and Przybylak (Reference Przybylak2016), and for Central and Northern Europe (Stockholm) by Büntgen and others (Reference Büntgen2011). Mean thermal conditions in the LIA showed some variability, so Eriksson and others (Reference Eriksson, Omsted, Overland, Percival and Mofjeld2007) distinguished eight cooler and warmer periods of various durations. In extremely cold winters, the mean air temperature fell below −8°C. By contrast, summers were both warmer and cooler in the LIA than now. In northern Scandinavia, they were ∼1–2°C cooler than now only for the first hundred years of the LIA and in the 19th century (Büntgen and others, Reference Büntgen2011). In Europe, the average annual air temperatures were ∼1.5°C lower in the LIA than now, and even 2°C lower in extremely cold years (PAGES 2k Consortium, 2013, Supplementary information, 2013, see Figs S1 and S2). The prevailingly cold climate of this period led to crop failures and even famine in Europe, especially when the growing season was cooler than today (Abel, Reference Abel2013). In addition to these economic and social consequences, it also quite often increased the extent and thickness of sea ice in the Baltic Sea, as reflected in numerous historical sources (Speerschneider, Reference Speerschneider1915). During severe winters in the LIA, the Baltic Sea froze to such a degree that goods could be transported by sled over the ice between Hanseatic ports, and even inns could be built on the ice to serve this practice (e.g. Magnus, Reference Magnus1555; Fant, Reference Fant EM1818; Korner Reference Korner and Schwalm1895; Olaus, Reference Olaus1772). Sled journeys on the ice usually happened in the near-shore zone, especially on the southern Baltic coast, as they offered a safer alternative to trips on snow-covered land routes (North, Reference North2011).
Over longer periods, ice extent on the Baltic Sea is closely related to climate change. Changes in the maximum ice extent (MIE) in the Baltic Sea are a good indicator of the severity of winter in this area. The average air temperature for winter months in northern Europe significantly influences the formation of ice (Seinä and Palosuo, Reference Seinä and Palosuo1996; Omstedt and others, Reference Omstedt, Pettersen, Rodhe and Winsor2004). As shown by reconstructions of MIE in the Baltic since 1720 (Seinä and Palosuo, Reference Seinä and Palosuo1996 and https://www.ilmatieteenlaitos.fi/jaatilastot, 2021), the MIE of the Baltic Sea has been much smaller in recent decades than in the end of the LIA, which is undoubtedly related to the global warming we are witnessing (IPCC, 2021). In the Baltic Sea, changes in macroscale atmospheric circulation over the North Atlantic Ocean also significantly affected the variability of ice conditions (Koslowski and Glaser, Reference Koslowski and Glaser1999). The North Atlantic Oscillation (NAO) reflects the climatic variability of winters in the North Atlantic and Europe (Walker and Bliss, Reference Walker and Bliss1932). The greatest differences in atmospheric pressure between the Icelandic Low and the Azores High occur in boreal winter. The nature of this atmospheric circulation is described in the winter NAOdjfm index (Hurrell, Reference Hurrell1995; Cook and others, Reference Cook, Kushnir, Smerdon, Williams, Anchukaitis and Wahl2019).
The macroscale atmospheric circulation that strongly influences Baltic Sea weather conditions is controlled by the aforementioned permanent baric regimes and the seasonally developing Siberian High in winter and Asian Low in summer. Climatic conditions, including ice conditions in the Baltic Sea region, are shaped by four weather regimes. Two are related to the positive and negative phases of the NAO, while the third is a ‘blocking’ regime (an anticyclonic ‘ridge’ over Scandinavia) and the fourth is the Atlantic (high) regime to the west of continental Europe (Hurrell and Deser, Reference Hurrell and Deser2010; Lehmann and others, Reference Lehmann, Getzlaff and Harlaß2011).
Numerical models of sea ice can be used to simulate the degree of the sea-ice coverage of sea bodies under various climatic conditions. To estimate ice occurrence, one-dimensional models (Leppäranta, Reference Leppäranta1993; Yu and others, Reference Yu, Leppäranta, Zhijun, Cheng, Mengxi and Demchev2015) based on air temperature as the dominant factor in sea-ice formation were used. In these models, negative temperatures determined the formation and growth of ice. In recent decades, the Baltic Sea ice has been in many numerical studies, starting in the 1980s (Leppäranta, Reference Leppäranta1981; Omstedt, Reference Omstedt1990; Leppäranta and Zhang, Reference Leppäranta and Zhang1992; Haapala and Leppäranta, Reference Haapala and Leppäranta1996; Omstedt and Nyberg, Reference Omstedt and Nyberg1996; Haapala, Reference Haapala2000; Cheng and others, Reference Cheng, Vihma, Pirazzini and Granskog2006). In subsequent years, increasingly advanced models have enabled the simulation of Baltic Sea ice conditions, including CRO (Meier and others, Reference Meier, Döscher, Coward, Nycander and Döös1999), RCAO (Döscher and others, Reference Döscher, Willén, Jones, Rutgersson, Meier, Hansson and Graham2002), HAMSOM (Schrum and others, Reference Schrum, Hübner, Jacob and Podzum2003), RCA4-NEMO (Dieterich and others, Reference Dieterich, Schimanke, Wang, Väli, Liu, Hordoir, Axell, Höglund and Meier2013) and HIROMB (Löptien and Axell, Reference Löptien and Axell2014). Sea-ice models coupled with numerical models based on the Princeton Ocean Model (POM) have been used to analyse Baltic Sea ice conditions by Herman and others (Reference Herman, Jędrasik and Kowalewski2011) and Ryabchenko and others (Reference Ryabchenko, Dvornikov, Haapala and Myrberg2010) Recently, Pemberton and others (Reference Pemberton, Löptien, Hordoir, Höglund, Schimanke, Axell and Haapala2017) developed the new NEMO–LIM3.6-based ocean-sea-ice setup for the North Sea and Baltic Sea region for climate studies and forecasts of various ice parameters, i.e. concentration, extent and thickness of sea ice. Baltic Sea ice conditions have also been modelled with the help of coupled ice–ocean numerical models like CEMBS (Janecki and others, Reference Janecki, Nowicki, Kańska, Golenko and Dzierzbicka-Głowacka2018) and CESM (Jakacki and Meler, Reference Jakacki and Meler2019). The European Centre for Medium-Range Weather Forecasts (ECMWF) developed the Integrated Forecasting System (IFS), the first global coupled ice–ocean–atmosphere medium-range forecasting system (Keeley and Mogensen, Reference Keeley and Mogensen2018), which also covers the Baltic Sea.
Many studies have attempted to perform long-term simulations of various ice characteristics for the period of recent climate change Haapala and others, Reference Haapala, Meier and Rinne2001; Meier, Reference Meier2002; Meier and others, Reference Meier, Döscher and Halkka2004; Meier, Reference Meier2006; Luomaranta and others, Reference Luomaranta, Ruosteenoja, Jylhä, Gregow, Haapala and Laaksonen2014; Höglund and others, Reference Höglund, Pemberton, Hordoir and Schimanke2017). However, to date, few have covered earlier historical periods. Meier and Kauker (Reference Meier and Kauker2002) simulated MIE in the Baltic for the period 1902–98 using the Rossby Centre coupled ice–ocean model (RCO). Hansson and Omstedt (Reference Hansson and Omstedt2008) modelled ice conditions in the Baltic Sea for the period 1500–2001 using the PROBE-Baltic model. The results showed a close relationship between maximum extent of sea ice in the Baltic and average air temperature in the December to February period. Schimanke and others (Reference Schimanke, Meier, Kjellstrӧm, Strandberg and Hordoir2012) carried out model simulations of MIE in the Baltic region in the MWP and LIA based on the RCO. LIA was the coldest period in the history of the last millennium, so it is important to know the ice conditions of the Baltic Sea and the air temperature prevailing at that time in order to obtain, firstly, knowledge about the scope of their natural changes, and secondly, to determine the changes that have occurred in relation to modern times.
Knowledge about the climatic conditions (including the most important factor—air temperature, which significantly affects the course of all processes in both marine and terrestrial ecosystems) is fundamental information concerning the shaping of living conditions. Therefore, it is a pertinent question whether the previously estimated average air temperatures prevailing during this period can be verified using models based on historical data on ice cover. Our main research objective was to assess variability in air temperature and ice conditions in the Baltic Sea in the final stage of the LIA (1721–1860). This was achieved using a numerical model of sea ice and information from historical sources.
Model simulations employing an appropriately validated PM3D hydrodynamic model coupled with the sea-ice model were used to recreate and verify spatial changes in ice conditions. These simulations were performed for the four scenarios presented below. These scenarios differed in average annual air temperature. The reliability of analyses that use a model depends on the model type (which should be suited to the assumed goals), on an appropriate set of empirical data being used to calibrate it, and on the correct definition of the boundary conditions. For the analysis of parameters relating to ice conditions, the obtaining of reliable results depends on whether the model correctly simulates the freezing of a water body under various conditions, both during severe and mild winters. Such a model must therefore be validated against long-term observational data. In our study, long-term data not only from the modern period but also sea-ice reconstructions for the Baltic Sea and extensive knowledge of air temperature values from the second half of the LIA, which is our study period, were used to validate the model.
The modelled Baltic Sea ice characteristics were compared against data from other sources. This resulted in a reliable tool allowing spatial distributions of air temperature or Baltic Sea-ice cover to be analysed for historical periods similarly as for modern times.
2. Location, data and methods
2.1. Location
The analyses covered the Baltic Sea area, including the Kattegat (Fig. 1). The Baltic Sea is a shallow inland sea that connects to the North Sea via the Danish Straits. Including Kattegat, it covers 415 × 103 km2 (HELCOM, 2007). The depth of the Baltic Sea averages 54 m and is greatest at Landsort Deep (459 m) (Leppäranta and Myrberg, Reference Leppäranta and Myrberg2009). The large inflow of fresh water, which reaches 6670 m3 s−1 (HELCOM, 2011), combined with the restricted water exchange with the North Sea, cause the Baltic Sea to have low salinity. Salinity ranges from 2 to 25 psu (Leppäranta and Myrberg, Reference Leppäranta and Myrberg2009) and averages 7.36 psu (Meier and Kauker, Reference Meier and Kauker2003), which, combined with its shallow depth, increases the susceptibility to freezing. The maximum proportion of the Baltic Sea area to be occupied by ice varies from 12.5% to 100% by year (Leppäranta and Myrberg, Reference Leppäranta and Myrberg2009). In Bothnian Bay, the annual maximum ice thickness is usually 0.65–0.80 m, reaching 0.3–0.5 m even in mild winters. On the Polish and German coasts of the southern Baltic Sea, the annual maximum ice thickness rarely exceeds 0.5 m (Vihma and Haapala, Reference Vihma and Haapala2009).
Modelled area, bathymetry and location of stations used in analysis and model validation.
2.2. MIE and ice thickness
To analyse historical ice conditions in the Baltic Sea and to validate the model for the period 1958–2007, a 300 year series of MIEs in the Baltic published by Seinä and Palosuo (Reference Seinä and Palosuo1996) was used. This series was supplemented for subsequent years by the Finnish Meteorological Institute (FMI) (Fig. 2). For the 18th and 19th centuries, this series was reconstructed mainly using records of ice and weather conditions in the Baltic Sea taken from Finnish lighthouse logs, the Danish lighthouse in Skagen, Finnish newspapers and travel journals. Unfortunately, because they were closed in winter, the Finnish lighthouses didn’t provide information about mid-winter ice cover. It was not until the second half of the 19th century that regular observations of ice conditions along the Finnish coast and then in offshore areas began; these used information from ships and icebreakers, and later from aerial reconnaissances. It is worth mentioning that average air temperatures (December–March) from Helsinki, Stockholm and Mariehamn were used in reconstructing the MIE. Coastal data were also used to estimate ice conditions in the open sea (Seinä and Palosuo, Reference Seinä and Palosuo1996). According to Vihma and Haapala (Reference Vihma and Haapala2009), the greatest inaccuracies in estimating MIE in the first period of ice reconstruction concern severe winters, such as 1739/40, but there is also a relative uncertainty of 15–25% for mild ice winters in the 1760s. Currently, MIE is largely estimated using synthetic aperture radar imagery (Mäkynen and others, Reference Mäkynen, Karvonen, Cheng, Hiltunen and Eriksson2020). The MIE in the Baltic averaged 241 × 103 km2 in 1761–1860 but only 186.5 × 103 km2 in 1958–2007, while the respective medians were 245 × 103 and 190 × 103 km2.
Variability in maximum ice extent in the Baltic Sea in the years 1720–2020 (blue line) and 50 year moving average (black line); chart shows periods of model analysis (1721–60, 1761–1860) and model validation (1958–2007) (maximum extent of ice cover in the Baltic Sea | European Environment Agency’s home page).
To validate the sea-ice model, ice-thickness measurements from the Swedish Meteorological and Hydrological Institute (SMHI) stations at Jarnashamn, Ratan, Kalix and Furuogrund were used, along with weekly ice-thickness measurements from the FMI stations at Helsinki, Airisto, Virpiniemi, Hailuoto, Orrengrund, Utö, Susiluoto and Kumlinge (Fig. 1).
2.3. Documentary evidence
Various historical sources were used to reconstruct winter air temperatures in areas adjacent to the southern Baltic Sea in the period 1721–78. The fundamental ones were those that presented long-term, systematic observations. The longest such series on the Baltic Sea coast is the Gdańsk observations of Gottfried Reyger that extend from 1721 to 1786 (Reyger, Reference Reyger1770, Reference Reyger1788). He indicated various weather phenomena, including air temperature and precipitation. He usually summarised weather conditions for each month. He also used some instrumental data (for details, see Filipiak and others, Reference Filipiak, Przybylak and Oliński2019). This observation series was supplemented with data from other subdaily air temperature measurements in (i) Gdańsk (1739–72: Hanov, Reference Hanov1772), (ii) Wrocław and Lower Silesia (1717–26: Kanold, Reference Kanold1725; and 1727–30:; Büchner, Reference Büchner1731, Reference Büchner1732, Reference Büchner1733, Reference Büchner1734), (iii) Toruń (1760–64: Pospieszyńska and Przybylak, Reference Pospieszyńska, Przybylak and Bednorz2010) and Warsaw (1760–62: Rojecki, Reference Rojecki1965). In addition, for Toruń and surrounding area, weather descriptions by J. Richtsteig for 1704–30 (Richtsteig, Reference Richtsteig1704–Reference Richtsteig1730), D. Brauer for 1719–50 (Brauer, Reference Brauer1719–Reference Brauer1750) and J.H. Zernecke for 1700–25 (Zernecke, Reference Zernecke1727) were used, among many others. Weather notes from 18th-century journals and various administrative sources were also used in the analysis. Good use was also made of the recently published book by Girguś (Reference Girguś2022).
Weather records for winters taken from historical sources were indexed on the 7-point scale (extremely cold: −3, very cold: −2, cold: −1, normal: 0, warm: +1, very warm: +2, extremely warm: +3) proposed by Pfister and others (Reference Pfister, Kington, Kleinlogel, Schüle, Siffert, Frenzel, Pfister and Gläser1994). Winters of indices −3 and −2 were usually characterised by long freezes (∼2–3 months) resulting in the appearance of ice in the Baltic Sea. Next, the indexation of winters was also reduced to 5-point and 3-point scales to facilitate comparison against the classifications of winters and ice conditions proposed by Seinä and Palosuo (Reference Seinä and Palosuo1996) and Koslowski and Glaser (Reference Koslowski and Glaser1995), respectively.
For the period from 1779 to 1860, winters were assessed based on a series of homogenised monthly mean air temperatures for Warsaw, which is the longest series available for Poland (Lorenc, Reference Lorenc2000). It was reported that, especially in winter, temperature values are highly correlated across the entirety of Poland (Kożuchowski and Żmudzka, Reference Kożuchowski and Żmudzka2003; Przybylak and others, Reference Przybylak, Pospieszyńska, Wyszyński and Nowakowski2014), and that the averages for the entire country also correlate highly with the averages in most areas of Europe (including the entire Baltic Sea) (Luterbacher and others, Reference Luterbacher, Przybylak, Majorowicz, R and Kejna2010). The severity of winters was assessed according to the ranges of average winter air temperature used in works by Przybylak and others (Reference Przybylak, Wójcik, Marciniak, Chorążyczewski, Nowosad, Oliński and Syta2004, Reference Przybylak2005) and expressed there using the winter evaluation indices according to the aforementioned 7-point scale [+3, +2, +1, 0, −1, −2, −3]. For the purposes of this study, the three warm categories (warm [+1], very warm [+2] and extremely warm [+3]) and the three cold categories (cold [−1], very cold [−2] and extremely cold [−3]) of winters were reduced to the two categories also used in Seinä and Palosuo (Reference Seinä and Palosuo1996). The modified ranges of average winter air temperature in Warsaw corresponding to the five categories of winter severity are given below:
To assess the nature of winters according to the three-grade classification applied by Koslowski and Glaser (Reference Koslowski and Glaser1995), the five categories of winters were first distinguished as described above, and then ‘extremely mild + mild’ as well as ‘extremely severe + severe’ winters were combined and named ‘weak’ and ‘severe’, respectively.
2.4. Sea-ice model coupled to PM3D hydrodynamic model
Baltic Sea ice conditions were modelled using a sea-ice thermodynamics and dynamics model (Herman and others, Reference Herman, Jędrasik and Kowalewski2011) coupled with a 3-D PM3D hydrodynamic model (Kowalewski and Kowalewska-Kalkowska, Reference Kowalewski and Kowalewska-Kalkowska2017). The PM3D model is an enhanced version of the M3D model (Kowalewski, Reference Kowalewski1997) and is based on the POM model developed by Blumberg and Mellor (Reference Blumberg, Mellor and Heaps1987). Water exchange with the North Sea through the Skagerrak strait across the open border between Skagen and Gothenburg was taken into account. The so-called radiation boundary condition for open border flows used hourly observed sea levels in Gothenburg and long-term average monthly salinity values for the area. Calculations were made in a rectangular numerical grid with a resolution of 3 nautical miles (∼5.5 km). The vertical division used the sigma-transformation and included 18 layers, regardless of depth. Atmospheric forcing at the sea surface was set based on the results of reanalysis of the REMO model (REgional MOdel; Jacob and Podzun, Reference Jacob and Podzun1997; Feser and others, Reference Feser, Weisse and von Storch2001) with a 1-hour step in the period 1958–2007. The calculations take into account the average monthly discharges of the 167 largest rivers. More details on the model setup can be found in the work of Jędrasik and Kowalewski (Reference Jędrasik and Kowalewski2019).
A detailed description of the sea-ice thermodynamics and dynamics model and the methods used to couple it with the M3D/PM3D hydrodynamic model are given in detail in the work of Herman and others (Reference Herman, Jędrasik and Kowalewski2011). Meanwhile, the simulations carried out here employed data from a different atmospheric model; the coefficients parameterising heat fluxes at ice–water–atmosphere borders in the thermodynamic part of the model thus needed to be recalibrated (see Appendix). The calibration was based on an analysis of discrepancies between values of ice extent, ice thickness and the temperature and salinity of water obtained from a series of model simulations for the period 1958–2007 on the one hand, and, on the other, empirical values for the same parameters obtained from International Council for the Exploration of the Sea (ICES https://data.ices.dk/). As a result, new values of empirical coefficients defining the heat fluxes were established (Table A1). An independent set of observational data was used to assess the accuracy of the values of ice characteristics estimated by the calibrated model. The results of this validation are presented in Section 3.
For numerical model results, the MIE can be calculated in different ways. For example, Pemberton and others (Reference Pemberton, Löptien, Hordoir, Höglund, Schimanke, Axell and Haapala2017) assumed the area for which the concentration of ice was at least 15%. For the present work, the extent of ice for each day was determined as the actual ice area, which takes into account the concentration (A), understood as the ratio of ice-covered area to total area of the numerical grid. MIE was calculated as the sum of the areas of individual grids multiplied by the ice concentration:
\begin{equation}MIE = \sum\limits_i {(A_i} \cdot {P_i)} \end{equation} where 
${A_i}$ is the ice concentration in the ith computational grid node and 
${P_i}$ is the area of the ith computational grid node.
2.5. Cumulative ice thickness
Because ice thicknesses and concentrations vary throughout the season, cumulative ice thickness (CIT) can be used to facilitate ice analysis. This characteristic is calculated as the sum of the daily products of ice thickness and ice concentration for all days since the beginning of the ice season. This quantity was originally called the ‘accumulated areal ice volume’ (Koslowski and Loewe, Reference Koslowski and Loewe1994). It is a measure analogous to the accumulated freezing degree days used in climatology and meteorology. It is expressed in metres multiplied by day (m · d), so it in fact expresses the thickness of ice accumulated during the winter season, taking into account its concentration. For a single observation station where ice thickness and concentration measurements were taken CIT can be expressed as follows:
\begin{equation}CIT_k{\text{ }} = \mathop {\mathop \sum}\limits_j {(Ah)_{jk}}\end{equation}where A is the ice concentration, h is the ice thickness (m), j is the index of winter season days and k is the index of observation station.
For several measuring stations for a larger area, such as a bay or even for the entire sea, an averaged CIT can be determined:
\begin{equation}CIT{\text{ }} = \frac{1}{n} \sum\limits_k \sum\limits_j (Ah)_{jk}\end{equation}where n is the number of stations.
This allows for the analysis of ice accumulation differences between various regions or stations.
The model results can be used to calculate the CIT for a certain area, e.g. for the entire Baltic Sea or for each individual node of the computational grid:
\begin{equation}{CIT_i} = \sum\limits_j {(Ah)_{ij}}\end{equation}where i is the number of computational grid node.
This allows one easily to characterise the amount of sea-ice coverage in spatial terms on a map. The maps can be plotted for individual years or to present averages for longer periods. The CIT maps presented in Figure 3 show the spatial differentiation of Baltic Sea ice coverage for a very severe winter (1987) and an exceptionally mild winter (1989). In this work, the year of a winter is identified by the year in which it ended: thus, ‘the winter of 1987’ is the winter of 1986/87. During the winter of 1987, CIT >10 m · d covered the area north of Öland and Gotland. In other water bodies, such values occurred only near the shore and values in the middle of the Baltic Sea and Kattegat decreased to ∼1–2 m · d (Fig. 3a). During the very mild winter of 1989, CIT exceeded 10 m · d only in the inner waters of the Gulf of Finland and Bothnian Bay, in shallow lagoons off the coast of the southern Baltic Sea, and along the coasts of Estonia, Finland and Sweden (Fig. 3b). A CIT value of 10 m · d can therefore be interpreted as compact ice covering 100% of the surface with a thickness of 1 m occurs in a given location for 10 days, or compact ice with a thickness of 0.5 m occurs for 20 days in a given season.
Cumulative ice thicknesses in the Baltic Sea in the years (a) 1987 and (b) 1989.
The CIT values calculated based on the modelled ice-cover concentrations and thicknesses illustrate the spatial distribution of ice conditions in the entire Baltic Sea effectively, as well as their changes over time. Using spatially averaged CIT values, the severity of winters can also be assessed for an entire water body or any part thereof. The averaged CIT values for the entire Baltic area ranged from 2.2 to 30.9 m · d in the period 1958–2007 (Fig. 4). These values correlated strongly with maximum ice area in the Baltic Sea (r = 0.94). High values of both parameters corresponded to severe winters, while the lowest values corresponded with very mild winters. CIT was highest at 30.9 m · d in the winter of 1987, while it was lowest at 2.2 m · d in the winter of 1989. The corresponding observed ice areas amounted to 407 × 103 and 60 × 103 km2 (almost seven times less). This sequence of extremely severe and mild winters within a 2 year period indicates the significant dynamics of ice phenomena in the Baltic Sea.
Variability of the cumulative ice thickness averaged over the entire Baltic Sea (red line) and maximum ice extent (blue dots) in the Baltic Sea, 1958–2007.
In analysing ice in the Western Baltic, Koslowski and Glaser (Reference Koslowski and Glaser1995) used a reconstruction of the severity of winters for 1701–1993 based on CIT on the German and Danish coasts in the period 1897–1993. For the earlier years (1701–1860), the reconstruction is based on historical records on ice cover duration in the Danish waters of the Baltic Sea, notes on ice cover along the German coast, and archival descriptions of ice conditions on the River Trave between Lübeck and Travemünde in 1701–1848 and on the River Elbe and other water bodies near Hamburg. Public records from various archival sources were also used.
3. Results and discussion
3.1. Model validation for 1958–2007
Previous validations of the results of long-term simulations by the PM3D model presented in Jędrasik and Kowalewski (Reference Jędrasik and Kowalewski2019) have shown the model’s values of sea water salinity, temperature and level to be very close to measured values. However, these simulations were performed without the thermodynamics and sea-ice dynamics module. Thus, the validation of the PM3D model that was performed in this study regarded the modelled MIEs and thicknesses in the Baltic Sea. First, ice extents for each day were calculated by summing the area of the grids occupied by ice multiplied by the ice concentration of each grid node, and then the maximum value for the ice season was determined. The model results were compared against a reconstruction of the maximum extent in the period 1958–95 (Seinä and Palosuo, Reference Seinä and Palosuo1996) supplemented by FMI observations (https://en.ilmatieteenlaitos.fi/baltic-sea-ice-winters) for 1996–2007 (Fig. 5). The model simulations were consistent with observations of the MIE. This agreement was slightly less for the first three decades. However, this may be related to the lower accuracy of ice extent observations before the use of satellite techniques.
Maximum ice extent in the Baltic Sea for 1958–2007: model simulation compared to observation (FMI) and severity of winters according to Seinä and Palosuo (Reference Seinä and Palosuo1996).
The high agreement between the model results and observations is also shown in a scatter plot (Fig. 6). The correlation coefficient was 0.97, and RMSE was 20.4 × 103 km2. The model slightly underestimated the MIEs in the Baltic Sea (BIAS = −13 × 103 km2), which is associated with the aforementioned tendency for the first decades of the simulation.
Maximum annual ice extents in the Baltic Sea in 1958–2007: model values (PM3D) compared against observations (FMI).
The comparison of the MIE modelled with the model used in this study and other models of the Baltic Sea shows the PM3D model to be highly accurate. The validation of the RCO model (Meier and Kauker, Reference Meier and Kauker2002) showed that the correlation coefficients for various simulations for 1903–98 ranged between 0.85 and 0.89 with RMSEs between 55.6 × 103 and 74.7 × 103 km2. In the case of the PROBE-Baltic model (Hansson and Omstedt, Reference Hansson and Omstedt2008), a comparison of the modelled against observed maximum Baltic Sea ice extents for the period 1895–1999 showed that the model overestimated the ice extent by 5 × 103 km2, the correlation coefficient between simulations and observations was 0.9, and modelling accuracy was found to be greater for severe winters than for mild ones. On the other hand, RCO simulations by Schimanke and others (Reference Schimanke, Meier, Kjellstrӧm, Strandberg and Hordoir2012) showed that the RCO model underestimated the value of the annual MIE in the years 1961–2000 by over 20 × 103 km2 compared to data published by Seinä and Palosuo (Reference Seinä and Palosuo1993). Simulations of ice conditions with the NEMO-LIM3.6 model (Pemberton and others, Reference Pemberton, Löptien, Hordoir, Höglund, Schimanke, Axell and Haapala2017) in the similar period 1961–2006 were shown to be good reflections of sea-ice extent (understood as the area of sea-ice concentration of at least 15%) with the exception of Kattegat, where the model exceeded observed values. The correlation coefficient between observed and model values was 0.85. A comparison of winter types for this period (according to the work of Seinä and Palosuo, Reference Seinä and Palosuo1996) showed that this model tended to overestimate the number of average annual maximum Baltic Sea ice extents and underestimate the number of extreme cases.
To validate the model, its estimated ice thicknesses were also compared against observations at selected measuring stations mainly along the Bay of Bothnia and the Gulf of Finland (Fig. 1). Good agreements were obtained for most of the FMI measuring stations (Fig. 7). However, for SMH stations the compliance was lower. At the Jarnashamni and Ratan stations, the compliance was unsatisfactory, as shown in the data scatter plots (Fig. 7a and b).
Comparison of modelled and observed ice thicknesses for four SMHI stations: (a) Jarnashamn, (b) Ratan, (c) Kalix, (d) Furuogrund, and four FMI stations: (e) Helsinki, (f) Virpiniemi, (g) Hailuoto, (h) Susiluoto.
For all stations for which ice-thickness measurements were available, statistical parameters were calculated to characterise the model’s compliance with observation (Table 1). At the SMHI stations (Jarnashamn, Ratan, Kalix and Furuogrund), the errors (RMSE from 13 to 16.2 cm) and negative biases were larger, indicating that ice thicknesses were underestimated by the model. By contrast, at the stations where measurements were made by the FMI (Helsinki, Airisto, Virpiniemi, Hailuoto, Orrengrund, Utö, Susiluoto, Kumlinge), RMSEs were smaller (from 6.7 to 12.1 cm), while biases were both negative and positive. Correlation coefficients were relatively high, ranging from 0.8 to 0.94.
Statistical characteristics of the agreement between modelled and observed ice thicknesses for various measuring stations.
The model was found to agree much better with observations for FMI stations than for SMHI stations (the lowest compliance for Jarnashamn and Ratan), which may be related to methodological differences in how each research centre conducts measurements. This may also be related to the location of the stations. The SMHI stations are located along the western coast of the Gulf of Bothnia, while the FMI stations are located on its eastern side and along the northern coast of the Gulf of Finland, where the wind conditions are different. Model accuracy can be affected when the rheology of sea ice is modelled under different wind directions, such as when specific winds prevail that cause sea ice to pile up. However, it should be remembered that the resolution of the model (∼5.5 km) often does not allow for an accurate reproduction of prevailing conditions near the land, where ice-thickness measurements are usually taken. Differences between modelled and observed ice thicknesses may also result from differences in spatial scale. Field measurements determine ice thickness over a small scale of a metre, while the model uses a much larger scale of several kilometres. The PM3D model validation results for ice thickness compare favourably against the rare reports of validations of point ice thicknesses for other models in the Baltic Sea (Saloranta, Reference Saloranta2000; Dieterich and others, Reference Dieterich, Schimanke, Wang, Väli, Liu, Hordoir, Axell, Höglund and Meier2013; Löptien and others, Reference Löptien, Mårtensson, Meier and Höglund2013; Pemberton and others, Reference Pemberton, Löptien, Hordoir, Höglund, Schimanke, Axell and Haapala2017). Simulations of ice thickness at Swedish stations using the RCO-HELMI (Löptien and others, Reference Löptien, Mårtensson, Meier and Höglund2013) and NEMO-LIM3.6 (Pemberton and others, Reference Pemberton, Löptien, Hordoir, Höglund, Schimanke, Axell and Haapala2017) models showed a similar underestimation compared to observation. This was especially true for the Ratan station, for which the PM3D model, too, gave its greatest underestimations of ice thickness (BIAS = −10.55 cm) and showed its largest random errors (RMSE = 16.33 cm). Perhaps this is related to the specific local conditions at this station, which make the modelling of ice thickness difficult. It is worth mentioning that the NEMO-Nordic model (Pemberton and others, Reference Pemberton, Löptien, Hordoir, Höglund, Schimanke, Axell and Haapala2017) heavily underestimated the ice thickness in the range 1.2–2.8 m, whereas the ice thickness in the ranges 0.4–0.5 m and 3.0–4.0 m was overestimated. The validation of the PM3D model carried out in this work showed a good agreement between the simulated and observed ice extents and thicknesses in the Baltic Sea. It was carried out for a long-term, 50 year simulation of sea-ice dynamics, which confirmed the usefulness of the model for the analysis of historical ice cover.
3.2. The climate in the Baltic Sea region in the final phase of the LIA on the basis of model simulations and documentary evidence
During the LIA period, climate changes were associated with air circulation over the North Atlantic and Europe, and probably concerned various parameters. However, the lowering of the air temperature was the most important factor influencing the degree of ice cover in the Baltic Sea. To estimate the climatic conditions in the final phase of the LIA, a series of model simulations was performed in which the air temperature was changed. The Baltic Sea has no areas of permanent ice and its annual formation is caused by negative air temperatures. Therefore, unlike in terrestrial areas, precipitation, atmospheric humidity and other meteorological factors have limited influence on maximum sea-ice extent. The Baltic Sea is a relatively shallow brackish sea, which is influenced by many factors affecting the ice formation process, but air temperature is the most important one. Aside from the air temperature, all other meteorological parameters were assumed to be the same as derived from the REMO model for 1958–2007 and used to validate the model. In each simulation, the air temperature was dropped by a constant value ΔT, respectively: 1°C, 2°C, 3°C and 4°C. This involved subtracting ΔT from the instantaneous air temperatures over the entire 50 year simulation period, both in the cold and warm seasons. This procedure caused the average air temperature at each node of the calculation grid to be correspondingly lower by ΔT. The range of assumed temperature reductions was consistent with the range of differences in decadal mean winter air temperatures in Stockholm (Moberg and others, Reference Moberg, Bergström, Ruiz Krigsman and Svanered2002; Moberg, Reference Moberg2021) between the LIA period and recent decades (Fig. 8). The model results at the air temperature lowered by 1°C, 2°C, 3°C and 4°C were designated as versions v1, v2, v3 and v4, respectively. Reference calculations based on real conditions for the period 1958–2007 were designated as version v0. For all simulations, the monthly river discharges into the Baltic Sea were assumed to be the same as for 1958–2007.
Time series of decadal average winter air temperature (Dec–Feb) in Stockholm in 1760–2020. Source: MOBERG, Reference Moberg2021.
Based on the simulation results, MIE for individual years were calculated, followed by statistics on the occurrence of winters of different severities. This was done using a classification based on MIE (Seinä and Palosuo, Reference Seinä and Palosuo1996) that divided winters into extremely mild (<81 × 103 km2), mild (81 × 103–139 × 103 km2), average (139 × 103–279 × 103 km2), severe (279 × 103–383 × 103 km2) and extremely severe (>383 × 103 km2). In addition, means, standard deviations and medians were determined for the MIEs presented in the work of Seinä and Palosuo (Reference Seinä and Palosuo1996) (Table 2). The table presents analogous statistics for the four PM3D model simulation variants (v1, v2, v3 and v4). The percentage shares of winters of different severities in the reference period v0 (1958–2007) were identical between the reconstruction (Seinä and Palosuo, Reference Seinä and Palosuo1996, supplemented by FMI observations) and the model. The standard deviations were very similar, but the means and medians of the MIEs produced by the model were lower than those observed. Of the simulation variants, the statistical characteristics of the MIEs reconstructed for the period 1761–1860 were closest to v2, i.e. the simulation for which air temperature had been lowered by 2°C. The means of the MIEs were very similar to one another, as were the medians. Differences in the classification of severity of winters were greater. The percentage share of severe winters was smaller for the v2 simulation than in the reconstruction for 1761–1860, while for extremely severe winters it was larger.
Percentage share of winters with various degrees of severity and other statistical characteristics calculated based on the reconstruction (Seinä and Palosuo, Reference Seinä and Palosuo1996, supplemented by FMI observations for 1996–2007) and for model simulations: reference (v0) and simulations with air temperature lowered by 1°C, 2°C, 3°C and 4°C (v1, v2, v3, v4).
The results of the reference simulation (v0) were used to calculate the contemporary distributions of average CIT for the Baltic area (Fig. 9). This map shows the variation in the amount of sea-ice coverage in the Baltic Sea. The mean CIT is close to zero in the middle of the Baltic Sea, where sea ice is almost never present. In bays and lagoons, average CIT values are usually above 10 m · d, and in the northern part of the Bay of Bothnia they are up to 100 m · d. Next, average CIT maps were determined for the four simulations (v1, v2, v3 and v4) in each of which the air temperature had been lowered by a fixed value of 1°C, 2°C, 3°C and 4°C, respectively (Fig. 10). Lowering the air temperature extends the sea-ice period and increases the ice thickness, thereby increasing the CIT value. The spatial distribution of average CIT does not change much, but the most marked increase is in the middle of the Baltic Proper and in Kattegat. If the air temperature is lowered by 2°C, the average CIT reaches values over 1 m · d there. These values mean that ice would occur frequently in these regions, and not—as now—exceptionally. For even greater reductions in air temperature, by 3°C and 4°C, the increase in CIT is more visible in the middle of the Baltic Sea than in Kattegat.
Map of mean cumulative ice thickness in the Baltic Sea, 1958–2007.
Distribution of mean cumulative ice thickness (CIT) in the Baltic Sea: modelled for four hypothetical variants of mean annual air temperature in the end of the LIA relative to the average for the years 1958–2007; temperature values were lowered by (a) 1°C, (b) 2°C, (c) 3°C and (d) 4°C.
The maps depicting differences in the mean CIT value over the 50 year simulation period (Fig. 11) illustrate the variation in the mean CIT value calculated for simulations with reduced air temperature (Fig. 10) compared to the reference simulation (Fig. 9). A clear spatial differentiation is observed, showing an increasing impact of reduced air temperature on the thickness and duration of ice cover moving northward. A 1°C reduction in air temperature delineates three distinct regions (Fig. 11a). The first region extends from the Danish Straits to the Baltic Proper, where CIT increases were minor, up to 2 m · d, except for the eastern coastal zone and the northern part of the Baltic Proper, where increases reached up to 4 m · d. The second region includes the Gulf of Riga, the Gulf of Finland, and the Bothnian Sea, where CIT increases ranged from 4 to 8 m · d. The third region covers the Bothnian Bay, where maximum CIT increases ranged from 8 to 12 m · d. For a 2°C reduction in air temperature (Fig. 11b), the spatial distribution of CIT increases remained similar, but the values were approximately twice as high. CIT increased by 4–8 m · d from the Danish Straits to the Baltic Proper, 8–18 m · d in the Gulf of Riga, the Gulf of Finland, and the Bothnian Sea, and 17–23 m · d in the Bothnian Bay. A further simulation with a 3°C air temperature reduction (Fig. 11c) showed an approximate doubling of the mean CIT values across the entire Baltic Sea compared to the previous simulation (Fig. 11b), except for the Bothnian Bay, where CIT increased by ∼1.5 times, reaching 24–36 m · d. The results of the simulation with a 4°C air temperature reduction revealed another doubling of CIT values in the region from the Danish Straits to the Baltic Proper and an ∼1.5-fold increase in other areas, reaching 20–26 m · d in the Gulf of Riga, 22–34 m · d in the Gulf of Finland, 27–34 m · d in the Bothnian Sea, and 34–45 m · d in the Bothnian Bay.
Distribution of differences in mean cumulative ice thickness (CIT) in the Baltic Sea relative to the average for the years 1958–2007: modelled for four hypothetical variants of mean annual air temperature in the end of the LIA, temperature values were lowered by (a) 1°C, (b) 2°C, (c) 3°C and (d) 4°C.
Koslowski and Loewe (Reference Koslowski and Loewe1994) analysed amount of sea-ice coverage for the Western Baltic using observations from six stations in this region in the period 1879–1992. Then, based on the CITs presented in that work, Koslowski and Glaser (Reference Koslowski and Glaser1995) presented a reconstruction of the severity of winter for the period 1701–1993. An analogous analysis was made on the basis of PM3D model simulations. For each year of the two periods (1721–60 and 1761–1868), CIT was determined in the model grid nodes corresponding to the locations of these six stations. Then, using the classification of winter severity based on CIT, the percentage shares of winters of different severities were calculated for the reference simulation (v0) and for the simulations of lowered air temperatures v1, v2, v3, v4 (Table 3). Winters were determined by CIT to be weak (<2 m · d), moderate (2–6 m · d) and severe (>6 m · d). The classification based on reconstruction and model simulations is compared in Table 3. In the period 1721–60, the percentage shares of weak and moderate winters in the reconstruction were similar to those obtained from the v2 simulation, in which the air temperature had been lowered by 2°C (Table 3). However, the percentage share of severe winters (10%) was equal to that obtained for simulation v1, i.e. a decrease in air temperature of 1°C. In the second analysed period (1761–1860), the percentage share of weak, moderate and severe winters was closest to the reconstruction for simulations v3 and v4 (for which air temperature had been lowered by 3°C and 4°C).
Percentage shares of winters of various severities in the Western Baltic according to the reconstruction by Koslowski and Glaser (Reference Koslowski and Glaser1995) and according to PM3D v0 (reference simulation in 1958–2007) and PM3D v1, v2, v3, v4 (simulations of air temperature lowered by 1°C, 2°C, 3°C and 4°C). avg, var—mean value and variance of the ice winter index numerals.
For the winter severity classes they distinguished, Koslowski and Glaser (Reference Koslowski and Glaser1995) assigned ‘ice winter index numerals’ in the range from 0 to 3. The mean ice winter index numerals determined for the period 1721–60 and shown in Table 3 indicate a high convergence between the reconstruction (0.31) and the v2 model simulations (0.32). In the second period (1761–1860), the mean winter ice index (0.71) is between the values obtained from simulations v2 (0.58) and v3 (0.89).
As already mentioned in Section 2.4, the severity of winters in the period 1721–1860 (Tables 4 and 5) was reconstructed mainly for northern Poland (including the Baltic coast) based on weather descriptions (1701–78) and available average winter air temperatures from Gdańsk and Warsaw (documentary evidence). Based on the results obtained, the new classification presented in Tables 4 and 5 is worked out. Both the division of winters into five categories proposed by Seinä and Palosuo (Reference Seinä and Palosuo1996) and into three categories according to Koslowski and Glaser (Reference Koslowski and Glaser1995) were used.
Frequency of occurrence of various categories of winter in Northern Poland, 1721–60 and 1761–1860, reconstructed based on documentary evidence. Relates to the reconstruction of severities of winter according to the Seinä and Palosuo (Reference Seinä and Palosuo1996) classification.
Frequency of occurrence of various categories of winter in Northern Poland, 1721–60 and 1761–1860, reconstructed based on documentary evidence. Relates to the classification of winters according to Koslowski and Glaser (Reference Seinä and Palosuo1995).
In the period 1761–1860, the new classification of winters presented in this paper is significantly similar to the classification (presented in Table 2) of Seinä and Palosuo (Reference Seinä and Palosuo1996) in terms of frequency of extremely mild (8% new classification, 7% Seinä and Palosuo (Reference Seinä and Palosuo1996)), mild (17% in both cases) and extremely severe winters (12% vs 13%, respectively), while some discrepancies were noted for normal and severe winters. According to our reconstruction, there are more severe winters (44%) and fewer average winters (19%) than in the Seinä and Palosuo (Reference Seinä and Palosuo1996) classification (32% and 31%, respectively). The agreement between the two classifications is even greater for the period 1721–60. Slight differences can only be found for average winters (22.5% and 20%) and severe winters (27.5% and 30%) (see Tables 2 and 4). However, greater discrepancies were found with the classification of winters proposed by Koslowski and Glaser (Reference Koslowski and Glaser1995), regardless of research period (Table 4). High compliances were only obtained for the incidence of severe winters in 1761–1860 and moderate winters in 1721–60 (see Tables 3 and 4). On the other hand, the categories of weak/moderate (1761–1860) and weak/severe (1721–60) winters in our reconstruction significantly differed in frequency relative to the Koslowski and Glaser (Reference Koslowski and Glaser1995) winter classification. The reconstructed frequencies of the various categories of winters in Northern Poland and the southern Baltic coast correlate well with those calculated using the PM3D v2 model, and slightly less well with version 3 (compare Table 4 with Tables 2 and 3).
Schimanke and others (Reference Schimanke, Meier, Kjellstrӧm, Strandberg and Hordoir2012) conducted model simulations of the MIE in the Baltic for 50 year periods in the MWP and LIA. Those numerical experiments showed that the long-term variability within the observational data is larger than the difference between the modelled MWP and LIA when 50 year periods are considered for the observations and the simulations. The comparative analysis of the results of the four simulation variants in the PM3D model showed that it best recreates ice conditions and the severity of winters in the northern and central part of the Baltic in the LIA period when the mean annual air temperature is reduced by 2°C (v2) relative to observed values from the second half of the 20th century. A similar value for the LIA cooling in Europe relative to modern values was obtained by the PAGES 2k project (PAGES 2k Consortium, Supplementary information, 2013, see Figs S1 and S2 there) using available proxy data (mainly tree-rings) and, to a small extent, documentary evidence. This independent reconstruction also confirms that the choice of the PM3D model version assuming average annual air temperatures of 2°C below the present day better reflects thermal conditions during the LIA.
4. Summary and conclusion
The thermodynamic model of sea ice extended by the 3-D PM3D model showed very high agreement between modelled and observed distributions of the maximum extent of the Baltic Sea ice for the period 1958–2007. Slightly worse, but still satisfactory, was the agreement in the case of ice thickness. Taking into account the uncertainties associated with the observational data acquisition methodology, the validation allowed to conclude that the model effectively reconstructs the Baltic Sea ice cover. In particular, the agreement of the modelling results with the values recorded at the measuring stations confirms the reliability of the model-simulated spatial distributions of quantities characterising the sea ice. In the case of data from historical periods when there were no satellite observations, model simulations are the only way to obtain such information and make analyses that are difficult to carry out on the basis of point measurements, the number of which is limited both temporally and spatially.
Model simulations of spatial changes of the CIT across the Baltic Sea presented on maps can characterise average ice conditions better than can ice-extent maps, as they take into account both ice thickness, its concentration and the length of the ice period. It also allows CIT to be calculated for smaller areas, e.g. bays, and even for specific points. The simulations performed for air temperature reduced compared to the reference period (1958–2007) showed that each temperature reduction by 1°C increases CIT about twice in areas where sea ice occurs incidentally and ∼1.5 times in areas where ice appears annually (Bothnian Bay). The simulations show clear spatial differentiation, and allow to distinct the three regions with different impact of change in temperature on CIT: (1) quite minor response on area from the Danish Straits to the Baltic Proper, (with a slightly greater response to a temperature drop at the eastern coastal zone and the northern part of the Baltic Proper); (2) the Gulf of Riga, the Gulf of Finland, and the Bothnian Sea with more distinct response; and (3) the Bothnian Bay, where CIT increase is the sharpest.
The sea-ice simulation assuming an air temperature of 2°C lower than present was found to be highly similar to the reconstruction of the final phase of the LIA performed by Seinä and Palosuo (Reference Seinä and Palosuo1996) in terms of mean and median values of MIE. This result is consistent with the assessments of temperature anomalies during the LIA described in the Introduction. In the analysis conducted by the PAGES 2k team for the region of Europe, the average LIA air temperature reconstructed based on proxy data was ∼1.5°C lower compared to the second half of the 20th century (Pages 2k Consortium, 2013), while instrumental data available for Stockholm indicate even greater difference of value up to 2°C (Moberg and others, Reference Moberg, Bergström, Ruiz Krigsman and Svanered2002).
The model simulation of CITs in the western Baltic Sea for reduced air temperature values makes possible to use a winter severity classification method used by Koslowski and Glaser (Reference Koslowski and Glaser1995) and to compare the results of this classification against the reconstruction based on historical records from the LIA for this region of the Baltic Sea. The best compliance of statistical characteristics was obtained for simulations for which the air temperature had been lowered by 3°C and 4°C. These results indicate that climatic conditions in the LIA may have varied between regions of the Baltic Sea. In the western Baltic, temperatures were lower or winters were more severe than in other parts. This may explain the use of transport by ice sled, the construction of inns on the ice, military marches and other events described in historical sources. They attest to ice having occurred much more frequently and thickly during the LIA than today in this region.
The classification of winter severities according to five- and three-degree scales based on independent available documentary evidence from northern Poland (including the southern Baltic coast) and statistics on the frequency of occurrence of particular categories of winters calculated on their basis were shown to be greatly in line with analogous calculations for the entire Baltic Sea made based on an assessment of winter severities presented by Seinä and Palosuo (Reference Seinä and Palosuo1996) and much worse compliance with the Koslowski and Glaser classification (Koslowski and Glaser, Reference Koslowski and Glaser1995) for the western Baltic. These discrepancies are most likely related to the inaccuracy of the reconstruction of winter characteristics resulting from the paucity of available documentary evidence from the areas being compared.
Thanks to the use of a numerical thermodynamic model of sea ice, we confirmed the air temperature reconstructions in the Baltic region in the final LIA period made by other methods. The successful reconstruction of ice conditions in the LIA allows a PM3D model to be used to study the ice conditions of the Baltic Sea in other historical periods. The applied approach may open up new possibilities for the analysis of the spatial ice coverage of the Baltic Sea and the verification of hypotheses based on historical sources, which are mainly point-based.
Data availability statement
Daily ice thicknesses and concentrations calculated with the PM3D model in the period 1958–2007 and used in this work have been archived and are available under a CC BY-NC-SA license in the repository: 10.34808/cdd9-be22.
MIE was downloaded from https://www.eea.europa.eu/data-and-maps/daviz/maximum-extent-of-ice-cover-3 and https://en.ilmatieteenlaitos.fi/baltic-sea-ice-winters. Ice thicknesses from Jarnashamn, Ratan, Kalix, Furuogrund stations were obtained from the reports in A summary of the ice season and icebreaking activities in the winter seasons 1970/71–2006/07 published by the Swedish Maritime Administration, Norrköping and SMHI, Norrköping, and from Helsinki, Airisto, Virpiniemi, Hailuoto, Orrengrund, Utö, Susiluoto, Kumlinge of the Helsinki FMI.
Air temperature in Stockholm for 1760–2020 was obtained from: Anders Moberg (Reference Moberg2021) Stockholm Historical Weather Observations—Monthly mean air temperatures since 1756. Dataset version 3. Bolin Centre Database. 10.17043/stockholm-historical-monthly-temperature-3.
Acknowledgements
Our dear colleague and co-author Halina Kowalewska-Kalkowska, who passed away before the manuscript was submitted, will remain in our memory. We would like to express our gratitude for her important contribution to this article and our sadness that our collaboration ended so soon. The research work of PO and RP was supported by the grant ‘The occurrence of extreme weather, climate and water events in Poland from the 11th to 18th centuries in the light of multiproxy data,’ funded by the National Science Centre, Poland (Grant No 2020/37/B/ST10/00710). In addition, work of PO and RP was done within the Research University – Initiative of Excellence: the Emerging Field ‘Global Environmental Changes’ and the Climate Change Research Unit at Nicolaus Copernicus University in Toruń. The research work of HKK was supported through the funds of the Institute of Marine and Environmental Sciences, University of Szczecin, Poland.
Competing interests
The authors declare no competing interests.
Appendix: Parametrisation of heat fluxes at the water–atmosphere–sea-ice boundary in the PM3D model
The total heat flux across the water–atmosphere boundary 
${F_{tot,aw}}$ and the sea-ice–atmosphere boundary 
${F_{tot,ai}}$ is the sum of five components: (1) short-wave solar radiation passing through the water surface 
${F_S}$ or ice 
${F_{Si}}$; (2) long-wave radiation leaving the water 
${F_{Lw}}$ or ice 
${F_{Li}}$ and (3) incoming from the atmosphere 
${F_{La}}$; (4) the latent heat flux of the water 
${F_e}$ or the sublimation of ice 
${F_{e,i}}$; and (5) perceived heat (the sensible heat flux) on the surface of the water and ice: 
${F_{wt}}$, 
${F_{wt,i}}$:
\begin{equation}{F_{tot,aw}} = {F_{Sw}} + {F_{Lw}} + {F_{La}} + {F_e} + {F_{wt}}\end{equation}
\begin{equation}{F_{tot,ai}} = {F_{Si}} + {F_{Li}} + {F_{La}} + {F_{e,i}} + {F_{wt,i}}\end{equation}Short-wave radiation fluxes: 
${F_{Sw}}$, 
${F_{Si}}$ account for the flux of solar radiation reaching the sea surface 
${F_{S,in}}$ and the albedo of water 
${\alpha _w}$ or ice 
${\alpha _i}$, respectively, and:
\begin{equation}{F_{Sw}} = \left( {1 - {\alpha _w}} \right){F_{S,in}}\end{equation}
\begin{equation}{F_{Si}} = \left( {1 - {\alpha _i}} \right){F_{S,in}}\end{equation}
${F_{S,in}}$ was determined according to Kreżel (Reference Kreżel1997). The long-wave radiation of water and atmosphere is determined using the Stephan–Boltzmann law:
\begin{equation}{F_{Lw}} = - { \in _w}{C_{SB}}T_w^4\end{equation}
\begin{equation}{F_{Li}} = - { \in _i}{C_{SB}}T_i^4\end{equation}
\begin{equation}{F_{La}} = { \in _a}{C_{SB}}T_a^4\left( {1 + {d_a}e_a^{0.5}} \right)\left( {1 + {d_C}{N^2}} \right)\end{equation}where 
${C_{{\text{SB}}}} = 5.67{\text{ }} \times \;{10^{ - 8}}{\text{ W }}{{\text{m}}^{ - 2{\text{ }}}}{{\text{K}}^{{\text{ }} - 4}}$ (the Stephan–Boltzmann constant), 
${ \in _w}$ and 
${ \in _a}$ is the emissivity of the water and the atmosphere, 
${T_w}$ is water surface temperature, 
${T_i}$ is ice surface temperature, 
${e_a}$ is water vapour pressure at atmospheric temperature 
${T_a}$, N is degree of cloudiness expressed on a scale from 0 to 1 and 
${d_a}$ and 
${d_C}$ are empirical factors whose values were determined during model calibration. The latent heat flux of vaporisation:
\begin{equation}{F_e} = {c_e}{\rho _a}{L_e}{U_a}\left( {{q_a} - {q_w}} \right)\end{equation}
\begin{equation}{F_{e,i}} = {c_e}{\rho _a}{L_s}{U_a}\left( {{q_a} - {q_i}} \right)\end{equation}
\begin{equation}{q_a} = \frac{{0.622{\text{ }}{e_a}}}{{p - 0.378{e_a}}}{\text{ }}\end{equation}
\begin{equation}{q_w} = \frac{{0.622{\text{ }}E\left( {{T_w}} \right)}}{{p - 0.378{\text{ }}E\left( {{T_w}} \right)}}\end{equation}
\begin{equation}{q_i} = \frac{{0.622{\text{ }}E\left( {{T_i}} \right)}}{{p - 0.378{\text{ }}E\left( {{T_i}} \right)}}\end{equation}where 
${c_e}$ is an empirical coefficient, 
${\rho _a}$ is air density, 
${L_e}$ is latent heat of water vaporisation, 
${L_s}$ is latent heat of ice sublimation, 
${U_a}$ is wind speed at 10 m above sea surface (ms˗1), 
${q_a}$ is absolute humidity of the atmosphere and 
${q_w}$ is absolute humidity at water temperature. p is atmospheric pressure, and 
$E\left( {{T_x}} \right)$ expresses water vapour pressure at water or ice temperature calculated using Clapeyron relations.
The sensible heat flux depends on the temperature difference between the air and the water or ice:
\begin{equation}{F_{wt}} = \left\{ {\begin{array}{*{20}{c}}
{{c_{wts}}{\rho _a}{C_w}{U_a}\left( {{T_a} - {T_w}} \right){ }for{ }{T_a} \geqslant {T_w}} \\
{{c_{wtu}}{\rho _a}{C_w}{U_a}\left( {{T_a} - {T_w}} \right){ }for{ }{T_a} \lt {T_w}}
\end{array}} \right.\end{equation}
\begin{equation}{F_{wt,i}} = \left\{ {\begin{array}{*{20}{c}}
{{c_{wts}}{\rho _a}{C_i}{U_a}\left( {{T_a} - {T_i}} \right){ }for{ }{T_a} \geqslant {T_i}} \\
{{c_{wtu}}{\rho _a}{C_i}{U_a}\left( {{T_a} - {T_i}} \right){ }for{ }{T_a} \lt {T_i}}
\end{array}} \right.\end{equation}where 
${C_w}$, 
${C_i}$ are the specific heat of water and ice, 
${c_{wts}}$, 
${c_{wtu}}$ are empirical coefficients characterising heat transfer for a stable and unstable atmosphere, respectively.
The heat flux passing through sea ice is the sum of three components:
\begin{equation}{F_{tot,iw}} = {F_{h,i}} + {F_{s,b}} + {F_{t,b}}\end{equation} The heat flux through the ice of thickness 
${h_i}$ as a result of the thermal conductivity of the ice:
\begin{equation}{F_{h,i}} = {\kappa _T}\frac{{\partial T}}{{\partial z}} \approx \frac{{{\kappa _T}}}{{{h_i}}}\left( {{T_i} - {T_f}} \right)\end{equation} where 
${\kappa _T}$ is the thermal conductivity of the ice, 
${T_f}$ is the temperature of the lower surface of the ice, which is assumed to be equal to the freezing temperature of water of a given salinity. The flux of short-wave solar radiation passing through the ice was determined using the formula:
\begin{equation}{F_{s,b}} = {F_{Si}}{f_{sol}}{e^{ - \kappa {h_i}}}\end{equation} where 
${f_{sol}}$ and 
$\kappa $ are ice penetration and light attenuation, respectively. The turbulent heat flux under the ice depends on the temperature difference between and relative velocity of water and ice:
\begin{equation}{F_{t,b}} = {c_{is}}{\rho _i}{C_w}\left| {{{\mathbf{U}}_w} - {{\mathbf{U}}_i}} \right|\left( {{T_w} - {T_f}} \right)\end{equation} where 
${c_{is}}$ is the bulk heat transfer coefficient from the water into the ice, 
${\rho _i}$ is ice density, 
${{\mathbf{U}}_w}$ is the velocity of water under the ice, and 
${{\mathbf{U}}_i}$ is ice velocity.
Empirical coefficients used for calculations of heat fluxes at the boundary of water, ice and atmosphere.
* θ denotes the sun’s position above the horizon (in degrees).
For a complete description of sea-ice thermodynamics and dynamics, as well as how it is coupled to the hydrodynamic model, see Herman and others (Reference Herman, Jędrasik and Kowalewski2011). Here, only a description of the parameterisation of heat fluxes with variable values of some coefficients has been provided.
Open access
