*2.2. Hydrodynamical Model*

The hydrodynamical core of the forced system is based on the Nucleus for European Modelling of the Ocean (NEMO, version v4.0 [64]). The NEMO code solves primitive equations (derived through assuming hydrostatic and incompressible approximations), and we used the split-explicit free surface formulation. The vertical discretization consists of 31 unevenly spaced z-levels with partial step, with a stretching factor equal to 5, the maximum stretching at the 25th level and a minimum thickness equal to 5 m. A time step of 150 s is used, together with a barotropic time step equal to 1.5 s. The model was initialised in January 2014 using 3D temperature and salinity climatological fields, as provided in the framework of SeaDataCloud v1 for the Black Sea basin [65]. The model computes air-sea fluxes using bulk formulation as implemented for the Mediterranean Sea [66], designed to handle the ECMWF atmospheric forcing data for computing momentum, heat and water fluxes [67,68].

The model considers 72 rivers, whose discharge is estimated from monthly climatological datasets developed in the framework of the SESAME project [69]. The major rivers (e.g., the Danube, Dniepr, Dniester) have been represented as multiple distributed sources and the remainder as source points. The Danube river is represented accounting for the three main branches—the Chilia, the Sulina and the St. George (Figure 2a). Its discharge distribution is set according to [70], with the Chilia accounting for 52%, Sulina for 20% and the St. George 28%. The impact of the Bosphorus Strait on the Black Sea dynamics is accounted for in terms of a surface boundary condition, considering the barotropic transport, which has been computed to balance the freshwater fluxes on monthly basis [71,72].

**Figure 2.** (**a**) The Danube Delta: the Chilia, the Sulina and St. George arms are labelled. Small grey circles show the active grid points of the numerical grid in proximity to the Danube delta. Black squares illustrate the Danube arms mouth representation in the hydrodynamical model. (**b**) The Bosphorus Box is used in the hydrodynamical model. Small grey circles show the active grid points of the numerical grid in the Bosphorus area. Black squares illustrate the points in which the temperature and salinity solution were relaxed to [73] profiles. Black star represents the location in which the T and S were relaxed as in the squares, but in addition, the surface boundary condition of "inverse river" is applied.

The representation of saltier and warmer Mediterranean waters that enter the Black Sea is performed by damping the model solution in the area of influence of the Bosphorus exit (Figure 2b) to monthly 3D temperature and salinity climatologies, as computed from [73], at hourly frequencies. The NEMO configuration uses the Turbulent Kinetic Energy vertical mixing scheme for computing the vertical diffusivity and viscosity. Both horizontal eddy viscosity and diffusivity are constant over the domain and set equal to 1.2 × 10−<sup>5</sup> m2/s

and 1.2 × 10−<sup>6</sup> m2/s, respectively. The NEMO configuration with main parameters can be found in Zenodo [74].

## *2.3. Wave Model*

The wave model used in this work is the third-generation spectral WaveWatchIII [75] version 5.16, hereafter denoted as WW3. The model solves the wave action density balance equation for wavenumber-direction spectra *<sup>N</sup>*(*k*, *θ*, *x*, *t*), where, *θ*, is the wave direction, *k* is the wavenumber, *x* = (*<sup>x</sup>*, *y*) is the position vector, and *t* is time. The source terms considered in this implementation concern deep water processes. The wind input source term *Sin* represents the momentum and energy transfer from air to ocean waves, the wave dissipation due to white-capping *Sds* and the nonlinear transfer by resonant four-wave interactions *Snl*:

$$S = S\_{\rm in} + S\_{ds} + S\_{nl} \tag{1}$$

WW3 has been implemented following WAM Cycle4 model physics [76]. The propagation scheme used is a third-order scheme (Ultimate Quickest) with the "Garden Sprinkler Effect" alleviation method of spatial averaging. Wind input and dissipation are based on [77], in which the wind input parameterization is adapted from Janssen's quasi-linear theory of wind-wave generation [78,79], following adjustments of [80,81]. Nonlinear wave-wave interactions have been modelled using the Discrete Interaction Approximation (DIA) [82,83].

WW3 is a spectral wave model, and therefore requires a discretization of the wave spectra in bins of frequency and direction. The spectral direction over the full cycle is selected and is divided into 24 sectors of 15◦ width. A first (lowest) frequency of 0.05 Hz was selected, with a frequency increment factor of 1.1. and 30 frequencies in total. The WW3 configuration parameters are provided via Zenodo [74].
