An Evaluation and Implementation of the Regional Coupled Ice-Ocean Model of the Baltic Sea

A three dimensional, regional coupled ice-ocean model based on the open-source Community Earth System Model has been developed and implemented for the Baltic Sea. The model consists of 66 vertical levels and has a horizontal resolution of approx. 2.3 km. The paper focuses on sea ice component results but the main changes have been introduced in the ocean part of the coupled model. The hydrodynamic part, being one of the most important components, has been also presented and validated. The ice model results were validated against the radar and satellite data, and the method of validation based on 5 probability was introduced. In the last two decades satellite and model results show an increase in the ice extent over the whole Baltic Sea, which is an evidence of a negative trend in air temperature in recent decades and increasing of winter discharge from the catchment area.

In this work, the recently developed Community Earth System Model (CESM, Hurrell et al. (2013) has been used and implemented for the whole Baltic Sea.CESM is a descendant of the Community Climate Model (CCM, Williamson et al. (1987); Hack et al. (1993)), the Climate System Model (CSM) and the Community Climate System Model (CCSM version 3, Blackmon et al. (2001) developed by the National Centre for Atmospheric Research (NCAR).Up to the mid-1990s, the CCM was of limited use as a climate model because it did not include submodels of the global ocean and sea ice.
The first plan to develop and use a climate system model (CSM) including the atmosphere, land surface, ocean, and sea ice was proposed to the National Science Foundation (NSF) in 1994 (Blackmon et al. , 2001) The recent releases were also improved by the University Centre for Atmospheric Research (UCAR).The latest CCSM version 4 consists of four separate components, atmosphere (atm), ocean (ocn), sea-ice (ice) and land (lnd), coupled by a central coupler component.The main difference between current version of CCSM and CESM is that new biogeochemistry capabilities were added in the CESM model.It means that CESM is only a newer version of CCSM.It is also important, that CESM includes CCSM -it is possible to run simulations implemented in CCSM using CESM model.CESM has a set of configurations and selected ones have been already scientifically validated.All configurations are based on global grids, except the RACM/RASM project (Regional Arctic Climate/ System Model) that develops all active components in a coupled, regional Pan-Arctic climate model ( (Maslowski et al., 2009), (Cassano et. al., 2010), (Maslowski et al., 2012)).Since resolutions of the provided configurations are not suitable for the Baltic Sea, one of challenges in the presented implementation was to adequately increase both horizontal and vertical resolution.
For the purpose of this study the CESM was adopted for the Baltic Sea.The implemented version includes two active components, ocean (Parallel Ocean Program, POP) and ice (Community Ice CodE, CICE), and two data models representing atmosphere (datm) and land (dlnd).The active (dynamical) components are fully prognostic models, whereas the basic function of data models is reading external data (interpolated outside of the model), modifying that data and sending them to the central coupler.The central coupler (cpl7) is responsible for the coordination and the flow of information between all components.The coupler and other models have no fundamental knowledge of whether another component is fully active or just a data model.The implemented setup of the Baltic Sea model is shown in Fig. 1.We call this model B-CESM.
In the above configuration, different data sets have been used as external forcings.The river runoff data from the Baltic HYdrological Predictions for the Environment model (Balt-Hype, (Lindström et al., 2010)), covering the years  were provided by the Swedish Meteorological and Hydrological Institute (SMHI).The river runoff climatology was calculated from available period and employed in the model for missing periods.In total, data for 75 rivers have been implemented in applied approach is to treat a river runoff as precipitation, which changes salinity of the uppermost cells.Another way is to add the freshwater input as precipitation but instead salinity change, the freshwater volume is added to the model upper cells.This requires removal of the additional water volume at the lateral boundary of a regional model as was done in the presented Baltic Sea model.Our approach is non-standard and has been implemented only for the purpose of this model for the Vistula River.A special channel has been added to the model domain and precipitation has been released only within this channel.The sea level difference between the channel and the Gulf of Gdansk assured velocity of the freshwater flow, which could not be obtained, if the river runoff was implemented as precipitation over the sea.The implemented channel approach gives much better results for large rivers, such as Vistula or Lena, while for smaller rivers including the river runoff as precipitation is commonly used.
The presented results were obtained with the Baltic Sea model using atmospheric forcing from ERA Interim reanalysis (ERAi, ECMWF, 1979-present, 0.75 degree horizontal resolution, (Berrisford et al. , 2011)) and data from the local forecast model mean sea level pressure, short and long wave radiation downward, total snow and rain precipitation.
The availability of global atmospheric reanalyses permits to use them from different data sources.For example in the case of NCEP (National Centers for Environmental Prediction) (Kalnay, 1996;Saha et al., 2010), JMA (Japan Meteorological Agency) (Onogi et al., 2007), NASA (National Aeronautics and Space Administration) (Rienecker et al., 2011;Schubert et al., 1993), and ECMWF ( (Dee et al. , 2011;Berrisford et al. , 2011), it is possible to use data that cover the Baltic Sea area and all of them have resolution lower then one degree.Also regional analysis made by High Resolution Limited-Area Model (HIRLAM) are provided with 22 km horizontal resolution.This reanalysis has been done using ERAi data as boundary conditions.In our case we use ERAi because it covers time of integration (which is limited by observational sea level data for Goteborg station) and runoff data are available for the same period.Also data from ECMWF could be used in the Baltic Sea modeling with reasonable accuracy (Omstedt &Wesslander , 2005).
All models (components) are governed by the central coupler (cpl7), based on the Model Coupling Toolkit (MCT).It controls the exchange of data between the individual models and all time steps.The coupling time step for all models has been set to 600 seconds and the same value was also used as the ice model time step.The ocean model has its own system that controls its time steps, which are not uniform.The requested time step of the ocean model has been equal to 150 seconds and we can assume that it has been the longest step of the POP.

Ocean model description and adaptation to the Baltic Sea
The ocean component of the Baltic Sea model in the presented configuration is the Parallel Ocean Program, a well-known, 3-dimensional, z-coordinate general circulation model ( (Smith & Gent, 2004)).The model uses hydrostatic and Boussinesq approximations for the primitive equations of the fluid motion on the orthogonal spherical grid.A free surface is implemented as pressure at the top boundary layer.A detailed description of the POP model and improvements can be found in numerous papers, e.g.(Bryan, 1969;Semtner, 1974;Cox , 1984;Killworth et al., 1991;Mintz & Semtner, 1977).

Initial state of the model
Both active models (ice and ocean) adopt a rotated spherical grid with horizontal resolution of 1/48 degree (approx.2.3 km).To obtain the same size of all cells in a horizontal plane, the equator of the grid is located in the center of the domain.The model bathymetry is based on the ETOPO1 (Earth Topography Digital Dataset), a global 1 arc-minute relief model of the Earth's surface (Amante& Eakins , 2009).The B-CESM has 66 vertical levels, of which the first 50 levels are 5 m thick, whereas thickness increases in deeper layers.The number of levels has been chosen to optimize the computational cost in the model, at the same time covering most of the Baltic Sea with identical cells.The compromise resulted in 66 vertical levels shown in  Initial temperature and salinity fields for January 1990 were interpolated over the model domain.The model started from nomotion state, which is also called a 'cold start'.The model was spun up for three years (1990)(1991)(1992) and then went back to 1 January 1990.In the Baltic Sea strong autumn storms drive very deep turbulent mixing, therefore three years seems enough for the ocean model spin-up.In case of ice cover or heat balance, memory is closed to one year.The salinity memory of the Baltic Sea is over 30 years (Omstedt & Hansson, 2006).But in case of given observed (or reanalyzed) salinity, could be realistically modeled for the whole Baltic Sea.Located in the narrow and shallow connection between Skagerrak and Kattegat, Göteborg provides a convenient point for applying boundary conditions.Moreover, sea level data for this location are openly available for over twenty years (1990today).To apply sea level assimilation for the Goteborg area, it has been necessary to introduce a modification in the barotropic equation, originally formulated as (Dukowicz & Smith , 1994): where H is the total depth, f is the Coriolis parameter, (u, v) are the barotropic (vertically averaged) velocities, g is the gravity acceleration, G x and G y are external forcings and η is the surface pressure.Symbols ∂ t ,∂ x ,∂ y represent time and space derivatives.The RHS of Eq.( 1) represents x and y gradient components of the surface pressure.This part of the barotropic equation has been modified to incorporate sea surface height (SSH) at the boundary area by adding the gradient of difference between the model sea level and the SSH measured in Göteborg.The Cressman analysis scheme (Cressman, 1959), which defines the weighting function in terms of a radius of influence, is widely used for simple assimilation systems.In our case the Hamming window (w h (x, y)), which is very often used in signal processing, has been used instead of the Cressman weighting function.An advantage of this approach is that the Hamming window is represented by simple harmonic function that is the basic solution of the wave equation.The modified barotropic Eq.( 1) has now the form: where η , is the difference between the modeled and measured sea level in Göteborg.
The Sommerfeld radiation conditions (Sommerfeld, 1949) are most commonly used at the boundary.In our case the condition at the northern boundary is formulated as: Where ψ is velocity, barotropic component of velocity or sea surface height at the lateral boundary in the model, and c is the wave speed.Two approaches are often applied to find the c parameter.One takes into account the barotropic wave speed, which is constant and depends on the depth and gravitational acceleration.The second approach, introduced by Orlanski ( 1976 (3) provides a prognostic variable at the boundary (Kantha & Clayson, 2000): where In Eq.( 4) and ( 5) n represents time step and jM is the northern boundary index.In the B-CESM model Eq.( 4) and ( 5) have been applied in Kattegat.
Comparison of the sea level measured in Göteborg and obtained from the model after implementing the modified boundary condition is shown in Fig. 3.Even if the modeled sea level is lower than the measured one, the applied modification of the barotropic equation produces realistic sea level variability.The main reason of discrepancy in the sea level amplitude is pressure averaging, turned on in the POP model to avoid the model becoming unstable.Pressure averaging allows for longer time step in the ocean model and reasonable integration time.1990 1991 1992 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006  Because the North Sea is the source of salt in the Baltic Sea, the flow through the Danish Straits is the most important driver of the salt exchange between Kattegat and the South-West Baltic.The two main straits, Sund and Langeland Belt, are the main passages for the salt transport into the Baltic (in example Lehmann et al. (2004); Mohrholz et al. (2015)).The flow through Sund is mostly driven by the sea level pressure difference between the southern and northern entrances, and can be calculated from the semi-empirical expression: where ∆η is the sea level difference between the ends of Sund (Skanor and Hornbaek), K f is specific resistance, and Q is the flow through the strait.
The long-term sea level measurements are available for two stations shown in Fig. 4. Based on measurements, the specific resistance in Sund was estimated on approximately 2.03 × 10 −10 .Figure 4 shows that flow rate depends on the sea level difference for the Sund.The black line represents measurements and the red points are from the model described herein.The nonlinear regression provides the modeled specific resistance Kf = 4.05e − 10.For comparison GETM (Grode, 2004)  The insert shows position of two stations at the northern and southern end of Sund and the section.
The dynamical state of the Baltic Sea is the result of balance between dense and salty waters incoming through the Danish Straits and the sea level driven by winds, salinity difference, precipitation and rivers.Because of it, the annually averaged sea level is expected to increase from the Danish Straits towards the Gulf of Bothnia .Figure 5 shows adjusted annual mean sea level at seven stations, located along the Baltic coast from west to east, obtained from three different numerical models: BS01, NB03 and GETM.BS01 and NB03 are the implementations of the HIROMB model with horizontal resolution of 1 and 3 nm respectively.GETM is a model set up by the Royal Danish Administration of Navigation and Hydrography and in this case it has a resolution of 1 nm (Grode, 2004).
The HIROMB models use meteorological forcing from HIRLAM (High Resolution Local Area Modelling, SMHI), while the GETM model is also driven by data obtained from HIRLAM, but from the Danish Meteorological Institute (DMI).The monitoring data were obtained from the height reference systems: Baltic HRS and Nordic Height System (NH60), which are 15 years long time series transformed from the Swedish Height System RH70 ( (Carlsson, 1997)).with observations.The averaged sea level varies between the individual models because of the different atmospheric forcing (in this case we use data from the ERA40 interim, ERA40 and Interdisciplinary Modeling Centre (ICM) at the University of Gdansk) and different resolutions.Additionally, pressure averaging in time has been implemented in the other models, resulting in smaller amplitude of the sea level variations.

Ocean model upper boundary layer
The conservation of heat at the top layer of the open ocean model requires (Yu at al, 2007): where subscripts net, sw, lw, lat, sen indicate net heat, solar short waves radiation, long wave radiation, latent and sensible heat flux, respectively.It has been assumed that fluxes are directed to the ocean model (the ice model will be described separately).In our model, short and long wave radiations are provided by datm model as external forcing (see Fig. 1).Latent and sensible heat fluxes are sent to the active components and are calculated by the central coupler based on the bulk formula: where E and H are latent and sensible heat fluxes, ρ A is the air density, ∆v is the velocity difference between the ocean surface and atmosphere, ∆q and ∆θ are specific humidity and temperature differences (between the atmosphere and the ocean surface) and Cp is the specific heat capacity.C E and C H are transfer coefficients between the ocean surface and atmosphere (sensible heat transfer coefficient is the Stanton number, the latent heat coefficient is the Dalton number), and they depend on stability: where ψ m and ψ h are the integration forms of the Monin-Obukhov similarity functions,Z A /Z 0m represents stability parameter, where C d is drag coefficient, S is the wind speed relative to the ocean surface and u A and u S are wind and ocean surface velocities, respectively (vector quantities indicated by bold characters).
In consequence, the vertical shear is maximal at the surface of the ocean and decreases with depth.The interaction between atmosphere and ocean produces a mixed layer, occupying the upper part of the ocean, where temperature and salinity can be approximated as constant.The mixed layer depth (MLD) depends mostly on wind speed and direction, as well as its temporal and spatial scales and in shallow waters also on bottom topography.An important part of the energy transfer between atmosphere and ocean upper layers is through turbulences.In the presented B-CESM model k-profile parameterization (kpp) was used (Large et al., 1994).
Due to the fact that these schemes are not available in the POP model (MY is for sigma coordinate models only), kpp has been modified to be more similar to MY based on Durski et al. (2004).The comparison of the two different schemes is not  the upper 5-m thick layer, whereas the satellite data represent sea temperature at the surface (usually the skin temperature).
Therefore, the modeled sea surface temperature is lower than the measured one in summer and higher in winter.Validation of the modeled SST against satellite data is complemented by a comparison of one-year long time series of temperature and salinity profiles from the HuvudskarOst station (location shown on Fig. 1).Despite numerous gaps in in-situ measurements, they can be still used for validation of the model results in terms of the vertical distribution of temperature and salinity and their The second active component in the presented Baltic Sea model is the Los Alamos sea ice model -Community Ice CodE (CICE).CICE is a descendant of the basin-scale dynamic-thermodynamic and thickness distribution ice model (Hibler, 1979(Hibler, , 1980)).The physics in the uncoupled ice model is the same as for the ice model used in the fully coupled system.The CICE model is based on the energy conserving thermodynamics of (Bitz & Lipscomb , 1999).It contains multiple layers in each thickness category and takes into account the effects of brine pockets within the sea ice.The ice dynamics utilizes the elasticviscous-plastic (EVP) rheology of (Hunke & Dukowicz, 1997, 2002).Sea ice ridging follows Rothrock (1975)

Fluxes at the ice-ocean and atmosphere-ice interfaces
Generally, heat and momentum fluxes, which influence the sea ice are similar to the ocean-atmosphere fluxes at the open ocean surface.The main difference is due to the ice albedo and ice area.The net energy flux from atmosphere to ice is calculated based on the heat budget: where i o is ice fraction and α is ice albedo.
Momentum flux required by ice model is based on the nonlinear integral boundary-layer theories (Hibler & Bryan, 1987).
Ocean and air momentum fluxes are calculated with the following equations: where u represents the sea ice velocity, U g is the geostrophic wind velocity (assumed to be higher than the ice motion velocity), U w is the ocean surface velocity, c a and c w are air and water drag coefficients, ρ a and ρ w are air and ocean density and phi and θ are air and water turning angles.

Ice model validation and results
Sea ice is the key element of the whole Baltic Sea system.Numerous satellite products exist that can be used   Since daily ice charts (IC), produced by the Finnish Ice Service during each ice season, are available only as images, it is difficult to compare them with model results.To obtain more relaible validation, the sea ice thickness reproduced by the B-CESM model was compared to a combined product of SAR data and thermodynamic ice model HIGHTSI (HIGHresolution Thermodynamic Snow/Ice model).HIGHTSI is a coupled snow and sea ice one-dimensional physical model targeted to determine the snow/ice surface temperature, in-snow/ice temperature and snow/ice thickness at a selected site (Karvonen et al., 2003(Karvonen et al., , 2007(Karvonen et al., , 2008)).To obtain an optimum balance between validation of sea ice parameters in the single location (precise but not very reliable) and validation for the sub-basins of the Baltic Sea (reliable but not sufficiently precise), in the next step we compare sea ice 18  The best agreement between the modeled and measured sea ice variables was obtained in the box located in the Bothnian Bay (Fig. 14, left panels).Ice frequency and averaged ice thickness from the simulation A almost ideally agrees with climatological data, while in the second and third simulations (B and C) reproduced ice thickness slightly higher than observed.Furthermore, there is a very good agreement between the measured and the modeled probabilities of sea ice occurrence in different thickness classes.The shapes of probabilities reflect the natural variability of sea ice thickness.In the beginning of winter thin ice appears.
Then, ice thickness increases, which is visible in higher probabilities of ice thickness in the range of 20-50 centimeters.After some time lag ice thickness reaches 50 centimeters and more.These time shifts are visible in the probability distributions.
Towards the end of winter melting processes result in decreasing thickness of sea ice cover.It is reflected in a visible symmetry 10 (quasi-symmetry) of probability that ice thickness is in the range of 0-20 cm.
In The average ice thickness in the Bothnian Sea simulated by the model is lower than obtained from climatology.We think, the differences in the ice probability shown for Bothnia Sea and Gulf of Finland suggest there could be problem with ice dynamics.
Ice cover in Bothnia Bay is much stable and influence of sea currents is smallest comparing to Bothnia Sea and Gulf of Finland.Long term variability of ice cover in the Baltic Sea is shown on Fig. 15A as maximum annual sea ice extent since the 18th century (Seinä & Palosuo (1996); Seinä et al. (2001Seinä et al. ( , 2006))).The black line depicts yearly values, the red line shows 15-year moving average and the modeled maximum sea ice extent has been overlaid for the last two decades (green line).
The maximum Baltic Sea annual ice extent has been slowly decreasing for the last two centuries.After a strong drop observed in the beginning of 1990s due to unknown reasons, the sea ice maximum extent has been increasing again.Figure 15 B shows the detailed comparison between modeled and climatological sea ice extent for this period.The linear regression between modeled and measured data provides the correlation coefficient higher than 0.9, as shown by the insert on Fig. 15B.
Furthermore, an increase in the ice cover is also visible in the last three decades.The main factor that could cause a visible positive trend of ice extent is temperature.However, global warming has been observed in recent decades and it could suggest that the Baltic Sea ice extent should have a negative tendency.Cohen et. al. (2014) shown that for the years 1990-2013 DJF surface temperature trends (per decade) are negative or neutral over the Baltic Sea area (Fig. 16) which could explain the trend.
Figure 16.Winter temperature trends for the most recent period from 1990 (Cohen et. al., 2014).
Another, much less important factor that could have influence on the ice extent is fresh water from the catchment area.The Baltic Sea runoff is an important part of the fresh water budget.In general, freshwater flux minus evaporation is positive and has a strong influence on the Baltic Sea circulation.Meier (2003) showed that the mean winter runoff to the Baltic Sea (January-  1990 1991 1992 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006

Figure 1 .
Figure 1.The Baltic Sea coupled model diagram UM run by the Interdisciplinary Modeling Centre at the University of Warsaw (ICM, approx.4 km horizontal resolution, 2010-2014).Atmospheric forcing variables required by the model includes: temperature at 2 meters, specific humidity at 2 meters, wind speed and direction at 10 meters, 3 Earth Syst.Dynam.Discuss., doi:10.5194/esd-2017-35,2017 Manuscript under review for journal Earth Syst.Dynam.Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.

Fig. 2
Fig. 2 together with the model domain and bathymetry (colors represent vertical levels).

Figure 2 .
Figure 2. Model domain, bathymetry (color scale in model levels) and vertical resolution (insert).
boundary conditions and their verificationAlthough the Baltic Sea is a semi-enclosed sea, a strong salinity difference between the North Sea and the Baltic Proper results in a difference in sea level of the order of 10 cm between the Kattegat and the South-West Baltic.Being a transit area, the Danish Straits play the key role in water exchange between the North Sea and the Baltic Sea.Therefore a proper implementation of Earth Syst.Dynam.Discuss., doi:10.5194/esd-2017-35,2017 Manuscript under review for journal Earth Syst.Dynam.Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.boundary conditions has strong influence upon the barotropic balance between the North Sea and the Baltic Sea and also upon the salinity distribution in the Baltic Sea.Two different approaches were used to implement lateral boundary conditions: a) sea level measured in Göteborg was adapted as a boundary condition in Kattegat, b) the Orlanski open boundary conditions (OBC) were enforced at the boundary area.
), was implemented also in our model.It uses data from the previous time step to calculate phase speed, which is applied in the current 6 Earth Syst.Dynam.Discuss., doi:10.5194/esd-2017-35,2017 Manuscript under review for journal Earth Syst.Dynam.Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.time step.Taking into account the varying c, Eq.

Figure 4 .
Figure 4. Flow rate vs. sea level difference in the Sund (red points -modeled, black line -calculated from the measured sea level difference).

Figure 5 .
Figure 5.Comparison of the adjusted annual mean sea level from three different models, observations and current Baltic Sea model.
0m is the roughness length.Wind -natural motion of the air -always exists over the ocean.Wind blowing over the ocean surface produces tangential Earth Syst.Dynam.Discuss., doi:10.5194/esd-2017-35,2017 Manuscript under review for journal Earth Syst.Dynam.Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.stress at the ocean-atmosphere interface, represented by a vertical flux of the horizontal momentum.Similar to heat fluxes, the horizontal momentum flux (wind stress τ o ) can be also described by a bulk formula:

Figure 8 Figure 7 .
Figure8shows a good agreement between observational data and model results.The model provides average temperature in

5Figure 8 .Figure 9 .
Figure 8.A comparison of the modeled and observed sea surface temperature at four different locations (based on daily data, results are for year 2012).
and Thorndike et al. (1975).CICE has several interacting components: a thermodynamic model that computes local growth rates of snow and ice due to vertical conductive, radiative and turbulent fluxes, along with snowfall; a model of ice dynamics, which predicts the velocity field of the ice pack based on a model of the material strength of the ice; a transport model that describes advection of the a real concentration, ice volumes and other state variables; and a ridging parameterization that transfers ice among thickness categories based on energetic balances and rates of strain.The CICE also has multiple thickness categories and ice thickness distribution evolves in time (CICE Scientific Reference).CICE has been successfully used in many regional and global adaptations (for exampleRinke et al. (2003); McLaren et al. (2006); Dzierzbicka & Jakacki (2013); Herman et al. (2011)).The model does not require a spin up because sea ice disappears in the Baltic Sea every summer.Horizontal resolution of 1/48 degree (about 2.3 km) is the same as in the ocean model.It has been configured for 5 vertical categories.The time step in the ice model of 600 seconds is the same as the coupler time step.
to validate the Baltic Sea model but only few are suitable for sea ice validation.Under the MyOcean framework (currently known as the Copernicus Marine Environment Monitoring Service) two sea ice datasets are available -the first from the Danish Meteorological Institute (DM) and the second provided by Finnish Meteorological Institute (FMI).The FMI data include sea ice concentration and thickness over the entire Baltic Sea and are available as daily data with the horizontal resolution of 1 km.DMI provides only sea ice concentration with 2 km resolution.Both observational data sets are based on SAR images.The operational sea ice service at FMI provides ice parameters, produced on a daily basis during the Baltic Sea ice season.The ice thickness chart (ITC) is a product based on the most recent available ice chart (IC) and a SAR image.The DMI data are also based on the operational product.

Figure 10
Figure10shows a comparison between the modeled and the measured ice thickness and concentration for two different points, mostly covered by ice.The first one (Fig.10, left panel) is located in the Bothnian Bay ( 24.2944 o E, 65.2023 o N), while the second one lays in the Gulf of Finland (28.6988 o E, 60.1510 o N).It is clearly visible that variability of ice concentration and thickness provided as operational products is much lower than modeled values.The cause of it is unknown, yet it is suspected that SAR fails to sense minor changes in ice thickness and concentration.Even if the SAR data do not fully resolve the shortterm variability, the ranges of measured and modeled sea ice concentration and thickness are similar as well as their changes on longer times scales.

Figure 10 .
Figure 10.Comparison of the measured (FMI, DMI) and modeled ice concentration and thickness for the Bothnia Bay (left) and the Gulf of Finland (right).
Figures 11 and 12 present a comparison of sea ice thickness simulated with the presented B-CESM model with available ice charts based on SAR data combined with the thermodynamical ice model.The modeled and derived from SAR/HIGHTSI mean level ice thickness do not match ideally but spatial patterns are similar in both cases.In particular a distribution and shapes of the open water areas are well reproduced by the Baltic Sea model.Ice thickness in the Gulf ofFinland is in the range between 10 and 30 cm in both cases.In the Gulf of Bothnia sea ice is approximately 15 cm thicker than 10 the value obtained from the SAR/HIGHTSI data.Sea ice thickness is similar in both cases also in the Gulf of Riga.EarthSyst.Dynam.Discuss., doi:10.5194/esd-2017-35,2017   Manuscript under review for journal Earth Syst.Dynam.Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 11 .
Figure 11. .Modeled (left) and measured (SAR) mean level ice thickness for the most of the Baltic Sea area (06-07.01.2003, color scales represent ice thickness in cm).

Figure 12 .
Figure 12.Mean level ice thickness generation process based on the data from SAR (three right images) and modeled ice thickness for the same day (07.03.2004, color scale represents ice thickness in cm).
Earth Syst.Dynam.Discuss., doi:10.5194/esd-2017-35,2017   Manuscript under review for journal Earth Syst.Dynam.Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.parameters averaged in the selected boxes, representative for different parts of the Baltic Sea.Locations of the rectangles are shown on Fig. 13.

Figure 13 .Figure 14 .
Figure 13.Locations of rectangles taken into account for validation of the ice cover.
the beginning of November sea ice appears in the northern part of the Bothnian Bay.Then, the freezing process spreads 20 Earth Syst.Dynam.Discuss., doi:10.5194/esd-2017-35,2017 Manuscript under review for journal Earth Syst.Dynam.Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.southwards.The range of values and temporal evolution of sea ice frequency and average ice thickness are similar in all boxes.

Figure 15 .
Figure15.A) Maximum annual extent of ice cover in the Baltic Sea since 1720Seinä & Palosuo (1996);Seinä et al. (2001Seinä et al. ( , 2006)).B) Maximum annual sea ice extent in the Baltic Sea in the last three decades from observations and B-CESM (and correlation between modeled and observational data as an insert).
February-March (JFM)) has been increasing since around 1970.The influence of runoff on the ice extent is much lower then heat budget, but excessive runoff decreases sea surface salinity.It is well known that freezing temperature depends on salinity and is going down when salinity increases (the lowest freezing temperature of the oceans is about −1.8 o C. Additional fresh water amount can decrease salinity which is equivalent to increasing freezing temperature and as a consequence it could accelerate growth of the ice extent 4 Summary High resolution coupled ice-ocean model has been successfully developed and implemented for the small regional domain of the Baltic Sea area.The main modification has been done in the ocean component which represents the lower boundary layer of the ice component.The model is based on a new open-source product -CESM.Both components (ice and ocean) were validated and the validation confirms quite good agreement with the in-situ measurements.Furthermore, evidence confirming Earth Syst.Dynam.Discuss., doi:10.5194/esd-2017-35,2017 Manuscript under review for journal Earth Syst.Dynam.Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure B4 .Figure B5 .Figure B6 .Figure B7 .Figure B8 .Figure B9 .
Figure B4.Flow rate vs. sea level difference in the Sund (red points -modeled, black line -calculated from the measured sea level difference).The insert shows position of two stations at the northern and southern end of Sund and the section.

Figure B10 .
Figure B10.Comparison of the measured (FMI, DMI) and modeled ice concentration and thickness for the Bothnia Bay (left) and the Gulf of Finland (right).

Figure B11 .
Figure B11. .Modeled (left) and measured (SAR) mean level ice thickness for the most of the Baltic Sea area (06-07.01.2003, color scales represent ice thickness in cm).

Figure B12 .Figure B13 .Figure B14 .
Figure B12.Mean level ice thickness generation process based on the data from SAR (three right images) and modeled ice thickness for the same day (07.03.2004, color scale represents ice thickness in cm).

Figure B15 .
Figure B15.A) Maximum annual extent of ice cover in the Baltic Sea since 1720 (Seinä & Palosuo (1996); Seinä et al. (2001); Seinä et al. (2006)).B) Maximum annual sea ice extent in the Baltic Sea in the last three decades from observations and B-CESM (and correlation between modeled and observational data as an insert).