MOMBA 1 . 1 – a high-resolution Baltic Sea configuration of GFDL ’ s Modular Ocean Model

We present a new coupled ocean-circulation–ice model configuration of the Baltic Sea. The model features, contrary to most existing configurations, a high horizontal resolution of≈ 1 nautical mile ( ≈ 1.85 km), which is eddyresolving over much of the domain. The vertical discretisation comprises a total of 47 vertical levels. Results from a 1987 to 1999 hindcast simulation show that the model’s fidelity is competitive. As suggested by a comparison with sea surface temperatures observed from space, this applies especially to near-surface processes. Hence, the configuration is well suited to serve as a nucleus of a fully fledged coupled ocean-circulation–biogeochemical model (which is yet to be developed). A caveat is that the model fails to reproduce major inflow events. We trace this back to spurious vertical circulation patterns at the sills which may well be endemic to high-resolution models based on geopotential coordinates. Further, we present indications that – so far neglected – eddy/wind effects exert significant control on windinduced upand downwelling.


Introduction
The Baltic is a shallow sea with an average depth of less than 60 m, situated in central northern Europe between 53 and 67 • N and 10 and 30 • E. It is a brackish sea where intermittent salt water inflows from the North Sea mix with river runoff and precipitation which, on average, exceeds evaporation.In total, 85 million people distributed over 14 countries live in its drainage basin (e.g.Leppäranta and Myrberg, 2009), and there is growing concern about their effect on the pelagic ecosystem (e.g.Neumann et al., 2012).
A major problem is anthropogenic eutrophication (e.g.Pawlak et al., 2009;Wasmund and Uhlig, 2003).Nutrients, such as phosphate and nitrate, stemming from, for example, sewage and fertiliser usage in agriculture drive enhanced primary production by phytoplankton at the sun-lit surface of the Baltic Sea.The associated export of organic matter to depth drives increased oxygen consumption (e.g.Karlson et al., 2002).The implications of the resulting decreased oxygen concentrations at depth are manyfold.One major problem is that decreasing oxygen reduces the volume of seawater which hosts the specific combination of density and dissolved oxygen that is so indispensable for cod eggs to survive.This hampers the recruitment of cod (e.g.MacKenzie et al., 2000), which is an ecologically and economically important fish species in the Baltic.
An additional problem is that depleted oxygen concentrations in bottom waters can trigger fundamental changes in organic matter cycling when certain thresholds are reached.Among the changes, two in particular are thought to set the scene for a dominance of cyanobacteria: (1) denitrification and (2) the release of dissolved inorganic phosphorous (bound to iron oxides) from the sediments to the water column.This is of concern because cyanobacteria can release toxins and thereby constrain recreational use of the water body.
It is considered possible to reduce the nutrient load entering the Baltic Sea so that its assets (e.g.fisheries and tourism) are sustained for generations to come.But, such a reduction calls for either increased cost for sewage treatment facilities or reduced profits from crops and cattle breeding.Because these are large-scale economic decisions affecting highly industrialised countries, the development of projection tools appears mandatory.Further, because any decision causing significant impact must be taken in concert with all the countries in the drainage basin (not all of which actually border on the sea), the development of, and access to, these tools should be as transparent as possible.
This paper describes a cutting-edge Baltic Sea configuration of an ocean circulation model applied to simulate a hindcast for the period 1987 to 1999.With a competitive horizontal resolution of 1 nautical miles (1 nmi = 1.85 km), it is eddypermitting in most of the Baltic (Fennel et al., 1991).We envision to develop, based on this nucleus, a fully fledged coupled ice-ocean-circulation-biogeochemical model.The biogeochemical module will be added in forthcoming projects -step by step as advised by Arhonditsis and Brett (2004) (and bearing in mind the harsh criticism by Flynn, 2005;Anderson, 2005), and starting with simple approaches similar to the one described in Lehmann et al. (2014b).The underlying framework is the MOM4p1 release of NOAA's Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model (Griffies, 2009).The code was chosen because of its exceptional documentation, maturity, modularity and user community.
In this study we aim to demonstrate that our configuration of the coupled ice-ocean-circulation model is competitive by comparing our hindcast with observations.Further, we give guidance on the research questions which could be addressed with the model in its present state.The configuration is distributed, free of charge, by the Biogeochemical Modelling group at the Helmholtz Centre for Ocean Research, Kiel, Germany (http://www.geomar.de/en/research/fb2/fb2-bm/).
The project is dubbed MOMBA, a configuration of GFDL's "Modular Ocean Model version 4p1 for the BAltic Sea".It is accessible via http://www.baltic-ocean.org.

Model setup
MOMBA 1.1 is an eddy-rich configuration of a Baltic Sea coupled ocean-circulation-ice model.The model domain comprises part of the North Sea.At the western edge of the North Sea we restore temperature, salinity and sea surface height to prescribed values.There is no data assimilation in the Danish Straits (i.e.neither any restoring nor prescription of transports).There are no open boundaries and there is no tidal forcing.Air-sea fluxes in the entire domain are based on results from an atmospheric reanalysis.The remainder of this section provides the details.

Grid and bathymetry
The model domain covers the area from 4.2 to 30.3 • E and 53.8 to 66.0 • N (Fig. 1a).There are no open boundaries.The North Sea is represented by a simplified basin with solid walls in the north, south and west.The horizontal discretisation is an Arakawa B-grid (Arakawa and Lamp, 1977) and comprises 603 × 588 tracer grid boxes.The longitudinal resolution is a function of longitude, starting with 0.18 • at the western boundary, increasing to 0.0375 • at 9.36 • E and decreasing slightly thereafter to 0.03998 at the eastern boundary.The meridional resolution is 0.02 • .Expressed as the harmonic mean of longitudinal and meridional resolution, this corresponds to a horizontal resolution of 0.9 to 1.3 nautical miles in the Baltic Sea, and up to 2 nmi in the North Sea.This resolution is exceptionally high and resolves almost the full range of first baroclinic Rossby radii in the Baltic Sea (e.g.Fennel et al., 1991).
The vertical discretisation comprises a total of 47 levels.Figure 2 shows the nominal depth and thickness of each level.We use the zstar (z ) coordinate (Stacey et al., 1995;Adcroft and Campin, 2004) in the vertical, i.e. the depth and thickness of each level varies with time.z is calculated as a function of nominal depth (z), water depth (H ) and sea surface height (η) (which varies with time because we use a free surface): (Eq. 6.6 in Griffies, 2009).The approach overcomes the problem with vanishing surface grid boxes which appears in generic z-level discretisation when sea surface height variations are of similar magnitude to the thickness of the uppermost grid box.Hence, the z approach allows for higher vertical resolution at the surface than the traditional z-level coordinates.
The model bathymetry is interpolated from the highresolution database compiled by Seifert et al. (2001) (http: //www.io-warnemuende.de/topografie-der-ostsee.html).We use partial cells.The interpolation scheme is bilinear and takes the respective 10 nearest neighbours into account.Water columns shallower than three grid boxes are deepened to three grid boxes in order to ensure numerical stability.The model bathymetry is smoothed with a filter similar to the Shapiro filter (Shapiro, 1970).The filter weights are 0.25, 0.5 and 0.25.The filtering procedure can only decrease the bottom depth, i.e. it essentially fills rough holes.The filter is applied three times consecutively.The resulting bathymetry contained lakes which we filled after visual inspection.In addition, we filled narrow inlets which received large river runoffs and thus developed numerical oscillations (river mouths were shifted accordingly thereafter).
In total, MOMBA has 1 483 236 wet tracer grid boxes.This is still insufficient to fully resolve the complex bathymetry and coastline (Fig. 1b-d), especially in the narrow and shallow Danish Straits (Øresund, Great Belt, Little Belt) which control the exchange of water between the North Sea and the Baltic Sea (Fig. 1b).Aggravating in that respect is that the horizontal Arakawa B-grid discretisation necessitates a "no-slip" boundary condition for velocity grid points bordering the bottom.We suspect that this dampens our simulated transport through the Danish Straits.In order to compensate for this spurious effect, we deepen the bathymetry at several grid points in the Danish Straits as documented in Table 1.

Initial conditions and sponges
Situated in central northern Europe, the Baltic Sea is considered to be among the best-sampled regions worldwide.
Even so, because of its high temporal and spacial variability it is not straightforward to derive initial conditions for a highresolution ocean circulation model which features more than 1 million grid boxes.
We use temperature and salinity bottle data archived by the International Council for the Exploration of the Sea (hereafter ICES, http://ocean.ices.dk).The data were extracted upon request by E. J. Green at ICES in February 2009.All data, from 1961 to 2008, were binned into monthly boxes of 6 nmi horizontal resolution with a vertical thickness corresponding to the model grid (Sect.2.1).These monthly binned data were subsequently averaged over all months of the whole 1961 to 2008 period.As a third step we applied an objective mapping with a horizontal decay length scale of 10 nmi in order to smooth and extrapolate the data.Our approach results in a smooth three-dimensional, plausible temperature and salinity distribution throughout the whole Baltic Sea.In order to identify a date at which these smooth temperature and salinity fields resemble actual conditions, we compared the smoothed fields visually to observations of nearbottom salinity (as provided within the SHARK database (http://produkter.smhi.se/pshark/)at SMHI (Swedish Meteorological and Hydrological Institute, Norrköping, Sweden)).This comparison, based on the monitoring stations BY5, BY15 and BY31 (cf.Fig. 1a), suggests that our smoothed temperature and salinity fields are similar to actual conditions on 1 January 1987.Hence, we start our simulation with the smoothed fields on 1 January 1987.
West of 9.3 • E, temperatures (and salinities) are restored to prescribed temperatures (salinities) using an exponential law with a restoring timescale of 1 year (1 day).Salinity is restored to 35 PSU throughout the water column.Temperature is restored to the temperature field described above (i.e. to the initial conditions) below 80 m.Above 80 m we do not apply temperature restoring such that a seasonal cycle can freely evolve.

Sea surface height restoring
On average, the Baltic Sea exports fresh water at the surface to the North Sea, and receives deep, saline inflows, from the North Sea.The inflow events are intermittent and show a strong correlation with atmospheric circulation patterns (e.g.Schinke and Matthäus, 1998;Wyrtki, 1954) because the winds drive differences in sea surface height between the Baltic Sea and the North Sea.Their associated pressure gradients, in turn, drive the exchange between the Baltic Sea and the North Sea.
As described in Sect.2.1, the North Sea is represented by a simplified basin.This retards the simulation of sea surface height (SSH) in the North Sea, which, in turn, leads to an unrealistic exchange between the Baltic Sea and the North Sea.In order to compensate for this we restore sea surface height between the western edge of the model domain and 6.895 • E (the latter longitude refers to the centre of the tracer grid box).The respective restoring timescale is prescribed as a linear function of longitude, starting with 1 h at the western boundary and increasing to 12 h at 6.895 • E.
We do not restore directly to observations of sea surface height because this would limit the application of the model: coupled ocean-atmosphere-Earth system models are not (yet) up to the task of providing reliable sea surface height projections at the entrance to the Baltic Sea.Hence, if we were to restore directly to sea surface height, we would be constrained to the past where we have observations.Instead, we restore to SSH derived from a proxy which can reliably be simulated.Specifically we restore to SSH as calculated by a regression to the atmospheric sea level pressure gradient between Oslo in Norway and Szczecin in Poland.The ratio behind the approach is the high correlation between the normalised atmospheric pressure gradient, dubbed the Baltic Sea index (BSI), and the exchange between the Baltic and the North Sea observed by Lehmann et al. (2002).
We restore to SSH (η), calculated as The respective coefficients are from a regression of monthly mean Permanent Service for Mean Sea Level (PSMSL) (PSML, 2012) tide gauge data against the BSI (A.Lehmann, personal communication, 2013).The BSI used here is identical to the one presented in Lehmann et al. (2002).

River runoff
The river runoff applied in this study is a crossover from observations (Bergström and Carlsson, 1994) and a large-scale hydrological model (Graham, 2004).It is identical to the one applied in Meier and Faxen (2002), Meier and Kauker (2003), Hordoir and Meier (2010) and Löptien and Meier (2011).Figure 3 shows the total discharge to the model domain.We resolve 30 discharge locations (rivers and groundwater inflows).In high-resolution Baltic Sea configurations, river runoff is prone to cause numerical problems because the salinity is close to zero in the vicinity of some rivers.This can, at times, due to the dispersive nature inherent to all state-ofthe-art advection schemes (other than the outdated upstream scheme), induce negative salinities for which the equation of state is not defined.We solved this problem by distributing the runoff over the six grid boxes closest to the actual position of the river mouth.The distribution among the six grid boxes is determined by weights (W i ), calculated with an exponential decay length scale τ = 3 km as a function of the distance of respective grid boxes to the actual river mouth (D i ).

Air-sea fluxes
The atmospheric boundary conditions driving our ocean-ice model are identical to those applied to the "RCO" oceanice model in Meier and Faxen (2002), Meier and Kauker (2003), Hordoir and Meier (2010) and Löptien and Meier (2011).They are based on a downscaled ERA-40 reanalysis (Uppala et al., 2005, and many more to follow) which, compared to the original ERA-40 reanalysis, features additional regional details which considerably affect the solution of stand-alone ocean models of the Baltic Sea (Meier et al., 2011).The downscaling procedure takes ERA-40 reanalysis data as boundary conditions for the Rossby Centre Regional Atmosphere model version 3 (hereafter RCA3).RCA3 features an enhanced (relative to ERA-40) horizontal resolution of 25 km (Jones et al., 2004;Samuelsson et al., 2011).As shown by Samuelsson et al. (2011), the approach provides a very realistic climatology.A shortcoming, however, is that RCA3 does not employ techniques such as spectral nudging (von Storch et al., 2000) which constrain simulated cyclone paths to match actual tracks.This is of some concern because a systematic analysis of cyclone tracks (e.g.similar to analysis carried out by Löptien et al., 2008) simulated with RCA3 has not been carried out yet.We use cloud cover, precipitation, near-surface humidity, air pressure, temperature and wind speed from the RCA3 model at a temporal resolution of 3 h as provided by the Rossby Centre at SMHI, Norrköping, Sweden.As described in the following, the cloud cover is used to calculate downward longwave and net shortwave radiation, while the other RCA3 variables, combined with the simulated sea surface properties, are used to derive actual air-sea fluxes following Large and Yeager (2004) (their Sect.2.1): Solar net shortwave radiation entering the ocean is calculated as described by Meier et al. (1999) and resolves daily cycles.Because the Meier et al. (1999) report is not available in electronic form, we repeat (keeping their notation) in the following their calculation of the net solar shortwave radiation (Q SW ): where Q 0 SW is incoming radiation and α w is sea surface albedo.The incoming radiation is with atmospheric turbidity T u = 0.95, solar constant S 0 = 1353 W m −2 , zenith angle θ, transmission function T r , absorption function A w , cloud function F a , and cloud cover C a .We define θ is a function of latitude φ, the Sun's declination angle δ and the Sun's hour angle α n : where d is the day of the respective year.
with t s the hour of the respective day.m in Eq. ( 8) is the optical path length, calculated as The albedo α w is calculated with Fresnel's formula: where the refraction angle ψ = arcsin (sin θ/1.333).
Q SW is absorbed in the uppermost wet grid box, which has a nominal thickness of 3 m (Sect.2.1).This is a pragmatic decision because a prognostic model of water turbidity is yet to be developed.Note that this induces an error since some of the incoming radiation reaches below the surface mixed layer depth under certain conditions (e.g.Löptien and Meier, 2011).
Longwave radiation entering the ocean (LWRAD down ) is calculated from cloud cover as in the RCO model used by Meier and Faxen (2002), Meier and Kauker (2003), Hordoir and Meier (2010) and Löptien and Meier (2011).The formulation is unusual since it contains neither the water vapour pressure nor the dew point temperature.However, given the high uncertainty of these empirical relationships in general (e.g.Zapadka et al., 2007) we apply, for legacy reasons, the following calculation: where σ = 5.6703 × 10 −8 W m −2 K 4 is the Stefan-Boltzmann constant, T air is the air temperature at 2 m height, and CLCO is the cloud cover expressed as a number between 0 and 1.
Longwave radiation leaving the ocean is calculated on-line by the model based on simulated sea surface temperatures following Large and Yeager (2004).We assume an emissivity of 1.
The turbulent air-sea fluxes -i.e. the wind stress, the evaporation, the sensible and the latent heat flux -are calculated with bulk formulas according to Eqs. (4a)-(4d) in Large and Yeager (2004), respectively.We account for ocean currents in the bulk parameterisation of the wind stress.As explained in the forthcoming Sect.3.2, this consideration of ocean currents significantly affects the model solutions.
The calculation of ocean-atmosphere fluxes which are actually applied to the ocean every time step is updated (based on linearly interpolated 3-hourly atmospheric boundary conditions described above and the contemporary ocean state) every 36 min, and kept constant in between.
Additionally, we account for the effect of atmospheric bottom pressure gradients on ocean model pressure gradients.

Results
To serve as a nucleus of a fully coupled biogeochemical circulation model of the Baltic Sea, the ocean circulation module has to perform with a certain level of quality.To this end we present a choice of model assessments: -Surface circulation (Sects.3.1 and 3.2), which, for example, affects the wind-driven upwelling of nutrients (as we will put forward) and subsequent pathways of nutrients.
-Sea surface temperature (SST, Sect.3.3), which can be observed from space and is therefore available in an www.geosci-model-dev.net/7/1713/2014/unrivalled (compared to in situ measured properties) spacial and temporal resolution.SST serves as a proxy for the realism of the surface mixed layer depth, whose dynamics are a major process involved in the supply of nutrients from depth to the sun-lit surface ocean.
-Inflow of saline, oxygen-rich waters from the North Sea to the Baltic and its subsequent spread at depth (Sect.3.4) is a major agent resupplying the Baltic Sea with oxygen at depth.This is of importance because it determines whether organic matter sinking to depth is denitrified rather than remineralised under oxic conditions.
-Sea ice (Sect.3.6), which caps the ocean from the atmosphere throughout almost half of the time up north and thereby controls the air-sea gas exchange of climaterelevant species such as carbon and N 2 O.

Average surface circulation
The circulation in the Baltic Sea is driven by the wind exerting stress at the sea surface, sea surface tilt effected by, for example, river runoff, thermohaline horizontal gradients of density and tidal forces.Overall, this results in a predominantly wind-driven, weak, cyclonic circulation with surface speeds of some 5 cm s −1 which are locally quite persistent (Leppäranta and Myrberg, 2009;Sjöberg, 1992).To this end our simulation is consistent with the observations, as is suggested by Fig. 4. Unfortunately a more quantitative model assessment is hindered by the lack of current measurements available to us, and thus we are restricted to a comparison with sea surface salinity and other model simulations.
Sea surface salinity (SSS) is determined by the surface circulation mixing river runoff and salty water of North Sea origin.Hence, as put forward by Eilola and Stigebrand (1998) and exploited by Lehmann et al. (2002), a realistic SSS is a proxy for a realistic surface circulation.Figure 5 shows a comparison between simulated and observed SSS which suggests that the simulated surface circulation is realistic.
As concerns other simulations, we find that patterns of the magnitudes of mean surface velocities are comparable, although weaker in amplitude, with recent results from Jedrasik et al. (2008) based on a model that covers the entire Baltic Sea.The most pronounced features in our simulation which are different to Jedrasik et al. (2008) are (1) lower velocities in the Archipelago Sea east of Åland, (2) a distinct band of high velocities in the northeastern Bothnian Sea  west of Ulvö Deep and (3) a distinct cyclonic gyre southeast of Gotland.
1.As regards the Archipelago Sea, it is safe to assume that both models are deficient because their spacial resolution does not allow for the multitude of islands east of Åland to be resolved (Fig. 1d).
2. Regarding Ulvö Deep, our model is consistent with the simulation of the Gulf of Bothnia by Myrberg and Andrejev (2006), which also features high southwards velocities.Furthermore, these authors report, consistent with our results, a high persistency, R, of these currents.
The persistency R is defined as the ratio between vector and scalar mean speeds: where vector U is mean vector speed and U the mean scalar speed.Figure 4 shows a map of simulated R: in summary and consistent with results from Jedrasik et al. (2008), Lehmann et al. (2002), Myrberg and Andrejev (2006), and Andrejev et al. (2004), we find currents of a persistency exceeding 80 % along closely packed isobaths (Fig. 4) -a feature that can be described in terms of a depth-integrated vorticity balance (e.g.Lehmann et al., 2002).

Transient surface circulation
Even though the mean circulation can provide insight into tracer dynamics (such as, for instance, transports of salinity, pollutants and nutrients from river mouths), the concept is limited by the highly transient nature of the Baltic Sea circulation.Figure 6 shows a measure of this transiency, the  simulated eddy kinetic energy (E) defined as where u and v are deviations of the zonal and meridional surface velocities from their time mean over the 1987 to 1999 period, respectively.The overbar denotes a temporal average over the same period, and u and v are calculated from daily means.Typical eddy kinetic energies are around 50 cm 2 s −2 , with hotspots exceeding 200 cm 2 s −2 .For comparison, this is an order of magnitude less than the up to 2000 cm 2 s −2 found, for example, off Japan, where the Kuroshio separates from the coast; however it is still substantial when compared with typical open-ocean conditions away from world's most powerful boundary currents (Dietze and Kriest, 2012;Liu et al., 2010).
A comparison between persistency (Fig. 4) and eddy kinetic energy (Fig. 6) suggests that the term eddy kinetic energy is somewhat misleading because some regions of high eddy kinetic energy are characterised by high current persistency (e.g.west of Ulvö Deep).This is apparently inconsistent with the generic concept of whirling eddies.The solution to this conundrum is that eddy kinetic energy (Eq.16) can be effected by intermittent unidirectional currents which score a high persistency (Eq.16).Hence, our model predicts that some of the regions which host the most powerful surface current variability are characterised by intermittent but unidirectional currents which are steered along isobaths (Figs. 6  and 4).
As concerns the small-scale transient eddying circulation we report, in the remainder of this subsection, on a subtletythe so-called eddy/wind effects (Martin and Richards, 2001): our model configuration explicitly includes the effect of the  16) in units of cm 2 s −2 .ocean surface circulation on the stress exerted by the winds on the sea surface.To this end our configuration is more realistic than those studies which neglect this effect (e.g.Jedrasik et al., 2008, to name one of the very few studies which document this aspect of their wind forcing).The eddy/wind effect introduces a non-linearity into the dynamics, since the wind stress is a function of the difference between the wind speed and the movement of water at the surface, which in turn is determined by the wind stress.Recent open-ocean studies have suggested that this non-linear effect is significant (e.g.Duhaut and Straub, 2006;Zhai and Greatbatch, 2007;Xu and Scott, 2008) even though surface currents are typically an order of magnitude smaller than the winds.The reason is that the spacial scales of transient circulation are much smaller in the ocean than in the atmosphere (which is a consequence of the stronger stratification in the atmosphere resulting in larger Rossby radii).Thus, even though small in magnitude, ocean currents can significantly alter the curl of the wind stress (which is a function of spacial derivatives of stress) exerted on the sea surface.As put forward by Martin and Richards ( 2001 etze, 2009).Note that the Ekman pumping features variability correlated with the surface currents on spacial scales much smaller than those that could possibly be effected by the ≈ 25 km resolution wind fields.Ledwell et al., 2008;Eden and Dietze, 2009).It is straightforward to assume that this holds even more so in the Baltic Sea, where shallow water depths reduce the Rossby radii to a couple of kilometres (compared to the several tens of kilometres typical for open-ocean conditions).Figure 7 supports this hypothesis: Ekman upwelling (downwelling), w e , calculated by taking into account the surface velocity in the calculation of the wind stress, is correlated with small-scale anticyclonic (cyclonic) features.Since these features are smaller than the resolution of the atmospheric forcing, they have to be significantly influenced by the surface circulation modulating the wind stress exerted on the ocean's surface.Further, we find, consistent with the Rossby radii being an order of magnitude smaller, that w e is of an amplitude of several metres per day -almost an order of magnitude larger than the ≈ 0.5 m day −1 observed in the open ocean by, for example, Ledwell et al. (2008).

Sea surface temperature
We use the KPP boundary layer parameterisation (Large et al., 1994) with the same parameters applied in eddypermitting global configurations of Dietze and Kriest (2012), Dietze and Löptien (2013) and Liu et al. (2010).Remarkably, the simulated SSTs are very realistic, which is a real asset of MOMBA. Figure 8 shows a comparison with SHARK (Swedish Oceanographic Data Centre at SMHI) data which suggests that the only flaw of the model is, at times, an underestimation of maximum summer temperatures.We speculate that this is the consequence of the following: 1.The vertical discretisation assigns a nominal (cf.Eq. 1) thickness of 3 m to the uppermost grid box (Fig. 2).This (implicitly) imposes a lower bound of 3 m on our simulated surface mixed layer depths which can be unrealistically high during low-wind summer conditions.
2. The neglected effects of phytoplankton blooms, which can shift the absorption of solar radiation in the water column upwards, which in turn increases the SST (e.g.Löptien et al., 2009;Löptien and Meier, 2011).
Observations from space confirm the extraordinary quality of our SST hindcast: Fig. 9  MetOp satellites distributed by the Bundesamt für Seeschifffahrt und Hydrographie, Germany, via http://www.bsh.de/de/Meeresdaten/Beobachtungen/Fernerkundung.The data are referred to as BSH SST data in the following.For the following model-data comparison we set all BSH SST data below −0.5 to −0.5 • C in order to prevent accidently comparing simulated SSTs with observed sea ice temperatures.Typically, we find that the model explains more than 90 % of the variance in the BSH SST data.Exceptions with lower correlations (which will be discussed in the following) are the (1) Skagerrak, (2) the Bothnian Bay, (3) the Archipelago Sea east of Åland and (4) those regions hosting frequent upwelling (Lehmann et al., 2012) such as the Swedish south and east coasts, south of Gotland and off Finland in the Gulf of Finland: 1. We speculate that simulated temperatures in the Skagerrak are deteriorated by restoring to an annual mean climatology (Sect.2.2).
2. The reason for the low correlations in the Bothnian Bay is probably associated with sea ice melting too late in spring (Sect.3.6).
3. The low correlation in the Archipelago Sea is not surprising given the unrealistic land-water distribution owing to scant resolution (Fig. 1d).
4. The low correlations in the upwelling regions are not necessarily indicative of a poor representation of associated processes.In fact, visual inspection of SST maps such as those shown in Fig. 10 reveals striking similarities between the model and the BSH SST.However, on a grid point to grid point basis there is a mismatch.First of all, this is because we compare weekly composites (assembled during cloud-free periods) with simulated daily averages of upwelling events, some of which have a lifetime of a few days only.Second, upwelling is highly non-linear (also because of the subtle circulation/wind effects described at the end of Sect.3.2).This non-linearity makes it difficult (or even impossible) to retrace the dynamics of single upwelling events on a one-to-one basis.

Deepwater renewal
The assimilation of inorganic carbon by algae in the sun-lit surface ocean and its subsequent export as sinking particulate organic carbon fuels oxygen consumption at depth.This consumption is balanced by resupply of oxygen-rich waters of North Sea origin which enter the Baltic intermittently via the Danish Straits.When this resupply falls short of the consumption, the oxygen concentrations can be depleted below certain thresholds.In turn, denitrification and the release of dissolved organic phosphorous bound to iron oxides from the sediments may set in.This is (although discussed controversially; Pahlow et al., 2013) of concern, because the associated change in the N / P ratio is thought to promote the growth of potentially harmful cyanobacteria (e.g.Vahtera et al., 2007).
Unfortunately, MOMBA features a spuriously hampered deepwater renewal: Fig. 11 shows a comparison between simulated and observed near-bottom salinities at stations BY2, BY5 and BY15.The further into the Baltic, away from the Danish Straits, the weaker the simulated salinity signature of the inflow events.This holds especially for the massive 1993 event, which is not reproduced at all at station BY15.
The underestimation of deep salinities is an infamous feature of z-level models which do not apply data assimilation techniques (e.g.Fu et al., 2012).However, the especially bad performance of MOMBA came as a surprise because of the model's the high and competitive spacial resolution and its -as suggested in the following sea surface height analysisrealistic exchange of water through the Danish Straits.
Sea surface height variations are effected by, for example, (1) variations in air pressure, wind and density in the Baltic; (2) internal eigenoscillations and standing waves (so-called seiches); and (3) "external" forcing by varying sea level outside the Baltic (and, to a lesser extent, by freshwater supply to the Baltic).The relative contribution of these mechanisms depends on the actual location and the frequencies considered, also because the seiches have distinct spacial characteristics both in space and in the frequency domain.To this end, the area around Stockholm/Landsort is special, because it is situated at the node of major seiches, and therefore sea level readings are rather unaffected by internal oscillations.Samuelsson and Stigebrandt (1996) report that the variability at the Stockholm/Landsort gauges is highly correlated with "external" forcing.For this reason, a realistic simulation of the sea surface height at Stockholm/Landsort is generally understood as proof of a realistic water exchange through the Danish Straits (e.g.Neumann and Schernewski, 2008).
Figure 12 shows the simulated daily mean sea surface height along with daily mean observations (which comprises the world's longest continued series of sea level observations; Ekman, 1988;Woodworth and Player, 2003;PSML, 2012).Specifically we use data retrieved on 5 June 2013 from the Permanent Service for Mean Sea Level at http: //www.psmsl.org/data/obtaining/(Holgate et al., 2013).Expressed in terms of a correlation, we reach 0.78 and 0.85 at Stockholm and Landsort, respectively.The number of observations is 4748 in both cases.Compared to previous studies, which typically report correlations higher than 0.9 (e.g.Neumann and Schernewski, 2008;Jedrasik et al., 2008), this is on the lower end of published models, but not extraordinarily poor.To this end it is worth noting that our rather low correlations are not necessarily indicative of poor model dynamics but may well be preliminarily caused by our boundary conditions: in contrast to other studies, we do not restore to observed sea levels in the Kattegat/Skagerrak but prescribe sea surface height derived from a correlation with a normalised atmospheric pressure gradient (Sect.2.3) further west, at the edge of our model domain.Our approach is adopted from Lehmann et al. (2002), who report a weekly (monthly) correlation of 0.83 (0.82) of their simulation with weekly (monthly) sea level data from Landsort.Compared to their performance, we reach a competitive 0.86 (0.88) calculated based on weekly (monthly) data at Landsort.In line with the results from the sea surface height analysis, a more direct comparison with an observational estimate from the 1993 inflow in Table 2 confirms that the simulated transport through the Danish Straits is indeed realistic.As concerns the associated salinity transport, it is safe to assume that it is not biased low because we restore salinity to a rather high (and unrealistic) 35 PSU in the North Sea (Sect.2.2).Overall, the combination of a rather high salinity transport to the Baltic and its missing salinity signatures at depth leaves us with the question as to where the salty and dense inflow gets lost on its way to the Baltic Proper.It is straightforward to assume that the problem of an apparently hampered deepwater renewal is associated with our choice of the bottom boundary layer parameterisation (Döscher and Beckmann, 2000).However, excessive tuning exercises (not shown) suggest that this is not the case.This is probably not so surprising, because Fig. 13a suggests that most of the topographic slopes are actually well resolved, which in turn should prevent much of the infamous mixing at sills in Boussinesq-approximated geopotential coordinate models.
Analysis based on temporal high-resolution output finally revealed the cause of the spurious dilution of salty inflows and their subsequent absence in the deeper basins: Fig. 14 shows a section along the Hamrare Strait (as indicated in the inlayed map), which is typical for the situation during strong inflows at the sills in MOMBA.On timescales of hours and spacial scales of several nautical miles, vertical circulation cells with alternating up-and downward velocities develop during stronger inflows.The magnitude of these circulation patterns is of the order of several tens of metres per day and they intermittently cross the surface mixed layer.Thus, the upward branches of these patterns transport deep saline bottom water to the surface (and, likewise, lighter surface waters to the bottom), where it is spread horizontally by the prevailing turbulent circulation (Fig. 6).This can cause vertical instabilities as dense bottom water is shifted horizontally over lighter subsurface waters, which in turn drives further vertical mixing.The overall effect is one of enhanced vertical mixing at the sills during strong inflow events.The salt associated with the inflow events is mixed upwards to the surface and eventually flushed out of the Baltic by the net positive freshwater supply rather than reaching the deeper basins further away from the Danish Straits.
We speculate that the vertical circulation cells during strong inflows at the sills (as the one shown in Fig. 14 for the Hamrare Strait) are numerical artefacts caused by the divergence of horizontal currents interacting with the steplike representation of topography.At the very least, even if effects like that exist (although we have found no supporting evidence in the literature), their mixing effect is obviously unrealistically strong in our model.Further, we speculate that the problem is more pronounced in high-resolution configurations like ours, because the spacial scale of up-and downwards jets is several grid boxes only (Fig. 14).Previous Boussinesq-approximated configurations on geopotential coordinates with resolutions of several nautical miles might just have been fortunate enough to not "resolve" these spurious cells.
The remaining question is how to dampen theses cells.To this end we performed an extensive set of tuning experiments including modified winds, surface mixed layer depth and bottom boundary layer parameterisations, and initial conditions, all branching off the reference simulation on 1 January 1993.As regards improved realism, all were to no avail.Figure 15 shows the most promising approach:  increased viscosity dampens the spurious circulation cells.Hence, higher viscosity reduces the spurious mixing and comes along with increased simulated near-bottom salinities at station BY5.When the viscosities get too high, however, the saline inflows are decelerated, and even less salt finds its way to the deep Baltic.Note that the apparently optimal solution with an added vertical viscosity of 10 −2 m 2 s −1 , which reaches realistic values in January 1993 at station BY5 in Fig. 15b, features an overall worse performance than the reference simulation.Similar behaviour in fact applied to all our tuning experiments: while quite a few of them reduced local misfits to observations, none of them achieved global improvement.
We conclude that spuriously damped deepwater renewal is endemic to MOMBA.The underlying mechanism is apparently related to the high spacial resolution, a finding which is consistent with results from Ezer and Mellor (2004) and tuning exercises with the BaltiX model configuration (Hordoir et al., 2013): R. Hordoir (personal communication, 2012) reports that increasing viscosity at the surface only increases near-bottom salinities in the deep Baltic basins up to realistic values, leading to an improved global performance.We speculate that this is the case because the higher near-surface viscosity dampens that part of the spurious circulation pattern which effectively mixes across the surface mixed layer and, at the same time, does not slow down the near-bottom overflows which transport the salt from the Danish Straits far into the Baltic.
Additional and pragmatic remedies other than applying non-homogenous viscosities could be adopted from previous efforts to model the Denmark strait overflow: Kösters et al. (2005), for example, parameterise the transport across the (Denmark Strait) sill as a function of the large-scale conditions and then -directly -impose this transport, thereby overriding the effects of resolved spurious circulation on the sill.It is worth noting that their improved and apparently more realistic model features significantly differing climate projections (Köller et al., 2010).We anticipate a similar behaviour in our case and conclude that a comprehensive study on this subject should preclude any attempts to make projections with the current modelling framework.

Halocline
A permanent halocline which separates the dense saline nearbottom waters of North Sea origin from the surface layer, which comprises fresh river runoff, is a characteristic feature of the Baltic.Its depth is determined by the balance between inflowing saline waters and its subsequent diapycnal upwards transport (via upwelling and mixing processes).Thus, a realistic simulation of the halocline depth is a necessary precondition for simulating vertical nutrient and oxygen fluxes.
The simulated median salinity profile at station AE (Fig. 1a) in the Kattegat is of a shape very similar to the observed one (Fig. 16a).However, there is an offset, probably caused by the rather high (and unrealistic) salinity value to which we restore in the North Sea (35 PSU, Sect.2.2).Further into the Baltic, away from the Danish Straits, for example at station BY5 (Fig. 16b) and BY15 (Fig. 16c), salinities are no longer biased high.In fact, as outlined in Sect.3.4, near-bottom salinities are low because of spurious vertical circulation patterns at the sills during major inflow events.Even so, we conclude, based on Fig. 16 (and corresponding analysis at all other stations indicated in Fig. 1, which we omit here for the sake of brevity), that the transition from the surface mixed layer to the stratified halocline is simulated realistically (i.e.within ≈ 10 m of the observations).This puts confidence in the simulated diapycnal transport mechanisms at the base of the surface mixed layer, which is an important asset of the model since the associated nutrient fluxes drive primary production and subsequent respiration at depth (e.g.Dietze et al., 2004).

Sea ice
The Baltic Sea is seasonally ice-covered.To first-order approximation, this cover acts as a lid, hindering the exchange of heat and moisture between the ocean and the atmosphere, thereby essentially decoupling the systems.Further, it prevents the exchange of climate relevant trace gases such as CO 2 and N 2 O.To this end, sea ice dynamics are of interest, maybe even more so in a warming world of ever-retreating ice cover.Here we compare our model results with digitised ice charts from the operational ice service provided by SMHI.The charts are based on information from ship crews, monitoring stations blended and complemented with measurements from space.The major consumer of this information is commercial shipping, and it remains to be shown that it provides a realistic ground truth to which models can be compared.That said, we assume in the following that the operational ice charts are, to date, the most reliable large-scale source of information for the respective time period.
Figure 17 shows that the simulated ice extend is biased high compared to the observations.This spurious feature is further highlighted by exemplary snapshots for particularly warm and cold winter conditions in mid-February (typically the time of maximum extent; Leppäranta and Myrberg, 2009) in Fig. 18.
Figure 17 suggest that, although biased high, the simulated inter-annual variability of the ice extent is realistic.This is confirmed by an empirical orthogonal function (EOF) analysis of both simulated and observed mean seasonal mean ice concentrations: the leading EOF explains 71 % (68 %) of the variability of observed (simulated) ice concentrations and thus covers the major signal.The corresponding leading principal components (PCs) are almost identical (corre-0.92)and show (Fig. 19), furthermore, a strong correlation (0.8) with the PC-based North Atlantic Oscillation (NAO) index (as applied, for example, in Hurrell, 1995;Hurrell and Deser, 2009;Löptien and Ruprecht, 2005).This is in agreement with published state-of-the-art model studies (e.g.Tinz, 1996;Jevrejeva and Moore, 2001;Jevrejeva et al., 2003;Koslowski and Loewe, 1994;Löptien et al., 2013).
As regards spacial patterns of variability, the similarity between the leading EOFs of the simulation and the observations in Fig. 20 suggests that the underlying dynamics which drive the ice are -overall -realistically reproduced by the model.However, as a consequence of the model's notorious overestimation of ice concentrations, the corresponding EOF patterns are more pronounced in the simulation.This holds particularly in the Bothnian Sea.
Additional model deficiencies are apparent in Fig. 20, which shows spuriously reversed phase relations off Germany and up north in Bothnian Bay and the Bothnian Sea.Off Germany this is the consequence of the simulated bias in ice extent.As for up north, we speculate that the spuriously reversed sign in the EOF pattern is related to unresolved wind-ice interaction: cold winters are often related to persistent northerly winds, which can cause frequent leads (Löptien et al., 2013).The reproduction of these areas of open water is probably hindered by model's spacial resolution (even though we run on a competitive 1 nmi).

Grid-scale noise
Ocean circulation models which use an Arakawa B-grid spacial discretisation (Arakawa and Lamp, 1977) (like MOMBA) are prone to catch numerically induced gridscale noise in the simulated velocity field and sea surface height (Mesinger, 1973;Killworth et al., 1991;Pacanowski and Griffies, 1999;Griffies et al., 2001).These so-called checkerboards are especially pronounced in coarsely represented enclosed embayments (Griffies, 2009) when the applied viscosity is rather small (e.g.Jochum, 2008).Given that (1) the Baltic Sea is almost enclosed, (2) our spacial discretisation is still rather coarse even though we have a resolution of 1 nmi (compare Fig. 1b, c and d) and (3) that we aim to keep the explicit viscosity as low as possible such that we have no unrealistic damping of the circulation, special attention was paid to the checkerboard problem, especially because Löptien and Meier (2011) report checkerboards in simulated temperature in an Arakawa B-grid model of the Baltic Sea with a horizontal resolution of 2 nmi.
We use the dissipative predictor-corrector barotropic time stepping scheme with a dissipation parameter of 0.2.In addition we smooth the sea surface height spacially (1) at each barotropic time step with a Laplacian operator and (2) at each longer, baroclinic time step with a scale-selective biharmonic operator.In both cases (i.e, the Laplacian and the biharmonic operation) the associated parameter is a velocity of 0.5 cm s −1 .The applied Laplacian (biharmonic) diffusivity is then calculated via multiplication by the (cube of the) harmonic mean of the zonal and meridional grid length.Details of the sea surface height filtering applied in MOMBA are described in Sect.8.3 of Griffies (2009).Visual inspection of snapshots of all prognostic variables in MOMBA revealed no evidence of checkerboard whatsoever.

Computational uncertainty
Model integrations are not necessarily reproducible.This applies when switching from one type of hardware to another, changing the compiler, or changing compiler switches.Among the reasons is that round-off errors generally violate commutative laws.Hence, the order of execution -which is tailored by the compiler for optimal performance on the respective hardware -can influence the results.Note that reproducibility is not even guaranteed if the identical executable is executed on the identical hardware (as is the case with aggressive floating-point optimisations set by Intel's fp-model compiler option).
Reproducibility issues with numerical integrations of highly non-linear systems have been reported elsewhere (e.g.Song et al., 2012) and can be problematic: Fig. 21 shows results from two "identical" integrations applying the default (but rather aggressive) Intel fp-model compiler option.The difference in surface properties seems to be insignificant (compare panel a with panel b in Fig. 21).However, the effects of the highly non-linear vertical circulation patterns above the sills during inflow events (Sect.3.4) amplify to a striking difference in near-bottom salinities at station BY5 in early 1993 (panels c and d in Fig. 21).

Computational cost
MOM4 runs on a wide range of machines.For instance, we successfully tested coarser configurations on a eightcore MacPro (Intel Xeon, 2.5 GHz, 12 MB L3 cache, 16 GB RAM) before we switched to higher-performing hardware.
The configuration presented in this study was benchmarked on three different high-performance computing clusters (HPCs), all of them interconnected with InfiniBand.The respective underlying computational nodes were (1) Intel Xeon Harpertown E5472, 8 core; (2) Intel Xeon Gainestown X5570, 8 core; and (3) Intel Xeon Sandybridge, 2 × 8 core.In line with Schmidt (2007), we find that the code is suitable for parallel processing up to several hundred tasks (Fig. 22).Unfortunately, however, we found comparably poor scalability on our main workhorse, the Xeon Sandybridge (RZ-Kiel/GEOMAR): as shown in Fig. 22, performance drops dramatically at around 100 cores, probably as a consequence of an IO bottleneck during initialisation.
In addition to tests on HPCs, we ran preliminary benchmarks (albeit running the MOM5 kernel, which may or may not be comparable) on a contemporary, inexpensive (≈ EUR 10 000) system consisting of two 32-core workstations, based on a 6320 AMD Opteron (Abu Dhabi) (8-core CPU, 2.80 GHz, 16MB L3 Cache, DDR3 1600) and interconnected with a QDR InfiniBand.Wall-clock times for the integration of 1 model year are less than 4 days and the associated cost (write-off and power consumption) are below EUR 100 (based on a machine lifetime of 3 years).If processing speed and cost remain to be linked to Moore's law, we expect a similar cost for a fully coupled biogeochemicalocean-circulation model in 2 years' time from now.

Summary and conclusions
We set out to develop a competitive Baltic Sea oceancirculation-ice model configuration, dubbed MOMBA, which is intended to serve as a nucleus of a coupled oceancirculation-biogeochemical model.To this end, our main focus was a realistic simulation of surface mixed layer dynamics and diapycnal near-surface transport processes, because these processes are key controls on pelagic biogeochemistry (e.g.Sverdrup, 1953 andEppley andPeterson, 1979, respectively).
MOMBA features a horizontal resolution of ≈ 1 nautical mile and a total of 47 vertical layers.It is based on a recent release of GFDL's modular ocean model.The computational cost associated with integrating 1 year is substantial, but yet decades or even centuries are computationally feasible: expressed in terms of wall-clock time the cost ranges, depending on the hardware ("low cost" or high-performance computing), from 4 days to several hours.
The high spacial resolution of MOMBA permits the resolution of eddies.Overall, we find an impressively eddying surface circulation.The associated eddy kinetic energy is, locally, even comparable to open-ocean conditions, which showcases the highly non-linear characteristics of the circulation field in the Baltic Sea.As a peculiarity, we find strong evidence that this non-linearity is further amplified by eddy/wind effects which modulate the wind-induced up-and downwelling.Rough scaling suggests that these eddy/wind effects may well be among the major processes transporting deep nutrient-rich waters to the surface.Despite the overall transient characteristics of the flow field, we simulate local, highly persistent circulation features which are correlated with topographic slopes.Whether this topographic steering of surface currents is indeed as dominant as indicated by the model is, however, yet to be investigated.
A real asset of MOMBA is its realistic simulation of sea surface temperature, which we interpret as indicative of a realistic simulation of surface mixed layer dynamics.Additional evidence indicating that the simulation of near-surface diabatic processes, in general, is realistic is deduced from the model's realistic halocline depths.Overall, this suggests that the simulated physics of MOMBA is of a sufficient fidelity to provide a sound base for research focused on environmental near-surface pollution (such as in, for example, Lehmann et al., 2014a), or on biogeochemical questions such as, for example, the controls of cyanobacteria (which are discussed controversially; Pahlow et al., 2013;Ward et al., 2013).
Caveats are (1) the simulated sea ice extent, which is reasonable with regard to temporal and spacial variability, but is apparently biased high relative to observations, and (2) the failure of MOMBA to reproduce major inflow events.The latter problem is associated with spurious vertical circulation cells developing at the sills during strong inflows.These cells effectively dilute the saline deep inflows with lighter surface waters and thereby inhibit deepwater renewal in the deep Baltic basins.A remedy for this model artefact, which is probably a consequence of the high spacial resolution, is still to be developed.We suggest to investigate approaches similar to the one developed by Kösters et al. (2005) for the Denmark Strait overflow.

Figure 1 .
Figure 1.Horizontal model grid and bathymetry.Panel (a) shows the entire model domain, with red crosses indicating internationally renowned monitoring stations.Panel (b), (c) and (d) are close-ups of the Danish Straits, the Gulf of Finland, and the region between Åland and Finland, respectively.The coloured bins denote the water depth of each model water column.The white areas are considered as land by the model.The actual coastline is denoted by a black line.

Fig. 2 .
Fig. 2. Vertical model grid.The nominal box thicknesses (∆z) as a function of depth.

Fig. 4 .
Fig. 4. Simulated mean (1987 to 1999) surface circulation.(a) The colour and arrows denot the magnitude and direction of the velocity, respectively.Note that the actual spacial resolutio is 20 times higher than indicated by the arrows, both in zonal and meridional direction.(b shows the persistency R of the surface circulation (in colour, as defined by Eq. 15) along wit bathymetry (black contours every 25 m).Note that the actual model domain extends further int the North Sea (c.f.Fig. 1 a).

Figure 4 .
Figure 4. Simulated mean (1987 to 1999) surface circulation.(a)The colour and arrows denote the magnitude and direction of the velocity, respectively.Note that the actual spacial resolution is 20 times higher than indicated by the arrows, both in zonal and meridional direction.(b) shows the persistency R of the surface circulation (in colour, as defined by Eq. 15) along with bathymetry (black contours every 25 m).Note that the actual model domain extends further into the North Sea (cf.Fig.1a).

Fig. 5 .
Fig. 5. Sea surface salinity (PSU) averaged over 1987 to 1999.Panel (a) shows a composite based on ICES bottle observations (upper 10 m, binned monthly into boxes of 6 nm resolution and subsequently averaged over time).Panel (b) refers to MOMBA.Note that the actual model domain extends further into the North Sea (c.f.Fig. 1 a).

Figure 5 .
Figure 5. Sea surface salinity (PSU) averaged over 1987 to 1999.Panel (a) shows a composite based on ICES bottle observations (upper 10 m, binned monthly into boxes of 6 nmi resolution and subsequently averaged over time).Panel (b) refers to MOMBA.Note that the actual model domain extends further into the North Sea (cf.Fig. 1a).

Fig. 7 .Figure 7 .
Fig. 7. Daily mean Ekman pumping calculated from the curl of the wind stress in units day −1 (in color) overlaid by surface currents (arrows) for an exemplary time slice.Strong upwelling is denoted in red.The calculation of the wind stress takes the ocean currents into account (as explained, e.g., inEden and Dietze, 2009).Note, that the Ekman pumping features variability correlated with the surface currents on spacial scales much smaller than those that could possibly be effected by the ≈ 25 km resolution wind fields.

Fig. 8 .Figure 8 .
Fig. 8. Sea surface temperature at station AE, BY15 tively.The positions of the respective stations are i the Swedish Oceanographic Data Centre at SMHI (S

Figure 9 .
Fig. 9. Assessment of simulated sea surface temperature with weekly composites observed from space throughout 1990 to 1999 (BSH SST, Sect.3.3).The color shows the variance explained by the model in units %.

Figure 10 .
Figure 10.Observed (left panel, weekly composite observed from space, BSH SST; Sect.3.3) and simulated (right panel, daily mean) sea surface temperature in degrees Celcius.Note that areas of missing observations are flagged white in both the left and right panel.

Fig. 11 .Figure 11 .
Fig. 11.Near-bottom salinities at station BY2, BY5 and BY15 in panels (a), (b) and (c), respectively.The positions of the respective stations are indicated in Fig. 1.Observations are from the Swedish Oceanographic Data Centre at SMHI (SHARK).

Figure 13 .
Figure 13.Topographic slope and thickness of lowermost wet grid box in MOMBA.Panel (a) shows the difference in the number of grid boxes representing neighbouring water columns.Of the respective differences (up to four; calculated from all neighbours), panel (a) shows the maximum difference.Differences smaller than two grid boxes indicate "resolved" topographic slopes and are masked by white areas.Panel (b) shows the thickness of the lowermost wet grid box (note that we use partial cells).

Figure 14 .
Fig. 14.Simulated deep saline inflow into the Bornholm Basin through Hamrare Strait along t transect indicated in the inlayed map.The color denotes salinity.The thick black dashed li is the surface mixed layer depth as calculated by the KPP scheme.The white (grey) contou denote upwards (downwards) velocities in mday −1 .All properties are averaged over one hou

Fig. 15 .Figure 15 .
Fig. 15.Maximum near-bottom salinity at station BY5 in the period 9 January to 18 April 1993 Panel (a) and (b) show maximum simulated salinities as function of horizontal and vertic Laplacian viscosity, respectively.The dashed black line in panel (a)shows the maximum salini simulated with the "standard" horizontal Laplacian viscosity of 0. Observed salinities peak a 18.5 PSU (Fig.11b).

Fig. 16 .Figure 16 .
Fig. 16.Median salinity profiles for the period 1987-1999 with 10th-and 90th-percentiles at stations AE, BY5 and BY15 in panels (a), (b) and (c) respectively.The locations of the stations are indicated in Fig. 1.Observations refer to SHARK.

Fig. 17 .Figure 17 .
Fig. 17.Ice covered area in the Baltic Sea.The red line denotes an estimate based on digitise ice-charts from the operational ice service provided by SMHI.The black line shows daily mea simulated ice extent.

Figure 18 .
Figure 18.Snapshots of observed (left column) and simulated (right column) ice concentrations for exemplary days in mid-February, during a particularly warm (1989, upper panels) and cold winter (1996, lower panels).Units are in percent.

Fig. 19 .Figure 19 .
Fig. 19.Leading principle components of observed (red line) and simulated (black line) seasonal mean ice concentrations.The blue line (and right y axis) depicts the (unitless) NAO-index for comparison.

Fig. 21 .
Fig. 21.(a) Simulated surface salinity (color) and sea surface height (black contours in unit m).(b) Supposedly identical to (a) but slightly different due to "computational uncertainty".(c temporal evolution of sea surface salinity at station BY2.The black and red line correspon to identical setups which differ due to "computational uncertainty".(d) Temporal evolution o salinity at 80 m depth at station BY5.Black and red lines are coded like in (c).The location o station BY2 and BY5 are shown in the upper panels and in Fig. 1a, respectively.

Fig. 22 .Figure 22 .
Fig.22.Computational performance (wall clock hours required to integrate one year) as a function of allocated computational cores.Black crosses, grey circles and grey x refer to architectures based on Sandybridge, Gainstown and Harpertown CPUs, respectively.The red line depicts a theoretic, ideal parallel performance of the code on the Sandybridge architecture.

Table 1 .
List of model grid boxes where the bathymetry was deepened to increase the salt water supply to the Baltic Sea.

Table 2 .
Volume transport [km 3 ] through the Danish Straits from 5 to 26 January 1993 to the Baltic.