Circulation of the Turkish Straits System under interannual atmospheric forcing

A simulation of the Turkish Straits System (TSS) using a high-resolution, three-dimensional, unstructured mesh ocean circulation model with realistic atmospheric forcing for the 2008-2013 period is presented. The depth of the pycnocline between the upper and lower layers remains stationary after six years of integration, indicating that despite the limitations of the modelling system, the simulation maintains its realism. The solutions capture important responses to high frequency atmospheric events such as the reversal of the upper layer flow in the Bosphorus due to southerly severe storms, i.e., blocking 5 events, to the extent that such storms are present in the forcing dataset. The annual average circulations show two distinct patterns in the Marmara Sea. When the wind stress maximum is localised in the central basin, the Bosphorus jet flows to the south and turns west after reaching the Bozburun Peninsula. In contrast, when the wind stress maximum increases and expands in the north-south direction, the jet deviates to the west before reaching the southern coast and forms a cyclonic gyre in the central basin. In certain years the mean kinetic energy in the northern Marmara Sea is found to be comparable to that of the 10 Bosphorus inflow.

ern coast and three east-west deep basins separated by sills, connected to the two shallow, elongated narrow straits providing passage to the adjacent seas at the two ends, as shown in Fig. 1.
The TSS mass and property balances are mainly controlled by the Black Sea in the upstream. At the Bosphorus Black Sea entrance, the long-term salinity budget implies a ratio of about two between the upper and lower layer volume fluxes (Peneva et al., 2001;Kara et al., 2008). The net flux is estimated to be comparable to the Black Sea river runoff, as the annual average 5 precipitation and evaporation over the sea surface are roughly of the same order (Özsoy and Ünlüata, 1997). Daily to seasonal variations in net fluxes through the TSS are driven by changes in Black Sea river runoffs, barometric pressure and wind forcing.
The hydrodynamic processes of the TSS extend over a wide range of interacting space and time scales. The complex topography of the straits and property distributions have resulted in hydraulic controls being anticipated in both straits (Özsoy 15 et al., 1998; Özsoy et al., 2001), which can only partially be demonstrated by measurements at the northern sill of the Bosphorus (Gregg and Özsoy, 2002;Dorrell et al., 2016). Hydraulic controls have since been found by modelling at the southern contraction-sill complex and the northern sill, confirming a unique maximal exchange regime adjusted to the particular topography and stratification (Sözer and Özsoy, 2017a;Sannino et al., 2017). These findings support the notion that the Bosphorus is the more restrictive of the two straits in controlling the outflow from the Black Sea to the Mediterranean. The analysis of 20 moored measurements by Book et al. (2014) demonstrated this, and indicated a more restrained sea level response transmitted across the Bosphorus than in the Dardanelles. Improvements in modelling have provided a better scientific understanding of the TSS circulation, and they can now address the complex processes characterising the system. The initial step in this formidable task is to construct separate models of the individual compartments of the system, which are the two Straits and the Marmara basin. The first simplified models 25 of the Bosphorus were by Johns and Oguz (1989) who solved the turbulent transport equations in 2D, and found a twolayer stratification to develop. Simplified two-layer or laterally averaged models of the Dardanelles and Bosphorus were later developed by Oguz and Sur (1989), Stashchuk and Hutter (2001) and Oguz et al. (1990) respectively, while Hüsrevoglu (1998) introduced a 2D reduced gravity ocean model of the Dardanelles inflow into the Marmara Sea. Similar 2D laterally averaged models (Maderich and Konstantinov, 2002;Ilıcak et al., 2009;Maderich et al., 2015) and 3D models (Kanarska and Maderich,30 2008; Öztürk et al., 2012), which were of limited extent, have been used to construct simplified solutions for the Bosphorus exchange flows. Bosphorus hydrodynamics were extensively investigated by Sözer and Özsoy (2017a) using a 3D model with turbulence parameterisation under idealised and realistic topography with stratified boundary conditions in adjacent basins, demonstrating the unique hydraulic controls in the maximal exchange regime that are to be established in the realistic case.
The combined effects of the Bosphorus and the proposed parallel channel known as Kanalİstanbul have been investigated by 35 Sözer and Özsoy (2017b), indicating weak coupling between the two channels, which have very different characteristics, but this has been found to be of climatic significance in modifying the fluxes across the TSS.
Very few studies have attempted to model the circulation in the Marmara Sea, even as a stand-alone system excluding the dynamical influences of the straits and atmospheric forcing. Chiggiato et al. (2012)  Similarly, the interannual variability of the Marmara Sea has been examined by Demyshev et al. (2012) using open boundary conditions at the strait junctions in the absence of atmospheric forcing. They reproduced the S-shaped jet current traversing the basin under the isolated conditions of a net barotropic current, which with appropriate parameterisation successfully preserved the sharp interface between the upper and lower layers when the model steady-state was reached after 18 years of simula-10 tion. The S-shaped upper layer circulation of the Marmara Sea predicted by Demyshev et al. (2012) appears similar to what Beşiktepe et al. (1994) found in summer, when wind forcing is at its minimum or at least close to being in a steady-state.
An anti-cyclonic pattern has generally been identified in the central Marmara Sea, like the cases reported by Beşiktepe et al. (1994).
The challenges of modelling the entire TSS domain were recently undertaken by Gürses et al. (2016). The effects of at-15 mospheric forcing were considered, excluding the effects of the net flux through the TSS. The study used an unstructured triangular mesh model, the Finite Element Sea-Ice Ocean Model (FESOM), with a high horizontal resolution reaching about 65 m in the straits in the horizontal. The water column is discretized by 110 vertical levels.
The study of Sannino et al. (2017) used curvilinear coordinate implementation of the MITgcm (Marshall et al., 1997) with a non-uniform grid in the horizontal, a minimum of 65 m resolution in the narrowest part of Bosphorus and 100 levels in the 20 vertical. The model was used to investigate the circulation of the TSS under varying barotropic flow through the system in the absence of atmospheric forcing. The overall circulation in the Marmara Sea was found to differ significantly in each of their experiments with variations in the net volume transport at the Bosphorus. The circulation changes from a large anticyclonic circulation at the centre of the basin at low flux values, to different gyres wrapped around an S-shaped jet as the net flux is increased. A cyclonic central gyre is eventually generated as the increased flux leads to the lower layer in the Bosphorus being 25 blocked. The most significant finding was the nonlinear sea level response, which deviated widely from the linear response predicted by the stand-alone Bosphorus model of Sözer and Özsoy (2017a). Stanev et al. (2017) approached the challenge by using an unstructured mesh model. The model covers the entire Black Sea, and is seamlessly linked to the TSS and the northern Aegean Sea with open boundaries in the south and uses realistic atmospheric forcing. The focus is on the transport at the straits and the resulting dynamics of the Black Sea. Regarding the 30 TSS, the results support the notion of multiple controls of the barotropic flow. They were also able to simulate short-term events such as severe storm passages. However, in a high-resolution model of the TSS, the Courant-Friedrichs-Lewy (CFL) condition (Courant et al., 1967) restricts the time step to a few seconds with explicit schemes. As a solution, Stanev et al. (2017) used an implicit advection scheme for transport to handle a wide range of Courant numbers (Zhang et al., 2016) while satisfying the stability of the solution. However, the computational burden of using an implicit scheme imposed a coarser model resolution with 53 vertical levels at the deepest point of the Black Sea. We note that this limitation may, particularly in the Bosphorus, lead to excessive vertical mixing or a widened interface thickness, which are crucial for the intrusion of the Mediterranean origin water into the Black Sea. This can be seen in their Figures 11c and 11d, for example.
Another recent study by Ferrarin et al. (2018) presents the tidal dynamics in the TSS which is seamlessly included in a model domain covering whole Black Sea -Mediterranean system. They use a barotropic version of the SHYFEM (Umgiesser, 2012) 5 to provide the hydrodynamical background to the tidal analysis.
In our work, we simulate the complete TSS system with a high vertical resolution unstructured grid model forced by complete heat, water and momentum fluxes. A long-term six-year simulation is used to analyse the combined response of the Marmara Sea to atmospheric forcing and strait dynamics. The questions addressed in this paper are: What is the mean Marmara Sea circulation and its variability in a long-term simulation? What are the effects of the atmospheric forcing on the Marmara Sea 10 dynamics and circulation?
The paper is organised as follows: In the next section, we document the model setup and the details of the experiment. In section 3, the validation of the water mass structure and sea level differences along the TSS are demonstrated. The resulting volume transports through the straits, and the kinetic energy and circulation in the Marmara Sea, are presented. Finally, in section 4, we summarise and discuss the results.

5
Models for solving the dynamical equations for an unstructured mesh using finite element or finite volume methods have been implemented for many idealistic applications (White et al., 2008;Ford et al., 2004), and in realistic coastal ocean studies (Zhang et al., 2016;Federico et al., 2017;Stanev et al., 2017). An obvious advantage of using an unstructured mesh model is the varying resolution, which allows for a finer mesh resolution in coastal areas than in the open ocean. The general ocean circulation model used in this study is the Finite Element Sea-ice Ocean Model (FESOM). FESOM is an unstructured mesh 10 ocean model using finite element methods to solve hydrostatic primitive equations with the Boussinesq approximation (Danilov et al., 2004;Wang et al., 2008). We use the initial implementation of Gürses et al. (2016).
The model domain extends zonally from 22.5 • E to 33 • E and meridionally from 38.7 • N to 43 • N, covering a total surface area of 1.52 x 10 11 m 2 (Fig. 1). The mesh resolution is as fine as 65 m in the Bosphorus and 150 m in the Dardanelles. In the Marmara Sea, the resolution is always finer than 1.6 km and is not coarser than 5 km in the Black Sea and the Aegean Sea. The Model equations and parameters are documented in appendix A. The current model implementation considers closed lateral boundaries. Therefore, volume and salinity conservations are imposed to prevent drifts in the tracer fields. Our approach for volume and salinity conservations is described in appendix B.

20
The experiment detailed below was conducted over six years, commencing on 1 January 2008 and continuing until 31 December 2013. The initial fields were obtained after a three-month integration of a lock-exchange case, which was initialised from different temperature and salinity profiles in the basins of the Black, Marmara and Aegean Seas (Gürses et al., 2016;Sannino et al., 2017). Fine mesh resolution and energetic flow structures in the straits require small time steps so the CFL condition is not violated. The time step is thus set to 12 s throughout the integration. The simulation was forced by atmospheric 25 fields provided by ECMWF with 1/8 • resolution. The forcing data cover the whole experiment period with a time frequency of six hours. Precipitation data are obtained from monthly CPC Merged Analysis of Precipitation (CMAP, Xie and Arkin (1997)) and interpolated to the ECMWF grid as daily climatology.
The annual mean of wind, wind stress and wind stress curl for the simulation period are shown in Fig. 2. Wind stress, τ , is calculated as where C d is the drag coefficient and u wind is the wind velocity. The annual mean wind fields are northeasterlies which were strongest in 2011. τ was higher than 0.04 N m −2 in the central-north Marmara Sea in 2009. It then expands in a north-south  A surface area of 22693 km 2 north of 42.5 • N in the Black Sea functions as a buffer zone (the grey shaded area in Fig. 1). This zone is utilised to provide required water fluxes for a realistic barotropic flow through the Bosphorus. The model is forced by a climatological runoff in the Black Sea, which is essential to generate realistic sea level differences between the compartments of the system (Peneva et al., 2001). The monthly runoff climatology for water fluxes was obtained from Kara et al. (2008) and are the same for all the six years. The surface salinity at the buffer zone was relaxed to a monthly climatology, computed from 5 a 15-year simulation by the Copernicus Marine Environment Monitoring Service Black Sea circulation model (Storto et al., 2016). The salinity relaxation time is approximately two days. Although this is a strong constraint, it is required to prevent the surface salinity from decreasing in the buffer zone due to the excessive amount of fresh water input. The climatological values used for runoff and salinity relaxation are shown in Table 1.
Here, we note again that the model having closed lateral boundaries requires a correction for the water and mass fluxes 10 which we discuss in the Appendix A and B. Very briefly, the access or deficit of the water fluxes is redistributed to the surface model nodes weighted by the corresponding element area. Since a runoff climatology is inserted in the Black Sea to maintain the sea level difference along the system, the correction will enforce the volume conservation. This method was also used by Gent et al. (1998) for a climate global model. Another method was explored by Gürses et al. (2016) where the correction was applied only in the Aegean area but salinities became too high and the salinity profile of the simulation became unrealistic. 15

Results
In this section, we present the results obtained from a simulation of the TSS between 2008 and 2013. The focus is on the Marmara Sea and the straits, but we also consider the adjacent basins when necessary. We provide details of the main water mass characteristics of the system and our validation against the observations. The sea level differences and the volume transports through the Bosphorus and Dardanelles are analysed. Although we will include results demonstrating the response of the 20 system to daily atmospheric events, we focus on the interannual changes in the TSS. Therefore, we consider only annual means in time averages. The computed time averages for the simulation are for the period 2009-2013, as the first months of 2008 are considered as an initial spin-up period.

Surface Heat and Water Fluxes
The monthly averages of water fluxes and net heat flux are shown in Fig. 3a for the whole domain and Fig. 3b for the Marmara Sea. The mean runoff, evaporation and precipitation in the whole domain are 278.1 km 3 /yr, -131.5 km 3 /yr and 70.8 km 3 /yr, respectively. The resulting net water flux into the ocean is therefore 226.4 km 3 /yr which is balanced by the correction term (see Appendix A) shown in red in Fig. 3a. As a result the model volume is conserved.  In the Marmara Sea, runoff is identically zero Fig. 3b. Mean evaporation and precipitation are -8 km 3 /yr and 6.8 km 3 /yr, respectively. The water flux correction applied in the Marmara Sea is -15.9 km 3 /yr. This means there is a net surface flux loss in the Marmara Sea about -17.1 km 3 /yr. The total water budget of the Marmara Sea will be further discussed in section 3.3.
The net heat flux in the Marmara Sea is calculated in the range of -123.3 W/m 2 to 138.4 W/m 2 with minimums and maximums in December-January and May-June, respectively.
where α and β are thermal and haline expansion coefficients (McDougall, 1987), Q H is the heat flux (positive towards the ocean), C w is the specific heat capacity, S 0 is the surface salinity and Q w is the water flux. 10 The Black Sea and the Marmara Sea gain buoyancy except for a small area near their western coasts (Fig. 4), while the Aegean Sea has a buoyancy loss except near the Dardanelles exit and the Anatolian coast. The average buoyancy flux changes between -7x10 −8 and 3.4x10 −8 m 2 s −3 over the domain. It does not show significant spatial differences interannually, but the gradient between the Aegean and Black Seas was stronger in 2011 than in the other years (not shown).

15
The surface salinity ranges from 16 to 38 psu over the whole domain (Fig. 5). The surface waters leave the Bosphorus and the Dardanelles with salinities of about 21 psu and 27 psu, respectively. The surface salinity in the northern Marmara Sea is less than 23 psu and increases to 25 psu in the south. Long-term measurements from 1986 to 1992 in the Marmara Sea (Beşiktepe et al., 1994) suggest a salinity ranging between 23±2 psu at the surface, which is satisfied by the simulation. The vertical structure of the time mean salinity and temperature along the thalweg (see Fig.1) is shown for the last year of the integration in Fig. 6. The depth of the interface between the upper and lower layers is stationary throughout the simulation and is located at a depth of around 20 m. The water column salinity is mixed below 25 m in the Marmara Sea until approximately 200 m. However, the deeper layers are colder and continuously modified by the Dardanelles inflow after the initialization as can be seen in Fig.7 in which the volume temperature has a negative trend and salinity a positive one. We argue that this is due 5 to initialization adjustment and the closed boundary conditions. Salinity is completely mixed in the Marmara Sea below 100 m depth.
In    the Dardanelles side of the Marmara Sea than the Bosphorus side. The RMS error is larger than those found in operational models of the Mediterranean Sea (Oddo et al., 2009), but this is due to a thermocline and halocline vertical shift, as shown in Fig. 9a. The vertical mixing and the missing interannual variability of the Black Sea runoff probably account for this. The model performs substantially better at the surface and below 30 m in depth than it does at the depth of the interface between the upper and lower layers of around 20 m (Fig. 9b). Temperature RMS errors are highest around the seasonal thermocline in June 2013, while at the same level as the halocline in the other two months.

Sea level and volume fluxes through the straits
The time mean of the sea surface height across the system is shown in Fig. 10   Moller (1928) measured the sea level differences between the two ends of the Dardanelles and the Bosphorus as 7 and 6 cm, respectively. Gunnerson and Ozturgut (1974) and Büyükay (1989) found the sea level difference in the Bosphorus to be 35 cm for the 1966-1967 period and 28-29 cm for 1985-1986 period, respectively. Bogdanova (1969) gave a sea level difference of 42 cm between the northern Black Sea (Ialta) and southern coast of Turkey (Antalya) with high seasonal variability. Alpar contraction, suggesting that the hydraulic control there has been lost. However, the thin interface layer on top of the northern sill and the flow north of it continues to preserve its shape, suggesting continued hydraulic control at the northern sill.
All the above features are consistent with the findings obtained from idealised and realistic models of the Bosphorus provided in Sözer and Özsoy (2017a) and the TSS model experiments of Sannino et al. (2017).
In the Dardanelles, the two-layered water mass and flow structure can be seen in Fig.13a. The simulated layered-structure in in Nov. 22, 2008 results in extensive mixing in the upper layer (Fig.13b) where the salinity increases substantially. Furthermore, even the stratification is broken partially in the Marmara Sea side of the strait.
Upper layer velocity in the southern exits of both the Bosphorus and Dardanelles are generally higher than their northern exits, due to the hydraulic controls exerted at the constrictions in the middle of both straits and the expansion area in the mouth (Sözer and Özsoy, 2017a). Conversely, the lower layer velocity has much higher maxima at the northern exits in both straits, 5 and are roughly 0.2 m/s less than the measurements by Jarosz et al. (2011aJarosz et al. ( , 2012. The upper layer maximum velocities are in accordance with most observations. The depth of the zero velocity level varies significantly in the exits of each strait and are listed in Table 3. These depths are consistent with the recent approximations reported in Jarosz et al. (2011aJarosz et al. ( , 2012 Table 4. Annual mean of net, upper layer and lower layer volume fluxes (km 3 /yr) for the whole simulation period. A negative value means the flux is in the Black Sea-Aegean Sea direction.
The net, annual mean upper layer and lower layer volume fluxes are given in Table 4 (Figure 14). In the northern Dardanelles, there is a large discrepancy between the simulated net fluxes and the estimates from observations. Sannino et al. (2017) and Özsoy and Altıok (2016) both have similar disagreements with the measurements of Jarosz et al. (2013) and have concluded that the discrepancies could be a result of measurement or computational inaccuracies at the wide 10 northern section of the Dardanelles, where the instrument data were located.
The historical estimates of the Dardanelles layer volume fluxes show much higher values than our simulation ( Table 5). The changes found in the layer flux from one end of the strait to the other have been attributed to turbulent entrainment processes, which transport water and properties across the hypothesised layer interface Özsoy et al., 2001;Özsoy and Altıok, 2016). The larger disagreement of baroclinic volume fluxes compared to the net fluxes, and the lower estimations of    Table 5. Annual means of volume fluxes (km 3 /yr) through the Dardanelles estimated by different studies. UL, LL and Net stand for upper layer, lower layer and net fluxes, respectively. A negative value means the volume flux is from the Marmara Sea to the Aegean Sea. boundary model by correcting the water flux every time step comes with the limitation of reducing the transport along its way from the Black Sea to the Aegean Sea. The total water flux over the Marmara Sea is −17.1km 3 /yr (Fig.3b) which is balanced by the difference between the southern Bosphorus and northern Dardanelles (Table 4). Due to the surface flux correction the net outflow in the Dardanelles is decreased by about 10% with respect what entered the Bosphorus. However, such correction allowed us to obtain realistic salinity profiles with respect to other correction methods such as Gürses et al. (2016).

5
Finally, the differences in transport between the transects at the two ends of each strait are 0.5km 3 /yr and 2.1km 3 /yr in the Bosphorus and Dardanelles, respectively. We attribute these differences partially again to the water flux correction at the surface. The rest is due to the numerical errors accumulating during the post-processing of the model outputs.

Marmara sea dynamics and circulation
Using the six-year simulation, it is now possible to estimate the kinetic energy input by the wind in the Marmara Sea. 10 The time series of monthly mean wind stress is shown in Fig. 15. It exhibits interannual differences with a mean of about 0.03 N/m −2 and a maximum of around 0.05 N m −2 in August 2008. The monthly mean is highly variable and there is not a well-defined seasonal cycle between 2008-2013. The wind work normalised by the density and surface area is computed as: where ρ 0 is the surface density, τ is the wind stress, u s is the current velocity at the surface and A is the surface area of the 15 Marmara Sea. The wind work is positively correlated with the wind stress. It is highest in 2011 and the maximum monthly mean wind work is 9.1x10 −6 m 3 s −3 in the autumn of 2011.
To compare with the other marginal seas described in Cessi et al. (2014), we normalise the integral in (3) by the volume of the Marmara Sea. The six-year mean of the wind work is computed as 1.09x10 −8 m 2 s −3 , one order of magnitude higher than the Mediterranean Sea. The wind work in the Baltic Sea was computed to be 9.15x10 −9 m 2 s −3 , which is comparable to the 20 Marmara Sea but is still lower.
The volume-mean kinetic energy in the Marmara Sea normalised by the unit mass is calculated as 0.006 m 2 s −2 for six years. The daily mean kinetic energy time series reveals that severe atmospheric events are able to energise the basin up to In the Marmara Sea, the buoyancy gain is mainly due to the Bosphorus inflow, and thus the latter competes with the wind work to change the kinetic energy of the basin, as explained by Cessi et al. (2014). However, at the surface the kinetic energy shows that the Bosphorus surface jet energises the northeastern basin in addition to the wind (Fig. 16). The kinetic energy of the Bosphorus inflow is always greater than 0.075 m 2 s −2 . In 2009, a kinetic energy maximum appears in the north of the Bozburun peninsula where the Bosphorus jet arrives. In the western basin, the kinetic energy is generally less than 0.025  Sea after turning south, merging with the flow circulating around the islands in the southwestern Marmara Sea and eventually exiting from the Dardanelles. This circulation pattern in the western Marmara Sea is persistent throughout the simulation but with different intensities. This type of circulation structure has been reported in other studies (Chiggiato et al., 2012;Beşiktepe et al., 1994). In 2011 (Fig. 18c), and the circulation in the middle of the Marmara Sea evolves into a single cyclonic structure.
The shift in the circulation can be explained by the shift of the wind stress maximum towards the north (Fig. 2). Sannino et al. 5 (2017) demonstrates a similar cyclonic pattern in the central Marmara Sea due to the potential vorticity input by the Bosphorus.
However, in our case, the main driver of the cyclone should be the wind as the volume transport through the Bosphorus is much lower than that of Sannino et al. (2017)

Summary and discussion
We present a six-year simulation of the Turkish Straits System, with a highly resolved unstructured mesh model of the TSS using realistic atmospheric forcing to identify the circulation structure and the interplay between the atmospheric forcing and the Bosphorus input. The results demonstrate that a realistic representation of the pycnocline from the two Straits to the This type of circulation has been observed and modelled by several earlier studies (Beşiktepe et al., 1994;Chiggiato et al., 2012;Sannino et al., 2017). The other circulation structure includes a cyclone in the central basin because of intensifying and The long-term simulation with atmospheric forcing made it possible to evaluate the wind energy input and compute the 15 kinetic energy in the Marmara Sea. The wind work in the Marmara Sea is shown to be even higher than the Baltic Sea. The high energy input from the wind significantly increases the kinetic energy in the Marmara Sea. During severe storms, kinetic energy can increase by 10 times the time-averaged value. The annual mean of kinetic energy in the regions under the influence of the wind forcing can also exceed that from the Bosphorus jet, depending on the wind stress structure.
In our modelling approach, we focused our attention on the TSS proper, which is the central domain including the Straits and the Marmara Sea, while assigning a limited storage role to the truncated exterior domains of the Black and Aegean Seas.
This represents a trade-off between a truthful reproduction of the TSS circulation and the ability to impose far field boundary 5 conditions in artificial closed basins at the two ends, using a similar strategy to that of Sannino et al. (2017). In future studies, we expect to obtain improved results by incorporating lateral open ocean boundary conditions in the Aegean and Black Seas.
The skill of the model predictions also appeared sensitive to the accuracy of the atmospheric forcing used in the simulation.
Within the relatively small domain of the TSS, an improved representation of the atmospheric forcing, particularly during the severe storms frequenting the region in winter, appears to be essential for improving its skill. The momentum equations are: where u = (u, v) and v = (u, v, w) are 2D and 3D velocities, respectively in the spherical coordinate system, ρ 0 is the mean density, g is the gravitational acceleration, η is the sea surface elevation, f is the Coriolis parameter andk is the vertical unit 20 vector. ∇ and ∇ 3 stand for 2D and 3D gradients or divergence operators, respectively. The horizontal and vertical viscosities are denoted by A h and A v . p is the hydrostatic pressure obtained integrating the hydrostatic relation (A3) from z = η. Atmospheric pressure is not considered not to excite basin-modes in a closed lateral boundary model.
The Laplacian viscosity is known to be generally too damping and strongly reduces the eddy variances of all fields compared to observations when the model is run at eddy resolving resolutions (Wang et al., 2008). Therefore, biharmonic viscosity is 25 used in the momentum equations. Here, A h is scaled by the cube of the element size with a reference value A h0 of 2.7 × 10 13 m 4 /s (Table A1), which is set for the reference resolution of 1 degree.
The continuity equation is used to diagnose the vertical velocity w: and the hydrostatic equation is: where ρ is the deviation from the mean density ρ 0 .
Tracer equations (A4) and (A5) are solved for the potential temperature, T , and salinity, S where K h and K v are the horizontal and vertical eddy diffusivities, respectively. Laplacian diffusivity is used for the tracer equations. K h is again scaled as A h but by the element size with a reference value of 2.0 × 10 3 m 2 /s. These values are set following the convergence study of Wallcraft et al. (2005).
The Pacanowski and Philander (1981) (PP) parametrization scheme is used for vertical mixing, with a background vertical viscosity of 10 −5 m 2 /s for momentum and diffusivity of 10 −6 m 2 /s for tracers. The maximum value is set to 0.005 m 2 /s. 10 The density anomaly ρ is computed by the full equation of state (A6). ρ = ρ(T, S, p) The surface and bottom momentum boundary conditions are, respectively: where τ and C d are the wind stress and the surface drag coefficient, respectively.
The surface kinematic boundary condition is: where E (m/s), P (m/s) are evaporation and precipitation, respectively. R (km 3 /yr) is runoff (see Table1). It is converted to 20 m 3 /s and normalised by the area of the buffer zone in the Black Sea (see Fig. 1) before entering to the equation. Finally, W corr is a correction applied to conserve the volume of the model, as described later.
The sea surface height equation can now be derived from equations A2 and A9 as: The upper limit of integration in (A10) is set to η in this version of FESOM and is different from Wang et al. (2008) to provide 25 a non-linear free surface solution.
The bottom boundary condition for the temperature and salinity are (∇T, ∂ z T ) · n 3 = 0 (A11) (∇S, ∂ z S) · n 3 = 0 (A12) where n 3 is the 3D unit vector normal to the respective surface.

5
The surface boundary condition for temperature is where C p = 4000J/(kg K) and Q (W/m 2 ) is the surface net heat flux into the ocean.
In global applications, surface salinity is generally relaxed to a climatology to prevent a drift. In our regional application, the water flux term in the boundary condition (A14) is applied over the whole domain whereas the relaxation term is prescribed 10 only in the Black Sea buffer zone.
In (A14), S 0 and S * are the surface salinity and the reference salinity, respectively, γ is the relaxation coefficient (Table A1).
Finally, S corr is the counterpart of W corr for salinity conservation corresponding to boundary conditions (A14), which will be defined in the following section.

Appendix B: Salt Conservation Properties
As our model domain is closed we need to enforce salt conservation. Volume salinity conservation requires the time rate of the change in the volume salinity term in equation (B1) to be zero. A balance must be satisfied between the two integrals.

∂ ∂t
In FESOM, this balance is achieved by applying a correction for each term separately. The amount of water flux by evap-20 oration, precipitation and runoff is integrated over the surface with every time step (Equation B2). After normalising by the domain area as in (B3), the surplus or deficit is added to or subtracted from the total water flux equally from each node of the mesh with the W corr terms in equations (A9) and (A10).
The salinity flux is corrected in a similar manner for the boundary condition

∆S =
x,y (S 0 (E − P − R) + γ(S * − S 0 ))dxdy (B4) After applying these corrections, we get a surplus of water corresponding to about 1 mm of sea surface height increase a year, which we believe is due to random numerical errors. Correspondingly, the volume-mean salinity decreases by an order of 10 −5 psu a year. Although these errors may be significant in climate scales, they are acceptable for our six-year long experiment.