Modeling of the Turkish Strait System Using a High Resolution Unstructured Grid Ocean Circulation Model

The Turkish Strait System, which is the only connection between the Black Sea and the Mediterranean Sea, is a challenging region for ocean circulation models due to topographic constraints and water mass structure. We present a newly developed high resolution unstructured finite element grid model to simulate the Turkish Strait System using realistic atmospheric forcing and lateral open boundary conditions. We find that the jet flowing from the Bosphorus Strait into the Marmara creates an anticyclonic circulation. The eddy kinetic energy field is high around the jets exiting from the Bosphorus Strait, Dardanelles Strait, and also the leeward side of the islands in the Marmara Sea. The model successfully captures the two-layer structure of the Sea of Marmara. The volume transport at the Bosphorus is around 120 km3/year which is consistent with the recent observations. The largest bias in the model is at the interface depth due to the shallower mixed layer.


Introduction
The Turkish Strait System (TSS) consisting of the Marmara Sea, the Bosphorus, and the Dardanelles Strait is the only connection between two fundamentally different marginal seas: the Black Sea and the Mediterranean Sea. While the Black Sea has the dynamics of a typical estuarine circulation i.e., precipitation plus river input surpass evaporation (E − P − R < 0) [1,2], the Mediterranean Sea is an example of an inverse-estuary (E − P − R > 0) [3,4]. The TSS has unique dynamics since the flow is both density-driven between the salty (≈38 psu) Mediterranean Sea and the brackish (≈17 psu) Black Sea, and also barotropic-driven because of the permanent sea level difference between these two marginal seas. The sea level height in the Black Sea is approximately 0.3 m higher than the Marmara Sea due to its lower density. Climatological northeasterly winds weakly influence the general circulation in the Sea of Marmara [5].
The Marmara Sea is a relatively small inland sea which covers an area of 11,500 km 2 with three deep basins (>1000 m) and an extended shelf in the south. The Bosphorus and the Dardanelles Straits share common physical and geographical properties. They are both relatively narrow (Bosphorus 0.7-3.5 km, the Dardanelles 1.2-7 km), long (31 km for the Bosphorus and 61 km for the Dardanelles), and shallow (30-100 m for the Bosphorus and 50-120 m for the Dardanelles) channels [6][7][8]. Both Straits enable a water exchange between two different marginal seas with significant density differences.
The fresh water inflow from the Marmara Sea into the Aegean is crucial for the northern part of the Aegean Sea. Convection in the northern part of the Aegean is regulated by the brackish waters of the Black Sea coming from the Dardanelles. In the Black Sea, the river input is balanced by the salty and warm flow coming from the Mediterranean Basin through the Bosphorus Strait. This salty dense water (with the Cold Intermediate Water) is also important for the ventilation process of the deep, oxygen depleted Black Sea waters. Salty Mediterranean Sea waters are equilibrated at the depth of the suboxic and anoxic layers of the Black Sea and play an important role in the redox potential of the chemistry of the Black Sea [8].
Most ocean circulation models use structured grids with finite difference/volume discretizations. These types of grids are particularly challenging in modeling the TSS system because of the narrow and long Bosphorus and Dardanelles Straits, and the deep basins of the Marmara Sea. The complex topography and two-layer density structure require multiple hydraulic controls in both Straits [9]. These narrow Straits thus require a high resolution (less than 100 m) to resolve their dynamics. In addition, the TSS salinity and heat balances are controlled by the Black Sea and the Mediterranean Sea. Hence, the entire TSS region needs to be part of the modeling effort.
Some studies have modeled the transports of individual straits using either 2D idealized reduced gravity or non-hydrostatic models or 3D regional idealized models [10][11][12][13][14]. Some of these studies have attempted to model the Marmara Sea with or without the Turkish Straits using structured grid ocean models. Demyshev et al. conducted a finite difference numerical simulation without atmospheric forcing and reproduced the S-shaped circulation of the jet current exiting the Bosphorus and crossing the Marmara Sea with a basin scale anti-cyclonic circulation [15]. However, the model does not contain strait dynamics. Similarly, the Regional Ocean Modeling System (ROMS [16]) is used to model the Marmara Sea using realistic atmospheric forcing but with open boundaries at the Straits nudged to observation fields [17]. ROMS is a terrain following a structured grid model which needs bathymetry smoothing due to pressure gradient error. Without including the Straits in the model, the authors assumed that the surface circulation depends solely on the strength and directional pattern of the wind force in the Sea of Marmara. Sannino et al. used the structured grid MITgcm model [18] with curvilinear coordinates for a high resolution around the Bosphorus Strait [19]. Their model also did not use any atmospheric forcing. They investigated the circulation of the TSS changing the barotropic flow between the Black Sea and the Aegean Sea.
Few studies have used unstructured grid models to simulate the three-dimensional baroclinic circulation of the TSS. Stanev et al. used the SCHISM model for interconnected basins including the Black Sea, the TSS, and the northern Aegean Sea [20]. They used realistic atmospheric forcing with lateral open boundaries in the south. Their aim was to accurately represent the transport at the straits and the resulting circulation dynamics of the Black Sea. An implicit advection scheme was used in SCHISM for larger time steps. This leads to a coarser model resolution, with only 53 vertical levels at the deepest point of the Black Sea.
Similarly, the Finite Element Sea-Ice Ocean Model (FESOM), another unstructured grid model, was used by Aydogdu et al. to study the circulation of the TSS [21]. They analyzed the combined response of the Sea of Marmara with atmospheric forcing and strait dynamics. Although the FESOM model has a high resolution of up to 65 m in the horizontal and 110 vertical levels, the setup has a closed boundary, and a volume correction is needed to maintain the sea level difference between the Black Sea and the Aegean Sea sides.
Our aim in this study is to simulate the entire TSS region using an unstructured grid model with high resolution in both horizontal and vertical directions. The model has been forced by realistic atmospheric reanalysis and open boundaries with ocean analysis data sets. The main goal is to develop a regional model with an adequate representation of the mean and variability of the TSS circulation. We plan to use the output from the new model as lateral boundary conditions for future Copernicus Marine Service (CMEMS) Mediterranean Sea and Black Sea models. This paper is organized as follows: the numerical model and details of the experiment are introduced in Section 2. The main results including mean circulation properties, validation of water mass structure, and volume fluxes across the straits are presented in Section 3. Finally, we summarize and conclude in Section 4.

Model Setup
We used the System of HydrodYnamic Finite Element Modules (SHYFEM) as the numerical model. SHYFEM is a free-surface, hydrostatic, primitive equations ocean model which uses an unstructured grid finite element scheme in the horizontal and vertically layered in the vertical [22][23][24]. A key advantage of an unstructured grid is to obtain a varying resolution which allows for a finer mesh in the coastal areas, which is also sufficiently accurate in the open ocean and a comparable resolution with parent models. The model domain covers the area between longitudes 22.54°E and 31.41°E, latitudes 38.79°N and 43.41°N (see Figure 1). In this new model, there is a high resolution mesh around 50 m in the shallow areas to resolve the Turkish Straits and 2500 m in the deep areas.
Given that the Rossby radius of deformation (R d = g∆ρ ρ H/ f ) is around 18 km in the Sea of Marmara, our new model can be considered as an eddy resolving model. The model has 93 geopotential coordinate levels in the vertical in order to resolve the complex hydraulics of the rapidly changing bathymetry of the Marmara Sea. The vertical resolution is 1 m in the first 50 m of depth and increases to 100 m at the bottom boundary layer in the deepest part of the model domain.
We have used the high resolution bathymetry and initial conditions provided by Aydogdu et al. [21]. We used April temperature and salinity profiles since they had a minimum bias compared to observations. A Smagorinsky-type dynamical momentum closure scheme with a non-dimensional constant of 2.2 was used in the horizontal closure for momentum. Scale aware momentum closures, especially Smagorinsky, significantly reduce the numerical mixing while keeping the energetics of the mesoscale eddies [25]. For the vertical mixing scheme, a k-epsilon type second order turbulence closure with the Canuto-A stability function was chosen. This scheme performs significantly better in density driven flows (i.e., gravity currents) than the other closures such as k-omega and K-Profile Parameterization [26].
A nonlinear equation of state is used to compute the density [27]. A quadratic bottom drag formulation with a drag coefficient of C d = 2 × 10 −3 is used to compute the bottom friction. No slip boundary conditions are applied in the horizontal except open boundaries. A total variation diminishing (TVD) scheme is used for tracer advection to ensure conservation properties.
We conducted a four-year simulation between 1 January 2016 and 31 December 2019. Atmospheric surface boundary conditions were provided by the European Center for Medium-Range Weather Forecasts (ECMWF) with 1/8°resolution. The forcing data were applied using bulk formulae with a frequency of 6 h. Clamped type open boundary conditions were employed at the southern and northeast boundaries for sea surface height and inflow active tracers. Total velocities were nudged at the open boundaries and zero gradient boundary conditions were used for outflow active tracers. We used daily averaged values of sea surface height, u-velocity, v-velocity, temperature, and salinity fields for the open boundaries. These data were provided by a CMEMS (https://marine.copernicus.eu/, analysis, in particular the Black Sea Forecasting System [28] for the northern boundary and Mediterranean Sea Forecasting System for the southern boundary [29].

Results
This section presents the results of the model from the SHYFEM-TSS setup between 2016 and 2019, including detailed analyses of the surface circulation, vertical water mass structure, and validation against the observations. All the time mean results used here are averaged over the whole period of the model simulation unless stated otherwise. Time series are presented of the volume transports at the two gateways of the Marmara Seanorthern and southern parts of the Bosphorus Strait, and northern and southern parts of the Dardanelles Strait. Figure 2 shows the mean of the surface speed averaged over the simulation period. The Bosphorus Jet reaches the southern coast of the Marmara Sea (Bozburun Peninsula) and turns west towards the middle of the basin. The mean speed of the jet is around 0.6 m/s at the exit of the Bosphorus into the Marmara Sea. The jet separates into two branches due to the Imrali island at the west of the Bozburun Peninsula. The upper branch, which is also the main branch, also splits into two different branches. One part of the upper branch forms an anticyclonic gyre in the northern part of the Marmara Sea, whereas the second part crosses the Marmara Sea after turning south and exits the Marmara Sea through the Dardanelles Strait. The lower branch splits at the Imrali Island, circulates in between the Marmara islands and merges with the upper branch to exit the Marmara Sea.

Mean Surface Circulation and Water Mass Structure
Previous studies have shown similar circulation patterns [5,17,21]. One of the numerical studies which was lacking wind forcing demonstrated a cyclonic gyre in the central basin due to potential vorticity injection from the Bosphorus Jet [19]. This finding is similar to another previous study where a cyclonic circulation has been found in the Arctic Ocean using an idealized barotropic model [30]. However, the mean state in our model, with a more realistic forcing, is different from those idealized simulations. The mean circulation pattern in the central Marmara Sea was persistent with different strengths throughout the simulation. The mean kinetic energy (MKE) and eddy kinetic energy (EKE) fields are shown in Figure 3. We decompose the model velocity into low-pass mean (U) and fluctuation values (u ) with respect to this low-pass mean (i.e., u = U + u ), and compute MKE as 0.5 × (U 2 + V 2 ) and EKE as 0.5 × (u 2 + v 2 ). We used a 60 day low pass filter, which is approximately the mean residence time at the surface layer [5], to compute U. The mean kinetic energy field presents intensities on the rim current in the Black Sea, the Bosphorus Jet, the current out of Dardanelles, and the cyclonic gyre in the Marmara Sea. These regions are also evident in the mean speed field in the previous figure. The EKE, which represents the mesoscale eddy kinetic energy field, differs from the MKE, especially in the Marmara Sea and the northwest of the Aegean Sea. The magnitude of the EKE field is comparable to that of the MKE field. In the Marmara Sea, the EKE field is high around the Bosphorus Jet, Izmit Bay, and southern shelf. MKE decreases and EKE increases downstream of the Bosphorus Jet where it starts to turn to the west in the Marmara Sea. This indication of energy for the EKE field is provided by the release of energy of the mean flow through baroclinic instability processes. There is another high EKE region in the northwestern part of the Aegean Sea exit of the Black Sea water from the Dardanelles. Once again, there is an energy conversion from the MKE field into the EKE field. The sea surface temperature (SST) averaged over four years ranges from 15 • C at the northern part of the Black Sea to 21 • C at the eastern part of the Aegean Sea (Figure 4). The cold jet flow coming from the Bosphorus can be easily traced to the southern coast of the Marmara Sea. The jet is carrying the Black Sea water, approximately 16 • C, into the Marmara at the surface. This is consistent with the mean circulation mentioned above. One of the interesting features in the simulation is the eastern part of the Aegean Sea, outside of the Dardanelles Strait and Saroz Bay. The cold water upwelling process in the summer months dominates the averaged SST field in these regions.   Figure 6 shows the vertical section of time averaged salinity and temperature fields across the section shown in Figure 1 (red line). This cross section follows the deepest part of the Bosphorus and Dardanelles Straits and also passes through three deep basins of the Marmara Sea. The interface depth between the fresh-warm surface layer and salty-cold lower layer is at around 25 m in the central basin of the Marmara Sea. The salinity below 50 m becomes uniform in the Marmara Sea due to the influence of the Mediterranean Sea. However, during the simulation period, deeper layers (i.e., below 200 m) were getting modified with the relatively cold Aegean Sea water coming from the Dardanelles (temperature field at 150 km in Figure 6). There is a distinct cold intermediate layer (CIL) extending from the Bosphorus Strait into the Marmara Sea at the halocline depth in the model. Previous observations also reported a cold temperature anomaly around 11-12 • C at the same depth in June-July between 1996 and 2000 [13]. We conclude that the model successfully reproduces CIL, which is an important feature for the Marmara Sea water mass structure.
The two-layered water mass structure can also be seen in the Dardanelles Strait (first 100 km of the section shown in Figure 6). The mean salinity and temperature vertical structure are similar to those of Sannino et al. [19] and Aydogdu et al. [21]. In the Bosphorus Strait, salinity in the bottom layer is mixed with the ambient water when the layer flows into the Black Sea at the northern exit (section between 320 km and 400 km shown in Figure 6). Previous observations and numerical studies reported a similar salinity distribution [9,14,19,21,31]. The flow in the Bosphorus Strait has two distinct hydraulic control locations which increase the amount of mixing between two layers. The cold upper layer of the Black Sea water mixes with the relatively warm bottom layer of the Marmara Sea water and merges into the CIL at the southern side of the Bosphorus Strait. Surface layer velocities are higher at the southern exits of the Dardanelles and Bosphorus Straits compared to their northern exits. This is due to increasing exit widths in the straits and also contraction regions. For more details of the hydraulic controls in the Bosphorus and Dardanelles Straits, see Sozer and Ozsoy [14], Sannino et al., [19] and references therein.  In each cruise, conductivity, temperature, and depth (CTD) measurements were done from the surface up to the sea bottom at the pre-defined same 90 stations at the speed of 16 Hz. The accuracy rates of the sensors were 0.001 • C and 0.0003 S/m for temperature and conductivity, respectively. After quality control, the high-resolution data were averaged to 1 m bin size during the post-processing for later use. Figure 7a shows the spatial distribution of the interpolated sea surface temperature field in August 2017 (summer cruise). The black dots represent the station locations. The interpolation was performed using the Data-Interpolating Variational Analysis (DIVA) method. There is an approximately 3°C temperature gradient between the east and west of the Marmara basin. Relatively cold Black Sea water (≈23°C) has spread to the exit of the Bosphorus Strait (northeast side of the basin). However, the cold water at the coast of the Bozburun Peninsula (southeastern side of the basin) is probably due to an upwelling event, since the salinity at the same location is much higher than the Black Sea brackish water (Figure 8a). Observations show that the western part of the Marmara Sea is around 26-27°C. Figure 7b shows the model SST field in August 2017. Cold jet water can be seen at the eastern side of the domain. One of the prominent features of the model is the cold upwelling at the Saroz Bay, north of the Dardanelles. The northern coast of the Marmara Sea is one degree warmer in the model than in the observation.  Vertical salinity/temperature bias and root mean square (rms) error profiles compared to all observations are shown in Figures 9 and 10, respectively. The mean surface salinity bias improves over time. The largest bias in all seasons is around the mixed layer depth (20-30 m). This is probably due to the missing representation of mixed layer dynamics such as insufficient mixing. The surface temperature mean bias is around 0.1°C in the top 10 m in August 2018. Similarly to salinity, temperature bias and rms error are highest around the seasonal thermocline depth. We believe that the bias below 40 m is due to the initial conditions. The sea surface temperature at the interior of the Marmara Sea is close to the observed values. Figure 8a shows the sea surface salinity field from the observations. The surface salinity clearly shows a wide spread between 18 psu in the Black Sea entrance of the Bosphorus and 25 psu at the Dardanelles exit of the Marmara Sea. Once again, there is an east-western gradient in salinity in the basin. The exit of the Bosphorus and its surroundings have less than 21 psu of salinity, while the western side surface salinity is more than 24 psu. The low salinity jet flow can also be seen in the model simulation (Figure 8b). Interior salinity values of the model are similar to those in the observed field. The southern coast of the Marmara Sea is slightly saltier in the model due to the upwelling of the saltier deep waters to the surface. This is also likely due to relatively colder waters in the SST field.     Figure 11 shows the daily and monthly averaged net, upper and lower layers' volume fluxes at different exits of the Bosphorus and Dardanelles Straits. The means of each time series are also shown in Table 1. Positive values indicate flows towards the Black/Marmara Sea for the Bosphorus/Dardanelles sections. Significant temporal variability can be seen in all sections. Daily fluctuations may be three times more than the mean value of the corresponding time series. The upper layer variability is higher than the lower layer along the Bosphorus Strait sections, while layer fluctuations are comparable in the Dardanelles. The net flow at the Bosphorus Strait reverses in the winter because of weakening and/or reversing of the upper layer. The lower layer at the Bosphorus Strait is at its minimum in the summer. Field observations show similar fluctuation differences and a reverse of the upper layer in the Bosphorus Strait [32]. The net volume transports in our model are around 140 km 3 yr −1 and 170 km 3 yr −1 at the Bosphorus Strait and the Dardanelles Strait, respectively. The difference may be due to the evaporation-precipitation (E-P) balance in the Sea of Marmara and numerical errors during the post-processing of the unstructured grid. Nevertheless, the net values are approximately half those of the historical estimation which is around 300 km 3 yr −1 and which is based on the conservation of mass in the Black Sea. However, net transport using direct measurements obtained at the exits of the Bosphorus Strait is around 110-120 km 3 yr −1 [32]. We conclude that model transport values are within the uncertainty of the recent observations and other model studies.

Volume Fluxes through the Straits
The previous studies that obtained the net transport of 300 km 3 yr −1 assumed that maximum salinity values of the deep layer can be used to compute a two-layer budget analysis with (E-P) balance in the Black Sea [5]. Here, we use mean salinity values of the bottom and top layers at the northern and southern sections of the Bosphorus Strait and show that the net transport cannot be that high. The net transport (Q net ) is the difference between the top layer transport (Q 1 ) and bottom layer transport (Q 2 ) at the northern section This is also true for the top layer transport (Q 3 ) and bottom layer transport (Q 4 ) at the southern section (i.e., Q net = Q 3 − Q 4 ). We then follow the two-layer analytical model for the salt transport at each sections, where S i is the average salinity of the respected layer. The last two equations above are equal to each other, thus we get Change in transport of upper layer between sections (∆Q = Q 1 − Q 3 ) is around 50 km 3 yr −1 both in our model and in the analytical solution [5]. In addition, in the model, the averaged salinity differences between the upper and lower layers at a cross section are approximately equal to each other (i.e., ∆S = S 1 − S 2 ≈ S 3 − S 4 ). Hence, If we choose the maximum salinity values in a layer as used in [5], S 1 = 17 psu, S 2 = 35 psu, S 3 = 20 psu and S 4 = 37 psu, then Equation (3) satisfies Q net = 300 km 3 /yr −1 . However, the averaged values in a layer are significantly different to the maximum ones because of the mixing in the interfacial layer. The results of the model show that ∆S is around 8 psu rather than 17 psu. This ensures that Q net should be approximately 150 km 3 /yr −1 .
The layer transports at the Dardanelles are higher than shown by the model results provided by [21] but still lower than the historical estimates which are around 850/550 km 3 /yr −1 for the upper/lower layer in the northern section and 1200/900 km 3 /yr −1 for the upper/lower layer in the southern section, respectively [1,5,6,33]. This disagreement in layer fluxes may be due to an inaccurate representation of turbulent mixing and/or forcing along the open boundaries.

Summary and Conclusions
We have presented a new high resolution unstructured grid model for the Turkish Strait System using realistic atmospheric surface forcing and ocean lateral open boundary conditions. The SHYFEM-TSS model was run for four years and successfully reproduced vertical water mass structure, horizontal circulation, and volume fluxes through the Bosphorus and Dardanelles Straits.
The model captures the buoyant jet coming from the Bosphorus into the Marmara Sea on the surface. There is an anticyclonic gyre in the center of the Marmara Sea due to this jet and its potential vorticity injection into the semi-enclosed sea. There is also another fast current flowing out of the Dardanelles into the Aegean Sea carrying relatively fresh Black Sea water to the Eastern Mediterranean Sea. Previous models also observed similar circulation patterns [5,17,21]. Both mean and eddy kinetic energy fields are high around the jets. The mean kinetic energy also shows the main branch of the flow pattern in the Marmara Sea. The eddy kinetic energy field is high around the islands of the Marmara Sea, indicating the so-called island effect, eddies generated due to topographic features.
The temperature and salinity profiles were compared with the observations, and bias/rms fields were computed. Our model performs well compared to the observation, and reproduces the two-layer structure in the Marmara Sea. The SHYFEM-Marmara captures the realistic distribution of the water mass structure. The largest bias is around the pycnocline depth which indicates less mixing in the model. However, it is also possible that small changes at the interface depth show larger errors especially in a two-layer system.
The volume transports at the Straits were also computed. However, the net Bosphorus transport is smaller than the historical value of 300 km 3 /year estimated based on estuary circulation equations using the maximum salinity of the layers [5], and agrees well with recent observation values 120 km 3 /year [32] measured by multiple bottom-mounted acoustic Doppler current profilers. When we computed the estuary circulation equations using the averaged layer salinity values from the model results, we also get the 150 km 3 /year which is in line with net transport in the model. The upper and lower layer transports in both Straits display significant variability at different time scales.
We have demonstrated that the new SHYFEM-Marmara model using a high resolution unstructured grid with realistic atmospheric forcing and lateral open boundary conditions is capable of simulating the main circulation and water characteristics of the Marmara Sea. The future of Marmara Sea modeling relies on the data-assimilation for forecasting capability and also coupling with a biological model to study the impact of man-made pollution coming from the river discharge into the Sea of Marmara. Data Availability Statement: this work has been used CMEMS data from the Med-MFC and BS-MFC, in particular analysis and forecasting products for the physical component of the Mediterranean Sea and the Black Sea (https://marine.copernicus.eu/ (accessed on 9 July 2021)).