**2. 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–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 ( *Rd* = *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 *Cd* = 2 × 10−<sup>3</sup> 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].

**Figure 1.** Mesh and bathymetry of the SHYFEM-TSS model domain. Bathymetry of the Dardanelles (top left) and Bosphorus (bottom right) Straits are shown in detail in the small panels. The red line shows the location of the cross sections used in Figure 6. Northern and southern sections to compute volume transports are shown in black lines at the top left and bottom right bathymetry figures.
