**Appendix A**

Surface fluxes are iteratively computed by using bulk formulae parameterization, as proposed by [11,32], for the Mediterranean Sea and revised for the Black Sea basin for the purposes of the operational forecasting system. They are used to handle operational ECMWF analysis and forecast products at resolution of 3–6 h.

The total heat flux *QT* is computed according to:

$$Q\_T = Q\_s + Q\_b + Q\_h + Q\_c \tag{A1}$$

where *Qs* is the shortwave radiation flux, *Qb* is the net longwave radiation flux, *Qh* is the sensible heat flux, *Qe* is the latent heat flux. Such quantities depend upon the air temperature at 2 m ( *TA*), the sea surface temperature computed by the model itself ( *T*0), the total cloud cover ( *C*), the relative humidity computed from the dew point temperature at 2m(*rh*) and the 10 m wind velocity amplitude ( **<sup>V</sup>***w* ).

The *Qs* is computed by means of the [33]:

$$Q\_s = \begin{cases} \,\,\_2Q\_7 = \begin{cases} \,\,\_2Cr(1 - 0.62C + 0.0019\beta)(1 - \alpha) \,\, \dot{f} \,\, \, \text{C} \ge 0.3\\ \,\,\_2Q\_7(1 - \alpha) \,\, \dot{f} \,\, \text{C} < 0.3 \end{cases} \tag{A2}$$

where *β* is the solar noon altitude in degrees and *α* is the ocean surface albedo. The albedo is computed from [34].

The *Qb* is computed by means of the Brunt-Berliand formula as in [35]:

$$Q\_b = \epsilon \sigma T\_0^4 (0.39 - 0.05\sqrt{\varepsilon\_A})(1 - 0.8C) + 4\epsilon \sigma T\_0^3 (T\_0 - T\_A) \tag{A3}$$

where is the ocean emissivity, *σ* is the Stefan-Boltzman constant, *eA* is the atmospheric vapor pressure [36] given as function of the mixing ration of the air ( *wA*) and mean sea level pressure (*p*):

$$
\sigma\_A = \left(\frac{w\_A}{w\_A + 0.622}\right) p \tag{A4}
$$

In the Mediterranean Sea model, it is computed by [37]. The *Qh* is computed as follows:

$$Q\_{\rm h} = -\rho\_A \mathbb{C}\_P \mathbb{C}\_h |\overline{\nabla\_w}| \left( T\_0 - T\_A \right) \tag{A5}$$

where *ρA* is the density of the moist air, *CP* is the specific heat capacity, *Ch* is the turbulent exchange coefficient for humidity set as a constant and equal to 1.3·10−3.

The *Qe* is computed as follows:

$$Q\_{\mathfrak{c}} = -\rho\_A L\_{\mathfrak{c}} \mathbb{C}\_{\mathfrak{c}} \left| \overline{\mathbf{V}\_w} \right| (q\_0 - q\_A) \tag{A6}$$

where *Le* is the latent heat of vaporization, *Ce* is the turbulent exchange coefficient for temperature set as a constant and equal to 1.5·10−3, *q*0 is the specific humidity saturated at *T*0 while *qA* is the specific humidity of air.

The two constants, *Ch* and *Ce*, are computed according to empirical formulation as suggested by [38] and extensively described in [11].

The momentum flux is given by the wind stress components:

$$\left(\pi^{\lambda}, \pi^{\varphi}\right) = \rho\_A \mathbb{C}\_D |\overline{\nabla\_w}| \left(V\_{\mathbf{x}}, V\_{\mathbf{y}}\right) \tag{A7}$$

where *Vx*, *Vy* are the wind components, while *CD* is the drag coefficient given as a function of wind speed and temperature difference *TA* − *T*0 according to [39].

The bulk formulation for the Black Sea are implemented in NEMO and to be used by the model it requires the following list of atmospheric fields in specific units:


