A model of Black Sea circulation with strait exchange (2008–2018)

The Bosphorus exchange is of critical importance for hydrodynamics and hydroclimatology of the Black Sea. In this study, we report on the development of a mediumresolution circulation model of the Black Sea, making use of surface atmospheric forcing with high space and time resolution, climatic river fluxes and strait exchange, enabled by adding elementary details of strait and coastal topography and seasonal hydrology specified in an artificial box on the Marmara Sea side. Particular attention is given to circulation, mixing and convective water mass formation processes in the model, which are then compared with observations. Open boundary conditions relaxed to seasonal hydrology specified in the artificial box are found to enable Bosphorus exchange with a proper upper layer, lower layer and net fluxes comparable to the observed ranges. These improvements at the artificial boundary and in the interior evolution of the Black Sea allow the study to capture daily, seasonal to decadal climatic variability and change observed in the Black Sea in the last few decades.


A short review of Black Sea oceanography and recent modeling development
The Black Sea is a part of the old world seas that has limited communication with the World Ocean. Its unique oceanographic characteristics, natural and environmental history have been subject to a number of reviews by Özsoy and Ünlüata (1997, 1998); Oguz et al. (2005); Kosarev (2007); Sorokin (2002); Grinevetsky et al. (2002); Vespremeanu and Golumbeanu (2017). On the other hand, there is clear observational evidence of rapid environmental change in the Black Sea in the recent decades (Capet et al., 2016;Stanev et al., 2019). Fluxes of momentum, water and buoyancy at the sea surface and at coastal and open boundaries dominate the hydrodynamics of semi-enclosed seas. Lateral fluxes, especially those articulated in the Black Sea, by inputs from major rivers and the Bosphorus Strait exchange flows demand their adequate representation in this almost totally enclosed sea. Despite relatively small net flux on the order of 0.01 Sv across the strait, strong coupling between the Bosphorus and the Black Sea is suggested by observations, in terms of circulation and mixing in either region. Various modeling studies of the Black Sea circulation, often with objectives similar to ours, have considered realistic model configuration, initialization and time-dependent forcing elements aiming to generate results verifiable by observations.
To the best of authors' knowledge, modeling of the Black Sea circulation and hydrography allowing for a system of open boundary conditions applied at the Bosphorus has not been attempted in the available literature so far, although a number of stand-alone Black Sea models have attempted indirectly accounting for fluxes at the Bosphorus (Stanev et al., 2003Stanev, 2005;Miladinova et al., 2017Miladinova et al., , 2018. Previous modeling efforts in the Black Sea considered the Bosphorus exchange fluxes by using different open boundary condition techniques. For example, Staneva et al. (2001) expressed the Bosphorus fluxes by spreading the exchange over five horizontal and five vertical grid points with speci-fied upper and lower salinity in their DieCAST ocean model. (Kara et al., 2008) expressed the imbalance of water budget by adding negative river precipitation (i.e., a river evaporation) at the Bosphorus exit on their HYbrid Coordinate Ocean Model (HYCOM). Miladinova et al. (2017) described the Bosphorus inflow/outflow as a river flow by using the 3-D General Estuarine Transport Model (GETM). None of the previous studies included the Bosphorus exchange dynamics under the real-time atmospheric forcing and circulation, while in the present study we aim to achieve coupling with Bosphorus Strait under real-time forcing by extending the domain to include a portion of the Marmara Sea, as will be detailed in the following sections.
Our model strategy, on the other hand, could be compared to the addition of an "Atlantic box" preferred in the early phases of Mediterranean Forecasting System (MFS) (Oddo and Pinardi, 2008), which has only recently been updated to involve further refinements of coupled systems (CMEMS, 2017). Recently, however, there have been various efforts to couple the entire series of straits and basins of particular characteristics together, by making use of high-resolution unstructured meshes (Ferrarin et al., 2018;Federico et al., 2017;CMEMS, 2017) which have yet to survive the various obstacles to properly represent coupling for each of the straits (not only Bosphorus and Black Sea but others), realistically accounting for flux variability between the various coastal and basin-scale elements.
We skip further review of the voluminous literature on stand-alone modeling of the Black Sea, while in the present study we aim to achieve coupling with Bosphorus Strait, extending the domain to include a portion of the Marmara Sea and also the Sea of Azov, as will be detailed in the following sections.

Interaction with straits
A unique regime of hydraulically controlled, turbulent and strongly stratified flow is well documented by observations in the Bosphorus (Gregg et al., 1999;Özsoy and Ünlüata, 1998;Özsoy et al., 2001;Gregg and Özsoy, 2002).
Observations in the Marmara Sea have revealed two-layer stratified flow details, fluxes and circulation in this small internal sea (Ünlüata et al., 1990;Beşiktepe et al., 1993Beşiktepe et al., , 1994Özsoy et al., 2016), connected to the Black Sea and Aegean Sea, respectively, by the Bosphorus and Dardanelles straits.
Estimates of Bosphorus fluxes contributing to the Black Sea have been made based on satellite observations and numerical models (Peneva et al., 2001;Staneva et al., 2001;. Often, these indirect methods have yielded inaccurate estimates of the net and layer-averaged two-way fluxes through the Bosphorus by deforming the statistics of the short-term variability often displayed at coastal scales of the strait. In addition to estimates of annual and seasonal fluxes by in situ measurements obtained throughout the years, recent observations based on instrumental measurements continued for extended periods have been of great value to establish better statistical reliability in these estimates (Özsoy et al., 2001;Jarosz et al., 2011b, a). Collective re-evaluation of these recent observations has led to better estimates that are typical of the seasonal fluxes through the Bosphorus with relatively small trends and observational uncertainties within the last decades (Schroeder et al., 2012;Altıok and Kayışoglu, 2015;Özsoy and Altıok, 2016a, b;Jordà et al., 2017).
Recent advances in modeling (Ilıcak et al., 2009;Sözer and Özsoy, 2017a, b;Sannino et al., 2017) have provided details of the dynamics and mixing as well as the essential verification of the experimental synthesis of hydraulic controls, also providing evidence for "maximal exchange" special conditions satisfied in the Bosphorus. The Bosphorus Strait has been found to be a special environment that supports maximal exchange, even more strictly than the Strait of Gibraltar for which the theory has originally been developed Armi, 1986, 1988). Later developments in modeling have enabled resolving the coupled dynamics of the Turkish Straits System (TSS) including the Bosphorus and Dardanelles straits as well as the Marmara Sea in the model domain (Sannino et al., 2014Aydogdu et al., 2018a, b).
There is evidence that the interaction between the Black Sea and the TSS occurs in both ways and in different realms. For instance, in the Black Sea, which is located in a strong climate gradient between the warm Mediterranean and continental Eurasia, an upper-ocean cold water mass is formed in winter (the cold intermediate layer -CIL), which formerly used to prevail throughout the year, although recent evidence (additionally in the present study) shows it to be disappearing fast. The CIL resides just above the halocline and provides a fundamental contribution to static stability of the water column. At the same time, the CIL resides partially above the wide southern continental shelf and continental slope regions ( Fig. 1) where it is responsible for modifying the characteristics of the bottom Mediterranean plume exiting from the Bosphorus into the Black Sea. The CIL penetrates into the Bosphorus, coming into contact with the Mediterranean waters that it overlies, influenced by turbulent mixing and entrainment between waters of Mediterranean and Black Sea origin, transformed towards the sharply stratified surface and bottom waters that emerge, respectively, on the Marmara and Black Sea sides of the strait (Gregg et al., 1999;Özsoy and Ünlüata, 1998;Gregg and Özsoy, 2002). In this way, the CIL, which is a product of convective mixing in the Black Sea, influences water mixed on the shelf and returned back to deeper layers of the Black Sea, while also influencing the Marmara Sea.
In correspondence with the above example of Black Sea influence on the TSS, the TSS also predetermines what would later become of its waters flowing to the Black Sea.
The Mediterranean water plume exiting onto the Black Sea shallow continental shelf spreads along the shelf and canyon finally dips down the continental slope, where it creates intrusions (Latif et al., 1991;Murray et al., 1991;Özsoy et al., 1993;Özsoy and Beşiktepe, 1995) leading to double diffusive convection Özsoy and Beşiktepe, 1995;Kelley et al., 2003) and in extreme situations can penetrate throughout the Black Sea interior (Falina et al., 2017), interestingly found responsible for fast penetration of chemical signals such as the Chernobyl radiation fallout in the interior Black Sea (Rank et al., 1999;Özsoy et al., 2002;Delfanti et al., 2014). Novel features detected in the region of spreading along canyon topography of the Black Sea shelf adjacent to the Bosphorus have led to detailed description of density currents along undersea channels with new hypotheses proposed on their turbulent dynamics (Dorrell et al., 2016(Dorrell et al., , 2019. Based on the above literature, it is evident that the evolution of water properties and currents through the TSS influence mixing and water mass formation processes in the Black Sea shelf and continental slope boundary layers as well as the interior region. With regard to coupled modeling of Black Sea and the TSS, we can quote Sannino et al. (2017) results showing the changes in circulation of the Marmara Sea in response to net flows through the Bosphorus, finally leading to switching of the modes of circulation in response to extreme cases of flow blocking in the Bosphorus. Similarly, Stanev et al. (2017) showed significant changes in coupled circulations of the Black Sea and Marmara Sea in response to changes in flow through straits.

Model configuration
We aim to closely reproduce the circulation and hydrography of the Black Sea during the last decade by developing a unique high-resolution model prototype that relies on highfidelity forcing data sets available today on atmospheric reanalysis and climatological riverine fluxes, also considering coupling with the external ocean by including a "Marmara box" and the Bosphorus Strait in communication with the Black Sea. In this respect, the model aims at truthful reproduction of the interplay between thermohaline circulation, water mass formation dynamics and interbasin exchange in coupled fashion. Our renewed effort aims to improve the great number of previous modeling studies that have simulated some characteristics of the Black Sea circulation without reference to Bosphorus exchange, thereby to close gaps in predictability of the physical response of this complex basin on scales from days to a decade by introducing the essential strait coupling.
It is with worth noting that the original configuration of Nucleus for European Modelling of the Ocean (NEMO) version 3.3.1 we used earlier had a horizontal resolution of 1/12 • grid and 42 vertical levels. Having a coarser grid, us- ing an earlier version of General Bathymetric Chart of the Oceans (GEBCO) topography, Levitus gridded climatology for initial conditions, atmospheric fluxes based on ECMWF ERA-Interim reanalyses and only major rivers such as the Danube, Dniester, Dnieper and Don included, the strait exchange was not sufficiently representative, and an update seemed in order as we moved to the present model configuration for improved results.

Model domain, grid and bathymetry
The present Black Sea circulation model (BSEA) is based on NEMO version 4.0 (Madec, 2008) with advanced features of a community model accounting for complex physics of the marine environment. The model domain between 40.69-47.30 • N and 27.43-42.00 • E includes the deep Black Sea basin (max depth: 2178 m) together with the shallow Sea of Azov (depth: 10 m) in the north and the Marmara Sea in the south. The bathymetry for the model displayed in Fig. 1 is based on the General Bathymetric Chart of the Oceans (GEBCO08) grid version 20100927 released in 2010, which we find to be the most up-to-date and accurate data set for the Black Sea, with 30 arcsec resolution.
The model domain includes Sea of Azov in the north, the shallowest sea in the world, with a maximum depth of 14 m in the middle, separated from the main body of the Black Sea by the narrow Kerch Strait. The river Don joins Sea of Azov at its head, while the river Kuban joins near the Kerch Strait. The bathymetry and width of Kerch Strait has been artificially increased in the model ( Fig. 1) in order to have unobstructed exchange with the Black Sea, necessitated by the high river discharge into this relatively small, shallow domain with very sensitive response to wind stress forcing.
In order to represent the Bosphorus exchange flows coupled to the Black Sea, an artificial box representing the Marmara Sea and the Bosphorus Strait has been added to the domain. The Marmara box is bounded by the northern coast of the Marmara Sea including the Bosphorus entrance and three open boundaries: (i) along 40.74 • N, (ii) 28.403 • E and (iii) 29.28 • E, with the depth limited to 1200 m at the bottom and true bathymetry elsewhere.
The model nominal grid resolution is x ≡ y ≡ 2.5 km (0.03133 • longitude by 0.02286 • latitude) along both horizontal directions, resulting in an equally spaced horizontal grid of dimension 466×280 and a vertical grid with 60 z levels of partial step configuration. The resolution of the vertical grid has been adjusted so as to be compatible with the topography and density variations. The minimum vertical grid size of 2 m at the surface, increasing to about 3 m at a depth of 40 m, then to about 10 m at a depth of 110 m, is designed to adequately represent the surface mixed layer, which is typically shallower than 25 m, and the main halocline, which is usually shallower than 100 m in the Black Sea. There are 20 levels in the first 50 m of depth, followed by 10 levels in the next depth range until 100 m (Fig. 2). The wide shelf areas, especially of the western Black Sea typically limited by the 100 m depth contour at the shelf edge, are sufficiently resolved by this vertical grid. The model integration was carried out with a baroclinic time step of t i = 240 s.
The model bathymetry represented in Fig. 1 has been moderately smoothed using scale-selective filters, satisfying the slope criterion of the Regional Ocean Modeling System (ROMS) (Song and Haidvogel, 1994) so as to have rx0 = max( h/ h) < 0.3 between any neighboring grid points, where h is the depth.
In regard to the Bosphorus coupling, it should be clear that the present horizontal and vertical grid resolution needed to represent the complicated bathymetry at the strait is lower than what would be needed. It is therefore to be expected that the present configuration insufficiently resolving the strait topography would not be capable of reproducing the full range of nonlinear dynamical behavior of the Bosphorus Strait. The present choice is made in recognition of this fact in order to include strait coupling to the lowest order of approximation, saving the required effort for the future.

Parameterization
The numerical values of the horizontal viscosity length and velocity are taken at 1 km and 0.1 m s −1 , respectively. The horizontal diffusive length scale and velocity are 2 km and 0.01 m s −1 , respectively. Values for vertical eddy viscosity and eddy diffusivity are selected as A v m = 1.e−7 and A v t = 1.e − 8 m 2 s −1 . A turbulent kinetic energy (TKE) turbulent closure scheme (Blanke and Delecluse, 1993) is applied in the vertical with default settings in NEMO version 4.0 con- figuration. Bilaplacian diffusivity schemes are used for both tracer and momentum.

River input
The inputs of freshwater from major rivers, Danube in Romania, Dniester, Southern Bug and Dnieper in Ukraine, Don, Kuban in Russia, Rioni in Georgia, Ki˙zi˙li˙rmak and Sakarya in Turkey, and Kamchiya in Bulgaria have been accounted for by using climatological monthly discharge reported by the Global River Discharge (RivDIS) database. Climatological monthly discharge values of the rivers are given in Table 1.
The total annual river inflow represented with 11 rivers in the model amounts to 10 570 m 3 s −1 (= 333 km 3 yr −1 ). Kara et al. (2008) have found little difference in comparing different sources of river runoff, including a budget based on six major rivers of the Black Sea, totaling 9160 m 3 s −1 (= 289 km 3 yr −1 ).
More importantly, additional land-based sources of water other than rivers may be more important, though often neglected. If the overland runoff were to be taken into account, the total freshwater inputs apparently would increase by another 30 %, amounting to a total of 13 825 m 3 s −1 (= 436 km 3 yr −1 ) according to Black Sea data displayed on the Global Runoff Data Centre (GRDC) web portal (https: //www.bafg.de/GRDC/EN/Home/homepage_node.html, last access: 11 January 2020) "freshwater fluxes to the world oceans", although the monthly data for the overland flows could not be separately accessed and are therefore not included in our account.
According to the same source (GRDC), the total freshwater runoff has significant interannual variability. As it can be observed in Fig. 3, the total discharge of freshwater into the Black Sea by rivers plus the overland runoff varies by almost 2-fold from one year to the other. Considering sea-  sonal variability within each year, this variability is even higher on interannual timescales, as demonstrated by Özsoy and Ünlüata (1998), based on long-term data provided by Simonov and Altman (1991). The Danube river discharge monitored for last two centuries clearly displays this great interannual/interdecadal variability (Bondar, 1986;Sur et al., 1996). The air-sea exchange of water by precipitation and evaporation adds further variability to the water exchange. From the above discussion, it is obvious that the water budget of the Black Sea is highly variable on interannual and seasonal timescales because of climatological factors and is influenced by insufficient accounting of all freshwater sources such as overland runoff, implying that the model response to freshwater inputs is expected to vary greatly at any instance. For the moment, we use the climatological monthly river fluxes as lateral sources of freshwater in the model, although we fully recognize the future need to account for the net water exchange regulated by the Bosphorus Strait.

Open boundary conditions
Within the Marmara box domain external to the Black Sea, temperature and salinity were relaxed to climatological values, which were a function of depth only, along a buffer zone that is five grids wide in parallel to the open boundary, only leaving out a very small area near the Bosphorus entrance quickly adjusted to the rest of the domain shortly after initialization. The depth-and time-dependent values imposed at the open boundary were obtained from historical data sets as daily climatology, by averaging over corresponding dates of years with suitable smoothing and static stability checks among available archives of conductivity-temperature-depth (CTD) data at Institute of Marine Sciences -Middle East Technical University (IMS-METU). The averaging of the data that were finally applied at the Marmara box was done over the entire Marmara Sea for increased statistical reliability of the observational daily T , S profile time series, which are shown in Fig. 4. The data were used to initialize and later to relax for the properties within the Marmara box.

Initialization and atmospheric forcing
The initial conditions of the model are constructed from intercalibrated temperature and salinity data from the 1992 CoMSBlack cruise, a multi-national collaborative oceanographic survey which quasi-synoptically covered the entire Black Sea in maximum possible extent and detail (https://isramar.ocean.org.il/PERSEUS_ Data/CruiseInfo.aspx?criuseid=2318, last access: 11 January 2020), carried out during 2-26 July 1992 by five research vessels from the riparian countries, as quoted in Oguz et al. (1993) and Sur et al. (1994Sur et al. ( , 1996. The data obtained from a total of 394 hydrography stations at nominal spacing of 20 (= 1/3 • ) in longitude and latitude have been screened to eliminate few stations contributing to anomalous properties in the analyzed fields. The objectively analyzed temperature and salinity fields at the surface (model level 1) and at a depth of 89 m (model level 28), respectively, are shown in Fig. 5.
The model was initialized with the July 1992 analysis fields substituted for the 1 July 2008 initial conditions and run for about 10.5 years from 1 July 2008 to 31 December 2018, subject to the atmospheric forcing and surface and lateral flux boundary conditions as specified above. The substitution of the July 1992 analysis for initial conditions is misplaced by more than two decades from the modeling period, but it was justified in the absence of other reliable initial data, providing the best possible synthesis of the entire basin hydrography ever made. Although the model would need a few years of spinup time for adjustment to the initial and boundary conditions, the adjustment run was not carried out, as the objective of model performance testing prevailed. Results obtained so far illustrate rather fast adjustments of the model fields without much vertical displacement or oscillations encountered, although this remains to be finally tested in the future by persistence runs or multidecadal runs from actual initial conditions. We choose to start from a representative initial condition, albeit the fact that it is misplaced in time allows the model dynamics with realistic forcing alone to set up the circulation, without any climatological adjustments of the initial or evolving fields.
Surface fluxes are computed using the core algorithm of NEMO, based on the European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis product ERA5 data set. The ERA5 high space and time resolution atmospheric data at 1 h intervals and 31 km horizontal grid spacing of the atmospheric model are interpolated on the fly by the NEMO ocean model and applied at the surface using the Coupled Ocean-Atmosphere Response Experiment (COARE) 3.5 bulk formulae (Madec and the NEMO Team, 2012).
First model runs have been carried out with the ERA5 precipitation. However, the comparison of model salinity with Argo floats showed a difference in trends between model and observations at the end of model simulation. Detailed analyses have revealed that the ERA5 precipitation is significantly lower than the ERA-Interim precipitation in the Black Sea with a ratio of about 2. Therefore, we decided to replace ERA5 with ERA-Interim precipitation for the rest of model runs, keeping the source for other variables the same. The extension of the Black Sea to the Marmara "box" through the Bosphorus Strait aimed to represent strait and basin coupling in the model. This is achieved through relaxation of properties to daily climatological stratification in the box and Flather boundary conditions that Oddo and Pinardi imply that the sea surface height and normal velocity are adjusted to values specified at the boundary with some radiation of the disturbances, we set these to zero and let the model calculate them in accordance with wave radiation and with the help of upstream boundary conditions for outflowing components. As the temperature and salinity near the boundary (effectively in the greater part of the small box) are relaxed to climatological boundary conditions, the outflow and inflow at the artificial boundaries are verified to approach approximate volume conservation apart from a small inconsistency but ensuring twoway exchange at the Bosphorus approaching realistic values with a net barotropic component. The various historical estimates of Bosphorus fluxes based on mass budgets and direct observations in the literature, reevaluated by Özsoy and Altıok (2016a), are listed in Table 2, compared with the model-computed values in the last column.
It is noted in Fig. 6 that the longer-term model fluxes of the upper layer, lower layer and net flows in magnitude and direction are in better agreement with the observations, while the short-term components on the order of few days to a week have greater amplitude in the observations compared to the model, possibly because of the insufficient model resolution to fully duplicate hydraulically controlled, turbulent flow dynamics of the strait at this resolution, as remarked earlier in the model configuration section.

Comparison with Argo floats
Model-generated temperature and salinity profiles were compared with the available Argo float data in the Black Sea. For this purpose, the model results were sampled at each Argo float location and depth within each of the polygons shown in Fig. 1. Figure 7 shows profiles of salinity (left column) and temperature (right column) averaged over the three sub-   Ünlüata et al. (1990). (2) Özsoy et al. (1996); Özsoy and Ünlüata (1998).
(3) Altıok and Kayışoglu (2015). (4) (Jarosz et al., 2011b, a). regions. To construct this plot, a total of 8167 profiles were used for "east" region, 6230 profiles for the "west" region and 5820 profiles for the "rim" region. Model results are shown in color shading in the background, and Argo samples are contoured in black. There is good correlation between model results and Argo salinity for each of the three subregions, indicating that the model is able to reproduce salt conservation under the created mass fluxes and mixing events during the model integration period. The red line shows the 21 salinity contours obtained by the model. Although model results closely follow observations, with amplitude of the mixing events coinciding remarkably well, it seems that there is a slight increase of the depth of the 21 salinity profiles especially in the eastern subregion after 2017. The evolution of temperature in the model versus observational data shown on the right-hand side of Fig. 7 shows good agreement over a decade-long integration.
The model performance, evaluated on its ability to reproduce CIL water mass and in capturing the thermocline depth, seems also satisfactory. The model temperature is shown in the background color, while contour lines show incremented values of 7 to 8.5 • C for Argo (black line) and the model (red line). It is seen that the model is able to generate the CIL water for each of the investigated subregions. The model is also able to reproduce the interannual variation of the thermocline. The contour at the upper 50 m marks the 15 • C temperature (red line for the model; black line for Argo) for reference.
In Fig. 8, the time series of Argo float samples and model temperature at 25 m depth are compared for those samples in the east, west and rim regions. Considerably good visual agreement between the observed and modeled temperatures suggests that the surface heat flux in the model and mixed layer dynamics is appropriate for maintaining the surface temperature close to reality over the last decade. Figure 7. CTD profiles of salinity (left column) and temperature (right column) for the three subregions shown in Fig. 1: (a) "east" region, (b) "west" region and (c) "rim" region. Model results were shown as the background fill image for each figure. Argo and model results were also shown by a contour line (17 to 21 for Argo salinity; 21 for model salinity; 7 to 8.5 and 15 • C for the temperature), with the black line for the Argo data and the red line for the model. In Fig. 9, both the observed and the modeled salinity time series at 25 m depth display increasing trends over the 10year period, while the model trend appears slightly larger than the observed one in all regions but especially in the east. These trends are evidence that increased environmental change are in due course in the Black Sea, as has been recently shown by Stanev et al. (2017). We also note that the trend in model salinity has been significantly reduced by choosing ERA-Interim precipitation in place of ERA5 since the latter source was found to result in much smaller total precipitation and had produced a higher salinity trend as well as a much reduced net transport through the Bosphorus (not shown), as discussed in the earlier sections.

Comparison with the tide gauges
Model-generated sea surface height (SSH) was compared with the six tide gauges located along the western and southern Black Sea coasts. For this purpose, SSH anomaly was computed for model and observations, given in Fig. 10 for the run period. The model results are in general agreement with the observations especially at higher frequencies, ex-cept for the irregular seasonal variation of greater amplitude evident in the observations.
It is yet premature to explain this difference at present, given the possible deficit in closing the Black Sea water budget, e.g., by imposing climatological river discharge rather than actual inaccuracies in precipitation and evaporation estimates, as well as model approximations of the free surface. Especially in regard to the latter, first, we note that we have used the linear free surface assumption in the present runs, with plans to use the nonlinear free surface option in future runs. Secondly, the seasonal, spatially variable steric component of sea level in the Black Sea, reported to have an amplitude on the order of 2 cm by Tsimplis et al. (2004) but estimated to reach 6 cm at some locations by Arkhipkin and Berezhnoi (1995), often is not properly accounted for by ocean circulation models at present. All these factors could have contributed to the difference in modeled and observed SSH anomalies.
Samples of annual comparisons between modeled and observed SSH anomalies in Fig. 11 are more optimistic, however, as they display better agreement between the higherfrequency components in the range of daily to weekly and  fortnightly oscillations, leaving aside the longer-term and seasonal, possibly steric or freshwater-induced, components.

Circulation and mixing
The main characteristics of the Black Sea general circulation have been rather well established based on a multitude of experimental investigations, numerical modeling studies and reviews (e.g., Stanev, 1990;Oguz et al., 1993;Oguz and Malanotte-Rizzoli, 1996;Ünlüata, 1997, 1998;Sarkisyan and Sündermann, 2009). Our purpose is not to provide an exposé of the seasonal circulation expected but rather to briefly demonstrate model realizations of well-known circulation features through some examples.
In Fig. 12, we present an exemplary selection of monthly average circulation and temperature fields at the surface (left) and at depth of 50 m (right) developed from December 2011 to March 2012. It can be noticed by from the plots on the right-hand side that downwelling exists all along the continental slope region of the closed basin of the Black Sea, where higher temperatures near the periphery occur and  where the cyclonic rim current also roughly resides. In fact, a well-developed cyclonic surface circulation with meandering currents and anticyclonic eddies along the periphery of the deep basin is evident. On the left-hand side, the greater temperatures in December are reduced until March, where a cold patch surviving in the north central part of the western basin is advected to the southeastern basin near the Bosphorus, while at the same time cold water is formed in the shallow northeast shelf and the Sea of Azov, preserved there until March. The cold water on the shelf is later seen in February and March to leak out of the northwest shelf area to appear on the slope of the Danube Delta at 50 m depth and circulated around the basin. Well-noted circulation elements of the Sevastopol eddy trapped west of the Crimean Peninsula and Batumi eddy in the southeastern corner are in development.
In Fig. 13, we present an exemplary selection of monthly average surface currents for the months of July, September 2011, February, May 20112 (left) and July, September 2014, February andMay 2015 (right), where rapid changes occur in terms of the strength of the rim current, the creation and destruction of eddies, the propagation of eddies alon the rim current and across the basin, filaments of jets separating from the main rim current and their roles in weakening and re-organization of ordered fast currents following the periphery, which show a series of dynamic events that rapidly change the structure and sequence of circulation elements that are all too important in the short-and long-term changes in the environment under the evolving climate.

Time evolution at stations
We examine the temperature and salinity evolution in the system. In particular, we examine in Fig. 14 the upper-ocean properties at stations 1 and 2 of Fig. 1, to observe the typical seasonal/interannual response pattern of the forced system started from the initial condition during the relatively short period of about 10 years.
At station 1, at the very center of the basin, in Fig. 14a, we see cycles of summer warming that develops thermal stratification within the mixed layer region, up to a depth of about 25-30 m, followed by wintertime convective mixing. The CIL, conveniently defined as the layer of intermediate   depth water of temperature less than 8 • C, is present at depths of up to 80 m below the mixed layer in the initial conditions imposed on 1 July 2008. The CIL survives until the summer of 2012 and disappears afterwards at this station. While convection seems to occur every winter, fresh cold intermediate water contributing to the core with T < 8 • C principally occurs only in the winter of 2012, in fact, when cold water with T < 7 • C evidently creating convective overturning extends from the surface down to the core of the CIL at about 40 m. In the following years, the cooling seems to be reduced, with new cold water formed but evidently not as cold as the CIL has been characterized in the not-so-distant past, and intermediate waters in later years undergo warming. The salinity time series in Fig. 14a shows the halocline at depths of around 50 m, at just about the CIL lower depth limit. Both the temperature and the salinity records show many finescale features and oscillations associated with the eddying and small-scale motions.
At station 2 (Fig. 1), near the southern boundary of the basin, in Fig. 14b, we observe that the thermal and salinity stratification are much deeper than central basin station 1, evident by the increased depth of the halocline to about 80 m, and the depth of the mixed layer is also deeper compared to the central region. The deeper structure near the coast, in comparison to the central area of doming in the Black Sea, is also evident in Figs. 12 and 13. The CIL is thicker at the nearcoastal station (station 2), its lower limit reaching a depth of about 100 m in this mainly anticyclonic region. The convective overturning in the winter of 2012 therefore also becomes deeper, reaching the lower depths of the CIL core. The increased oscillatory behavior as compared to station 1 is a result of the mesoscale eddying motions of the rim current effective at the periphery of the basin.
The abundance of CIL in the initial conditions maintained for the first 2 years contrasts with the single event of cold intermediate water formation in 2012 and the weaker event in 2017. The reduced convection events in recent years, both in the deeper central basin and near the coast stand, as evidence that great changes are occurring in the Black Sea, very likely to be an amplified response to climate change in the isolated Black Sea basin that is severely limited in its communication with the Mediterranean Sea and eventually with the World Ocean.

Conclusions
In this study, a high-resolution Black Sea model was developed, with additional merits of coupling the basin hydrodynamics with exchange flows at Bosphorus Strait, with an artificial box on the Marmara side where the temperature and salinity are relaxed to climatology on a daily basis and open boundary conditions allowing a net flow in tune with the net flow across the strait and allowing two-layer outflows and inflows of the Marmara Sea. While previous modeling studies mostly used climatological inflow/outflow water transport specified at the strait entrance as a riverine-like input over the limited number of grid points, this study includes the highfrequency variations in the exchange flow while still satisfying the net water balance of the Black Sea achieved by extending the model domain to include part of the Marmara Sea and the Bosphorus Strait. The model results were rigorously evaluated by comparing them with the huge number of Argo float profiles and available tide gauge station observations in the Black Sea.
The detailed analysis of the monthly mean circulation revealed complex and rapidly changing eddy activity in the Black Sea. Model temperature and salinity results are in good agreement with observations during the investigated period. Slight increases of salinity detected by observations in the upper water column are also evident in the model results.
Overall, 10 years of simulation of temperature and salinity do not show any significant drift, which is an important feature that would allow the model to be used for long-term climate change simulations. The main features of the Black Sea oceanography such as CIL water mass, rim current and upwelling along the southern coast are all well detected by the newly developed model. The model was also able to reproduce the SSH variability in comparison to the tide gauge observations.
Although high-frequency sea-level variations seem to be well reproduced by the model, there appears some mismatch between the model and observation in the SSH estimates at seasonal timescales, generally amplified in the summer season. The differences in seasonal SSH pattern could be due to some deficits in water budget based on climatological river fluxes and steric effects often not considered in models.
In the future, the Marmara box domain should be extended to include the north Aegean Sea, which will enable specification of the SSH and baroclinic velocities from Copernicus Marine Environment Monitoring Service (CMEMS)-Mediterranean model at the open boundaries. It is expected that this will improve performance of the model. The current model has used the linear free surface option. Since we have only defined temperature and salinity over limited part of the Marmara Sea, the nonlinear free surface option was found to generate higher velocities and instabilities at strait entrances, eventually causing the model to blow up. Extending the model domain, with proper open boundary conditions, it is believed that the model could produce better re-sults in nonlinear free surface mode. The behavior of the model should also be investigated by conducting long-term model runs, starting with initial conditions from 1992 to the present to evaluate climate change aspects.
Code and data availability. BSEA is a regional configuration of NEMO at version 4 stable (Madec and the NEMO Team, 2012). Model code is freely available from the NEMO website (https: //www.nemo-ocean.eu/wp-content/uploads/NEMO_book.pdf, NEMO, 2016). After registration, the Fortran code is readily available using open-source subversion software (https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/ html/install.html#download-the-nemo-source-code, Madec and the NEMO Team, 2012). There is only a minor change from the original code. The friction was artificially increased over the Bosphorus Strait grid points.
Data availability. The data can be made available by contacting the authors.
Author contributions. MG, EO and RH set up the NEMO model for the Black Sea. MG and EO performed the experiment. MG, EO and RH analyzed the model results. MG compared the model results with the observations. MG and EO wrote the paper with the inputs from all coauthors.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. Murat Gunduz was supported by the Scientific and Technological Research Council of Turkey (TUBITAK) (project no. 114Y851). During the course of this study, we have obtained valuable help from various colleagues. An early version of the model configuration based on NEMO version 3.1 had been initially set up by Nathalie Toque; the objective analyses of wholebasin temperature and salinity data to set up Black Sea initial conditions were carried out with the help of Mayotte Patron and Hazem Nagy; help in designing the model grid and smoothing the topography was kindly provided by Adil Sözer. We have made efficient use of computer facilities at our affiliated institutions for the computations and those at ECMWF to obtain surface atmospheric data from MARS archives, in order to construct model forcing. Our special thanks are due to the few functionaries of the establishment that unintentionally strengthened the collaboration that made this study possible despite several breaks.
Review statement. This paper was edited by Robert Marsh and reviewed by three anonymous referees.