Why can modeling the Baltic Sea ice in studying air temperatures at the end of the Little Ice Age help predict future climate change? Are there other ways to model climate change? What role do ice and snow play in the Earth's atmosphere? What are the overall cycles of climate change in atmospheric chemistry from the beginning to the present?

Baltic Sea ice coverage was modelled using a sea-ice thermodynamics and dynamics model coupled with a three-dimensional PM3D hydrodynamic model. The validation for 1958–2007 showed the modelled maximum ice extents agree well with observations (r=0.97) and the ice thickness less so, but satisfactory for most stations (r>0.8). This enabled the production of Cumulative Ice Thickness (CIT) maps and the determination of the spatial variation in sea ice extent in the Baltic over the analysed period for four air temperature scenarios with a constant value reduction. This showed the spatial sensitivity of ice cover dynamics to temperature changes and allowed to distinct regions with different impact of change in temperature on CIT. The simulation for temperature of 2 °C lower than 1958–2007 was consistent with the reconstruction of maximum ice extents in the entire Baltic Sea for the end of the Little Ice Age (LIA) (1721–1860). For the western Baltic the compliance was highest for temperature reduced by 3 °C and 4 °C. This indicates that climatic conditions may have differed between This is an Open Access article, distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike licence (http://creativecommons.org/licenses/by-ncsa/4.0), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the same Creative Commons licence is used to distribute the re-used or adapted article and the original article is properly cited. The written permission of Cambridge University Press must be obtained prior to any commercial use. individual regions of the Baltic during the LIA, and the air temperature anomaly in the western Baltic may have been greater than indicated by previous studies. 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 approximately 200 years (the Transitional Period), which initiated the Little Ice Age (LIA) that predominated until the mid-19th century (Brázdil and others, 2005; Niedźwiedź and others, 2015). The mean winter air temperature in the Baltic Sea oscillated around ˗4 °C, i.e. about 2 °C colder than present conditions (Eriksson and others, 2007). Similar differences in winter temperature between the LIA and the modern period have been obtained for Poland by Przybylak and others (2005) and Przybylak (2016), and for Central and Northern Europe (Stockholm) by Brázdil and others (2010). Mean thermal conditions in the LIA showed some variability, so Eriksson and others (2007) 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 about 1 to 2 °C cooler than now only for the first hundred years of the LIA and in the 19th century (Büntgen and others, 2011). In Europe, the average annual air temperatures were about 1.5 °C lower in the LIA than now, and even 2 °C lower in extremely cold years (PAGES 2k Consortium, Supplementary information, 2013, see Figs S1 and S2 there). 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, 2013). 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, 1915). 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 (1555); Chronologia vetusta ab anno 1298 ad annum 1475,1818; Korner (1895) ; Olaus, (1772). 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, 2011). Over longer periods, ice extent on the Baltic Sea is closely related to climate change. Changes in the maximum ice extent 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, 1996; Omstedt and others, 2004). As shown by reconstructions of maximum ice extent in the Baltic since 1720 (Seinä and Palosuo, 1996 and https://www.ilmatieteenlaitos.fi/jaatilastot, 2021), the maximum ice extent 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 1999). The North Atlantic Oscillation (NAO) reflects the climatic variability of winters in the North Atlantic and Europe (Walker and Bliss, 1932). 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, 1995; Cook and others, 2019).The macro-scale 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, 2010; Lehmann and others, 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, 1-D models (Lepparanta, 1993, Yu and others, 2015) 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 modelled in many numerical studies, starting in the 1980s (Leppäranta, 1981; Omstedt, 1990; Leppäranta and Zhang, 1992; Haapala and Leppäranta, 1996; Omstedt and Nyberg, 1996; Haapala, 2000; Cheng and others, 2006). In subsequent years, increasingly advanced models have enabled the simulation of Baltic Sea ice conditions, including CRO (Meier, 1999), RCAO (Dӧscher and others, 2002), HAMSOM (Schrum and others, 2003), RCA4-NEMO (Dieterich and others, 2013) and HIROMB (Lӧptien and Axell, 2014). 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 (2011) and Ryabchenko and others (2010). Recently, Pemberton and others (2017) 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, 2018) and CESM (Jakacki and Meler, 2019). 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, 2018), 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, 2001; Meier, 2002; Meier and others, 2004; Meier, 2006; Luomaranta and others, 2014; Höglund and others, 2017). However, to date, few have covered earlier historical periods. Meier and Krauker (2002) simulated maximum ice extent in the Baltic for the period 1902–98 using the Rossby Centre coupled ice–ocean model (RCO). Hansson and Omstedt (2008) 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 (2012) carried out model simulations of maximum ice extent in the Baltic region in the MWP and LIA based on the RCO (Rossby Centre Ocean model). 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, 2009). The large inflow of fresh water, which reaches 6,670 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, 2009) and averages 7.36 psu (Meier and Kauker, 2003), 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, 2009). 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, 2009).2.2 Maximum ice extent 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 maximum ice extents in the Baltic published by Seinä and Palosuo (1996) 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 (Dec–Mar) from Helsinki, Stockholm and Mariehamn were used in reconstructing the maximum ice extent. Coastal data were also used to estimate ice conditions in the open sea (Seinä and Palosuo, 1996). According to Vihma and Haapala (2009), the greatest inaccuracies in estimating maximum ice extent 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, maximum ice extent is largely estimated using synthetic aperture radar imagery (Mäkynen and others, 2020). The maximum ice extent 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.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). Figure 2. 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)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 1770, 1788). 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, 2019). This observation series was supplemented with data from other subdaily air temperature measurements in: (i) Gdańsk (1739–72: Hanov, 1772), (ii) Wrocław and Lower Silesia (1717–26: Kanold, 1725; and 1727– 30: Büchner, 1731, 1732, 1733, 1734), (iii) Toruń (1760–64: Pospieszyńska and Przybylak, 2010), and Warsaw (1760–62: Rojecki, 1965). In addition, for Toruń and surrounding area, weather descriptions by J. Richtsteig for 1704–30 (Richtsteig, 1704–30), D. Brauer for 1719– 50 (Brauer, 1719–50) and J.H. Zernecke for 1700–25 (Zernecke, 1727) 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ś (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 (1994). Winters of indices ˗3 and ˗2 were usually characterised by long freezes (about 2–3 months) resulting in the appearance of ice in the Baltic Sea. Next, the indexation of winters was also reduced to five-point and three-point scales to facilitate comparison against the classifications of winters and ice conditions proposed by Seinä and Palosuo (1996) and Koslovsky and Glaser (1995), 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, 2000). It was reported that, especially in winter, temperature values are highly correlated across the entirety of Poland (Kożuchowski and Żmudzka, 2003; Przybylak and others, 2014), 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, 2010). The severity of winters was assessed according to the ranges of average winter air temperature used in works by Przybylak and others (2004, 2005) 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 (1996). The modified ranges of average winter air temperature in Warsaw corresponding to the five categories of winter severity are given below: Extremely mild: > 0.0 °C Mild: 0.0 °C – ˗1.9 °C Average: ˗2.0 °C – ˗4.6 °C Severe: ˗4.7 °C – ˗6.6 °C Extremely severe: < ˗6.6 °C To assess the nature of winters according to the three-grade classification applied by Koslovsky and Glaser (1995), 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, 2011) coupled with a three-dimensional PM3D hydrodynamic model (Kowalewski and Kowalewska-Kalkowska, 2017). The PM3D model is an enhanced version of the M3D model (Kowalewski, 1997) and is based on the POM model (Princeton Ocean Model) developed by Blumberg and Mellor (1987). 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 three nautical miles (approx. 5.5 km). The vertical division used the sigmatransformation 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, 1997; Feser and others, 2001) with a one-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 (2019). 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 (2011). 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 maximum ice extent (MIE) can be calculated in different ways. For example, Pemberton and others (2017) 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: ∑ (1) where is the ice concentration in the i th computational grid node, is the area of the i th 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 Loewe, 1994). 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: ∑ (2) where A is the ice concentration, h is the ice thickness [m], j is the index of winter season days, 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: ∑ ∑ (3) 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: ∑ (4) 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 Fig. 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 of greater than 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 approx. 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.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 km2 and 60×103 km2 (almost seven times less). This sequence of extremely severe and mild winters within a two-year period indicates the significant dynamics of ice phenomena in the Baltic Sea.In analysing ice in the Western Baltic, Koslowski and Glaser (1995) 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 (2019) 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 maximum ice extents 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, 1996) 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 maximum ice extent. 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.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 maximum ice extents in the Baltic Sea (BIAS = -13×103 km2 ), which is associated with the aforementioned tendency for the first decades of the simulation.The comparison of the maximum ice extent 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 Krauker, 2002) 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, 2008), 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 (2012) showed that the RCO model underestimated the value of the annual maximum ice extent in the years 1961–2000 by over 20×103 km2 compared to data published by Seinä and Palosuo (1993). Simulations of ice conditions with the NEMO-LIM3.6 model (Pemberton and others, 2017) in the similar period 1961–2006 were shown to be good reflections of sea-ice extent (understood as the area of seaice 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 Paluso, 1996) 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 (Figs. 7a and 7b).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),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 (approx. 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, 2000; Meier, 2000; Löptien and others, 2013; Pemberton and others, 2017). Simulations of ice thickness at Swedish stations using the RCO-HELMI (Löptien and others, 2013) and NEMO-LIM3.6 (Pemberton and others, 2017) 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, 2017) 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, 2002, Moberg, 2021) 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. Figure 8. Time series of decadal average winter air temperature (December-February) in Stockholm in 1760–2020 (source: Moberg, 2021).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 maximum ice extent (Seinä and Palosuo, 1996) that divided winters into: extremely mild (383×103 km2 ). In addition, means, standard deviations and medians were determined for the maximum ice extents presented in the work of Seinä and Palosuo (1996) (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, 1996, supplemented by FMI observations) and the model. The standard deviations were very similar, but the means and medians of the maximum ice extents produced by the model were lower than those observed. Of the simulation variants, the statistical characteristics of the maximum ice extents 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 maximum ice extents 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. Table 2. Percentage share of winters with various degrees of severity and other statistical characteristics calculated based on the reconstruction (Seinä and Palosuo, 1996, supplemented by FMI observations for 1996–2007) and for model simulations: reference (v0) and simulations with air temperature lowered by 1, 2, 3 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, 2, 3 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 and 4 °C, the increase in CIT is more visible in the middle of the Baltic Sea than in Kattegat.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 approximately 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 approximately 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.Koslowski and Loewe (1994) 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 cumulative ice thicknesses presented in that work, Koslowski-Glaser (1995) 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 cumulative ice thickness, 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 cumulative ice thickness to be weak (less than 2 m·d), moderate (2 to 6 m·d) and severe (more than 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). For the winter severity classes they distinguished, Koslowski and Glaser (1995) 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 (Table 4) 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 the Table 4 were worked out. Both the division of winters into five categories proposed by Seinä and Palosuo (1996) and into three categories according to Koslovsky and Glaser (1995) were used. 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 (1996) in terms of frequency of extremely mild (8% new classification, 7% Seinä and Palosuo (1996)), 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 (1996) 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 4a). However, greater discrepancies were found with the classification of winters proposed by Koslovsky and Glaser (1995), regardless of research period (Table 4b). High compliances were only obtained for the incidence of severe winters in 1761–1860 and moderate winters in 1721–60 (see Tables 3 and 4b). 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 Koslovsky and Glaser (1995) 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 (2012) conducted model simulations of the maximum ice extent 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 three-dimensional PM3 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 cumulative ice thickness 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 about 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) guite 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 (1996) in terms of mean and median values of maximum ice extent. 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 approximately 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, 2002). The model simulation of cumulative ice thicknesses in the western Baltic Sea for reduced air temperature values makes possible to use a winter severity classification method used by Koslowski and Glaser (1995) 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 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 (1996) and much worse compliance with the Koslowski and Glaser classification (1995) 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.

References :

Abel W (2013) Agricultural Fluctuations in Europe: From the Thirteenth to Twentieth Centuries. Taylor and Francis (doi: 10.4324/9780203816933) Blumberg AF and Mellor GLA (1987) A description of a three-dimensional coastal ocean circulation model. In: Three-dimensional coastal ocean models; Heaps NS (Ed.); American Geophysical Union, 1–16 (doi: 10.1029/CO004p0001) Brázdil R and 6 others (2010) European climate of the past 500 years: new challenges for historical climatology. Climatic Change, 101, 7 – 40 (doi: 10.1007/s10584-009-9793-y)Brázdil R, Pfister CH, Wanner H, Von Storch H, Luterbacher J (2005) Historical climatology in Europe—the state of the art. Climatic Change, 70, 363–430 (doi:10.1007/s10584-005- 5924-1) Büchner AE (1731, 1732, 1733, 1734) Miscellanea Physico-Medico-Mathematica Brauer D (1719–1750) Brauer, Geschichte der Stadt Thorn, vol. 1–3, Archiwum Państwowe w Toruniu, Kat. II, nr 54, 54a, 54b Büntgen U, Raible C, Frank D, Helama S, Cunningham L, Hofer D, Nievergelt D, Verstege A, Stenseth N, and Esper J (2011) Causes and consequences of past andprojected Scandinavian summer temperatures, 500–2100 AD. PLoS ONE 6:e25133 (doi:10.1371/journal.pone.0025133) Chen D, Hellström C (1999) The North Atlantic Oscillation and Swedish temperature variability from 1865–1994. Tellus A: Dynamic Meteorology and Oceanography, 51, 505–516 (doi: 10.3402/tellusa.v51i4.14086) Cheng B., Vihma T., Pirazzini R., & Granskog M.A. (2006) Modelling of superimposed ice formation during the spring snowmelt period in the Baltic Sea. Annals of Glaciology 44, 139–146 (doi:10.3189/172756406781811277 )Chronologia vetusta ab anno 1298 ad annum 1475 (1818) In: Scriptores rerum Svecicarum medii aevi, ed. Fant EM, T. 1, Upsaliae Cook ER, Kushnir Y, Smerdon JE, Williams AP, Anchukaitis KJ, and Wahl ER (2019) A Euro-Mediterranean tree-ring reconstruction of the winter NAO index since 910 C.E. Climate Dynamics 53, 1567–1580 (doi:10.1007/s00382-019-04696-2) Korner H (1895) Cronica novella de secundo opere bis 1423. In: Die Chronica Novella des Hermann Korner, ed. Schwalm J, Göttingen: Vandenhoeck & Ruprecht. doi: 10.17192/eb2010.0073 Dieterich C, Schimanke S, Wang S, Väli G, Liu Y, Hordoir R, Axell L, Höglund A, and Meier HEM (2013) Evaluation of the SMHI coupled atmosphere-ice-ocean model RCA4- NEMO. SMHI Report Oceanography, 47 (https://www.divaportal.org/smash/get/diva2:947918/FULLTEXT01.pdf) Döscher R, Willén U, Jones C, Rutgersson A, Meier HEM, Hansson U, and Graham LP (2002) The development of the regional coupled ocean-atmosphere model RCAO. Boreal

APPENDIX A: Parametrization of heat fluxes at the water–atmosphere– sea-ice boundary in the PM3D model The total heat flux across the water–atmosphere boundary and the sea-ice– atmosphere boundary is the sum of five components: 1. short-wave solar radiation passing through the water surface or ice ; 2. long-wave radiation leaving the water or ice and 3. incoming from the atmosphere ; 4. the latent heat flux of the water or the sublimation of ice ; and 5. perceived heat (the sensible heat flux) on the surface of the water and ice: , : (A1) (A2) Short-wave radiation fluxes: , account for the flux of solar radiation reaching the sea surface and the albedo of water or ice , respectively, and: (A3) (A4) was determined according to Kreżel (1997). The long-wave radiation of water and atmosphere is determined using the Stephan–Boltzmann law: (A5) ((A7) where (the Stephan–Boltzmann constant), and is the emissivity of the water and the atmosphere, is water surface temperature, is ice surface temperature, is water vapour pressure at atmospheric temperature , N is degree of cloudiness expressed on a scale from 0 to 1, and and are empirical factors whose values were determined during model calibration. The latent heat flux of vaporisation: (A8) (A9) (A10)(A7) where (the Stephan–Boltzmann constant), and is the emissivity of the water and the atmosphere, is water surface temperature, is ice surface temperature, is water vapour pressure at atmospheric temperature , N is degree of cloudiness expressed on a scale from 0 to 1, and and are empirical factors whose values were determined during model calibration. The latent heat flux of vaporisation: (A8) (A9) (A10)formula: (A17) where and 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: | |( ) (A18) where is the bulk heat transfer coefficient from the water into the ice, is ice density, is the velocity of water under the ice, and is ice velocity. Table A1. Empirical coefficients used for calculations of heat fluxes at the boundary of water, ice and atmosphere Symbol Original value or formula (Herman and others (2011) Values/formulas used in this work Units 2.4/(90−θ) ∗ 2.4/(90−θ) ∗αi 0.7 0.7 – 0.985 0.97 – 0.685 0.735 – 0.03 0.01 hPa−1/2 0.36 0.17 – – – – 2.0 3.0 Wm−1 K −1 0.1 0.1 – 1.5 1.5 m −1 – ∗θ 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 (2011). Here, only a description of the parameterisation of heat fluxes with variable values of some coefficients has been provided.