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 medium-resolution 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.

1 Introduction

1.1 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); Oğuz 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., 2003, 2004; Stanev, 2005; Miladinova et al., 2017, 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 specified 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.

1.2 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 et al., 1995, 1996; Ö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., 1993, 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; Stanev et al., 2004). 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ışoğlu, 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 (Farmer and 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., 2014, 2017; Aydoğdu 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 et al., 1995, 1996, 2001; Ö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 et al., 1991; Ö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, 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.

2 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 high-fidelity 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 $\mathrm{1}/\mathrm{12}{}^{\circ }$ grid and 42 vertical levels. Having a coarser grid, using 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.

2.1 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.

Figure 1 Model domain and bathymetry. The polygons show three subregions of the Black Sea (namely “east”, “west”, “rim”) used to sample the Argo floats. Six tide gauges along the Turkish coast are shown by red dots (İğneada, İstanbul, Şile, Amasra, Sinop and Trabzon). Contours in the range 100–2300 m (contour interval 200 m) correspond to the color scale on the right. Contours in the depth range 0–100 m (contour interval 10 m) are superimposed with different color scale (not shown) and 100 m contour marked in magenta. Two stations (S1 and S2) shown in white dots are used to construct time series plots.

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 $\mathrm{\Delta }x\equiv \mathrm{\Delta }y\equiv \mathrm{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 $rx\mathrm{0}=max(\mathrm{\Delta }h/h)<\mathrm{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.

Figure 2 Model vertical grid spacing Δz (meters) as a function of the vertical coordinate z displayed for (a) maximum depth of 2300 m and (b) in the upper 200 m.

2.2 Parameterization

The numerical values of the horizontal viscosity length and velocity are taken at 1 km and 0.1 m s⁻¹, respectively. The horizontal diffusive length scale and velocity are 2 km and 0.01 m s⁻¹, respectively. Values for vertical eddy viscosity and eddy diffusivity are selected as $A_{v}m=\mathrm{1}.e-\mathrm{7}$ and $A_{v}t=\mathrm{1}.e-\mathrm{8}$ m² s⁻¹. 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 configuration. Bilaplacian diffusivity schemes are used for both tracer and momentum.

2.3 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.

Table 1 Climatological monthly river runoff for major Black Sea rivers (m³ s⁻¹).

The total annual river inflow represented with 11 rivers in the model amounts to 10 570 m³ s⁻¹ (=333 km³ yr⁻¹). 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³ s⁻¹ (=289 km³ yr⁻¹).

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³ s⁻¹ (=436 km³ yr⁻¹) 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.

Figure 3 Interannual variations of mean annual freshwater inflow to the Black Sea – source: GRDC (https://www.bafg.de/GRDC/EN/Home/homepage_node.html, last access: 11 January 2020).

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 seasonal 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.

2.4 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.

Figure 4 Time series of daily climatological salinity (b) and temperature (c) profiles and locations of the CTD stations (a) applied for relaxation at the Marmara box shown in Fig. 1.

2.5 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 Oğuz et al. (1993) and Sur et al. (1994, 1996). The data obtained from a total of 394 hydrography stations at nominal spacing of 20^′ ($=\mathrm{1}/\mathrm{3}{}^{\circ }$) 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.

Figure 5 Model initial conditions of temperature (lhs) and salinity (rhs) at the surface (a, b) and at a depth of 89 m (c, d), based on the objective analyses of data from a total of 394 CTD measurement stations obtained during the July 1992 CoMSBlack cruise.

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.

3 Model evaluation

3.1 Bosphorus fluxes

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 (2008) applied at the Marmara box open boundaries. Although the Flather open boundary conditions 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 two-way exchange at the Bosphorus approaching realistic values with a net barotropic component.

Comparison of the time series for the upper layer, lower layer and net fluxes in the Bosphorus Strait obtained from the model and observations from August 2008 to January 2009 are superposed in Fig. 6. The experimental fluxes are based on Jarosz et al. (2011b) measurements. We observe that the model results are in the same range of mean fluxes obtained from observations, although there are differences in the oscillatory part of the signals. We also observe that the lower layer fluxes are in better agreement between the model and observations, while the upper layer fluxes outwards from the Black Sea are less compliant, most likely because the outflow is directly related to the surface water balance of the Black Sea.

The various historical estimates of Bosphorus fluxes based on mass budgets and direct observations in the literature, re-evaluated 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.

Figure 6 Time series of water transport at Bosphorus Strait. (a) Upper and lower fluxes from the model (red line) and from observations (black line) compared for the August 2008 to January 2009 period. (b) Net water transport from model (red line) and observation (black line). (c) Upper and lower layer fluxes from the model for the entire model run period.

Table 2 Water transport (km³ yr⁻¹) components at the two ends and mid-strait sections of the Bosphorus, based on analyses (1) and observations (2, 3, 4) compared with model results.

(1) Ünlüata et al. (1990). (2) Özsoy et al. (1996); Özsoy and Ünlüata (1998). (3) Altıok and Kayışoğlu (2015). (4) (Jarosz et al., 2011b, a).

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.

3.2 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 subregions. 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.

In Fig. 9, both the observed and the modeled salinity time series at 25 m depth display increasing trends over the 10-year 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.

Figure 8 Time series of model and Argo temperature (^∘C) at 25 m sampled over the three regions shown in Fig. 1: (a) east, (b) west and (c) rim. Black dots indicate data from Argo floats; red dots indicate data from model results.

Figure 9 Time series of model salinity at 25 m sampled over the three regions shown in Fig. 1: (a) east, (b) west and (c) rim. Black dots indicate data from Argo floats; red dots indicate data from model results.

3.3 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, except 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 higher-frequency components in the range of daily to weekly and fortnightly oscillations, leaving aside the longer-term and seasonal, possibly steric or freshwater-induced, components.

Figure 10 Time series of sea surface height anomaly (m) for the three tide gauge stations shown in Fig. 1: (a) Amasra, (b) Şile and (c) İğneada. Model results are indicated by the black line and the tide gauge record is indicated by the red line.

Figure 11 Annual time series of sea surface height anomaly (m) for stations (a) İğneada for the year 2014 and (b) Amasra for the year 2017. Model results are indicated by the black line and the tide gauge record is indicated by the red line.

4 Initial results and applications

4.1 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; Oğuz et al., 1993; Oğuz and Malanotte-Rizzoli, 1996; Özsoy and Ü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.

Figure 12 Monthly average currents and temperature at the surface (a) and at a depth of 50 m (b) developed from December 2011 to March 2012 (ordered down the page).

Figure 13 Monthly average current vectors and magnitude in color shading for the months of July, September 2011, February, May 2012 (a) and July, September 2014, February, May 2015 (b) (ordered down the page).

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 and May 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.

Figure 14 Time series of temperature and salinity versus depth at (a) station 1 and (b) station 2 of Fig. 1.

4.2 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 fine-scale 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 near-coastal 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.

5 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 high-frequency 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 results 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.