*2.4. Coupling Strategy*

In the last few years, several works [7,8,25,27,29,84] deeply described the major wavecurrents interaction processes. In this work we consider the wave-currents coupling processes implemented in a recent NEMO release (v4.0) to evaluate the impact of those processes on the Black Sea thermodynamic properties for the first time, with a primary focus on the SST, subsurface temperature and tracer advection due to vertical mixing. To a lesser extent, we gained information on the effect on water velocity, even if the exiguous number of available observations limited this analysis.

The reciprocally forced system between the hydrodynamic model (HM) NEMO and the wind-wave model (WM) WW3 is described in this section. WM and HM models share the same horizontal grid and bathymetry and are forced through a reciprocal hourly field exchange. Eight fields are exchanged between the hydrodynamic and wave models. HM sends to WM the surface currents and the Sea Surface Temperature (SST), which is used to compute the air-sea temperature difference. On the other side, WM provides to HM the following fields: significant wave height (hereafter Hs), mean wave period (Tm), wave peak frequency (Fp), Stokes drift at the surface (SD), wave energy (Φoc) and sea-state dependent momentum flux (τoc). An illustration of the coupling mechanism is given in Figure 3.

**Figure 3.** Sketch of the forcing mechanism used in this work. Free-runs are marked with FR, forced experiments are marked with OLCR. The fields provide occurs at a one-hour frequency. In the Wave OLCR numerical experiment, NEMO sends current fields (u, v) to WW3 and SST to a system (the light-orange coloured square in the figure), which computes air-sea temperature differences (ΔT) using 2 m temperature from the atmospheric forcing and sends ΔT to WW3. In the hydrodynamic OLCR, WW3 FR sends the Stokes Drift (SD), sea-state dependent momentum flux (τoc), wave energy to ocean energy flux (Φoc), significant wave height (Hs), mean wave period (Tm) and wave frequency peak (Fp) to NEMO.

Five physical processes related to wave-currents interaction are considered: three reflect the impact of waves on hydrodynamics (Figure 4) and two consider the effect of hydrodynamics on the wavefield:


**Figure 4.** Example of wave fields snapshot used to force the hydrodynamic model at a given time step. (**a**) significant wave height, (**b**) wave energy to ocean turbulence, (**c**) sea-state dependent momentum flux and (**d**) Stokes Drift.

#### 2.4.1. Sea-State Dependent Momentum Flux

In numerical ocean modelling, the momentum fluxes from the atmosphere to the ocean are traditionally calculated from the wind speed provided by an atmospheric model, using a drag coefficient relating the 10-m winds to the surface stress. With this formulation, the surface fluxes are dependent on the local wind speed only [85]. This formulation does not consider that surface stress is dependent on the sea state, that is, how the wave energy is distributed over the frequency range. Another problem related to this approach is that the net momentum flux is not necessarily conserved between the atmosphere and the ocean. For this reason, several authors continue to investigate the impact of the sea-state on wind stress for more than a half-century [25–27,77,80,86–89].

The air-sea momentum flux, or the total wind stress (*<sup>τ</sup>tot*), is the sum of the momentum flux into both surface waves and subsurface currents. A significant proportion of the momentum lost from the atmosphere is taken up by the wind waves as increased wave momentum. The appropriate stress (τoc) (Equation (2)) to use in an Eulerian ocean model is thus the proportion of the momentum flux that is not taken up by the waves (*<sup>τ</sup>in*) plus the momentum lost from the wind-generated waves through dissipation (*<sup>τ</sup>dis*) [90].

$$
\pi\_{\infty} = \pi\_{tot} - \pi\_w \tag{2}
$$

$$
\tau\_w = (\tau\_{in} + \tau\_{dis}) \tag{3}
$$

In a wave model, an important source input is the calculation of the momentum flux from the atmosphere to the ocean, which can be expressed as an integral of the wave variance spectrum multiplied by the wave growth rate (momentum-uptake rate).

$$
\pi\_w = \rho\_w g \int\_0^{2\pi} \int\_0^{\omega} \frac{k}{\omega} \left( S\_{in} + S\_{dis} + S\_{nI} \right) d\omega \, d\Theta \tag{4}
$$

where *ρw* is the water density, *k* is the wavenumber vector, *ω* is the wave absolute frequency (in radians), *Sin*, *Sdis*, *Snl* are respectively the input, dissipation and non-linear wave model source-term and Θ is the wave direction.

As in previous work e.g., [24,25,27,57,84] Equation (2), the atmospheric wind stress was corrected as follows 

$$
\pi\_{0\text{c}} = \pi\_{\text{tot}} \left( \frac{\pi\_w}{\pi\_{\text{tot}}} \right) - \pi\_w \tag{5}
$$

#### 2.4.2. Stokes–Coriolis Force and 3D-Stokes Drift Profile

In [9], it is demonstrated that the Stokes drift yields a force on the mean currents. In a rotating ocean, the along-wave crest velocity component of the wave motion is correlated with the vertical component, thus inducing wave stress on the mean Eulerian currents. This stress is proportional to *f* ·*vs* where *f* is the Coriolis force and *vs* is the Stokes drift. Following the notation in [91], this term is known as Coriolis–Stokes forcing (Equation (6)).

$$\frac{Du}{Dt} = -\frac{1}{\rho}\Delta p + (u + v\_s) \times f \,\mathcal{Z} + \frac{1}{\rho} \frac{\partial \tau}{\partial \mathcal{Z}}\tag{6}$$

where *Z*ˆ is the upward vector and *τ* is the stress.

Computing the Stokes drift profile is expensive as it involves evaluating an integral with the two-dimensional (2D) wave spectrum at every desired vertical level. It is also often impractical or impossible as the full 2D wave spectrum may not be available. Thus, the full Stokes drift velocity profile is commonly replaced by a monochromatic profile, matched to the transport and the surface Stokes velocity e.g., [86,92–96]. This is problematic, as it is clear that the shear under a broad spectrum is much stronger than that of a monochromatic wave of intermediate wavenumber, due to the presence of short waves, whose associated Stokes drift quickly vanishes with depth. The deep Stokes drift profile will

also be stronger than that of a monochromatic wave, as the low-wavenumber components penetrate much deeper.

Staneva et al. 2017 [25] developed an alternative approximate Stokes drift profile, which has a lower mean-square error deviation than the monochromatic profile for all tested spectra. It has a stronger shear in the upper part and does not tend to zero as rapidly as the monochromatic profile in the deeper part. It mimics the effect of a broader spectrum where the low-wavenumber components penetrate deeper than the mean wavenumber component, while the shorter waves (higher wavenumbers) only affect the upper part of the water column. In this study, the computation of 3D Stokes drift profiles is conducted using the approach in [26], which is included in NEMO v4.0. The Stokes Drift contribution to the water velocity components has been integrated into tracers and momentum equations as in [7,97].

#### 2.4.3. Wave Induced Turbulence

Following [11], the effects of breaking waves on upper ocean mixing are explicitly considered by the introduction of the energy flux from waves into the ocean water column. Thus, the transport of turbulent kinetic energy (TKE) must also be introduced. In [11], it is assumed that there is a direct conversion of mechanical energy to turbulent energy at the surface and therefore the turbulent energy flux is assumed to be given by the energy flux from waves to the ocean water column Φoc (Equation (7)), which follows from the dissipation term in the energy balance equation.

$$\Phi\_{\rm oc} = -\rho\_{\rm w} \mathcal{g} \int\_0^{2\pi} \int\_0^{\infty} S\_{ds} \, d\Theta d\omega \tag{7}$$

meaning that the injection of TKE at the surface is given by the dissipation of the wavefield via the sink term in the wave model energy balance equation (usually dominated by wave breaking) and converted into an ocean turbulence source term. In the absence of specific information on the sea state, the energy flux is parametrised in NEMO by considering a typical old wind-sea value. NEMO v4.0 does not include the wave-dependent TKE surface boundary condition, and thus we applied the code modification according to [9].

#### 2.4.4. Doppler Effect and Refraction Due to Currents

When surface currents interact with waves, phenomena such as the Doppler shift or refraction occur. The way the Doppler shift modifies the surface waves depends on the current speed relative to the wave propagation speed; thus, slow propagating waves are mainly affected by currents. The effects of currents and waves can merge constructively, creating single exceptionally large waves (rogue waves), or if waves are in opposition to strong currents, they become shorter and steeper (potentially hazardous for navigation) and they can also break. The propagation velocity in the various phase spaces of waves interacting with currents (Equations (8)–(10)) can be written as follows.

.

*x*

$$
\varepsilon = \varepsilon\_{\mathcal{S}} + \mathcal{U}\varepsilon \tag{8}
$$

$$\dot{k} = -\left(\frac{\partial}{\partial d}\sigma\right)\left(\frac{\partial}{\partial \mathbf{s}}\sigma\right) - K \cdot \frac{\partial}{\partial \mathbf{s}} \mathcal{U}\_{\mathbf{c}} \tag{9}$$

$$\dot{\Theta} = -\frac{1}{k} \left[ \left( \frac{\partial}{\partial d} \sigma \right) \left( \frac{\partial}{\partial m} \sigma \right) - K \cdot \frac{\partial}{\partial m} \mathcal{U}\_c \right] \tag{10}$$

where *cg* is the wave propagation velocity vector, *Uc* is the velocity of the current, *d* is the water depth, *s* and *m* are the directions along and perpendicular, respectively, to the wave direction.

#### 2.4.5. Effects of Air Stability on the Growth Rate of Waves

The difference between the sea surface temperature (SST) and the air temperature affects the stability of the lower atmosphere and thus the wind velocity structure. Tolman, 2002 [98] formulated a stability correction by replacing the wind speed with an effective wind speed so that the wave growth reproduces [99] stable and unstable wave growth curves.

The air-sea temperature difference is used to evaluate a stability parameter, *ST* (Equation (11)), which is written as follows:

$$ST = \frac{h\mathcal{g}}{\mathcal{U}h^2} \frac{T\_a - T\_s}{T\_0} \tag{11}$$

where *Uh* is the wind speed at *h* height, *Ta* is air temperature at *h* height, *Ts* is the surface temperature and *T*0 is the reference temperature. *ST* is used to compute effective wind speed, *Ue*:

$$\mathcal{U}L\_{\varepsilon} = \mathcal{U}\_{10} \left( \frac{c\_0}{1 \pm c\_1 \tanh[c\_2(ST - ST\_0)]} \right)^{\frac{1}{2}} \tag{12}$$

where *U*10 is the wind speed at 10 m, values *c*0, *c*1, *c*2 and *ST*0 are set to 1.4, 0.1, 150 and −0.01, respectively.

#### **3. Numerical Experiments Design and Validation Strategy**

Nine numerical experiments were conducted from January 2014 to December 2019, according to the numerical setups explained in the previous section. Five experiments concern hydrodynamics (labelled H) and 4 examine waves (labelled W). All of the numerical experiments are listed in Table 1.

**Table 1.** List of numerical experiments carried out from 2014 to 2019. H refers to hydrodynamics, W refers to waves.


#### *Validation Strategy and Observational Data*

In this work, the model evaluation was based on GODAE/Oceanpredict [100]. Standard statistics such as BIAS, Root-Mean Squared Error (RMSE), Scatter Index (SI), Pearson's correlation index () and Standard Deviation (SDev) are used to evaluate the performance of the ocean models by comparing the numerical results against observations (in-situ and/or satellite observed data):

$$\text{BIAS} = \frac{1}{n} \sum\_{i=1}^{n} (m\_i - o\_i) \tag{13}$$

$$\text{RMSE} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} (m\_i - o\_i)^2} \tag{14}$$

$$\text{SI} = \frac{\sqrt{\frac{1}{n-1}\sum\_{i=1}^{n} \left(m\_i - o\_i - BIAS\right)^2}}{o} \tag{15}$$

$$\rho = \frac{\frac{1}{n-1} \sum\_{i=1}^{n} (o\_i - \underline{\varrho}) \left(m\_i - \underline{m}\right)}{\sqrt{\frac{1}{n-1} \sum\_{i=1}^{n} (o\_i - \underline{\varrho})^2} \sqrt{\frac{1}{n-1} \sum\_{i=1}^{n} \left(m\_i - \underline{m}\right)^2}}\tag{16}$$

$$\text{SDev} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} (x\_i - \overline{x})^2} \tag{17}$$

where *o* and *m* stand for observed and modelled data, respectively. *x* applies for both *m* or *o*; the overbar over a variable denotes the temporal averaged value derived from the sample of length *n*.

Three sources of data, organised by platform, were used to evaluate the accuracy of the model results and described in the following.


Table 2 summarises data type, producers, variables and reference links of the observations used for validating our results.



Due to the lack of systematic in-situ measurements in the Black Sea (buoys or moorings), wave simulations of Hs are validated by comparing them with the J2 radar altimeter. Satellite Hs measurements [2016–2018] are used to compare the wave model results. Observations for more than 30 min in time or 2 km in geographical space are rejected. Validation was conducted using scatter plots, comparing observed against modelled Hs.

The bivariate probability density function was estimated using a 2D-Gaussian kernel on a squared grid in the variable space provided in [103]. The plots include summary statistics to describe the WW3 skills to evaluate the Hs, such as BIAS, RMSE, scatter index, correlation index, slope and standard error. Salinity (S) and Temperature (T) vertical profiles from the hydrodynamics simulation were validated against ARGO profiles, while SST validation was conducted using L4 satellite data provided in the framework of the Copernicus Marine Environment and Monitoring Service. BIAS and RMSE were computed using numerical results from January 2015 to December 2019. The Hovmöller diagram allows the evaluation of the model capability to reconstruct interannual-annual and seasonal cycles of analysed scalar quantities (e.g., temperature or salinity).

In this work Hovmöller T (S) diagram considered the uppermost 300 m to evaluate the impact of coupling on the temporal variability of water masses formation in the Black Sea, from 2016 to 2019. All the grid points with depth lower than 100 m have been excluded from the computation of the basin mean to avoid contamination from the coastal zones, then, a basin averaged daily mean has been computed.

**Figure 5.** Location of the observations used for validating this work. (**a**) Map of Jason-2 satellite tracks from 2016 to 2018 over the Black Sea, and moorings (red stars) location. (**b**) ARGO floats available between 2015 and 2019.

#### **4. Results and Discussion**

*4.1. Validation of Hydrodynamical Component*

4.1.1. T/S Profiles

Sub-surface temperatures and salinity from HM simulations, as reported in Table 1, were validated using available ARGO profiles from 7.5 m to 1000 m depth as described in Table 2.

Figure 6 shows the relevant metrics for temperature and salinity as given by the set of performed experiments, averaged over the basin for the period 2015–2019. Vertical averaged bias profile for temperature (Figure 6a, left) is negative up to 1000 m and specifically the free-run experiment (H0) is characterised by the largest cold bias in the thermocline layer, up to −0.6 ◦C, which is only slightly reduced when the Stokes–Coriolis term and wave breaking effects are accounted in experiments H1 and H3 respectively. Temperature bias decreases up to −0.25 ◦C in the layer 15–40 m when the whole sets of wave-currents interaction processes are considered (H4) and the wave stress (included in experiment H2) provides the largest positive impact in reducing the temperature cold bias in the basin thermocline. In terms of vertical averaged RMSE profile, the envelope of experiments (Figure 6a, mid) shows quite similar values, with higher error in the thermocline, exceeding 2 ◦C.

**Figure 6.** Domain averaged 2015–2019 validation of Temperature profiles (**a**) and Salinity (**b**) between 7.5 and 1000 m deep. For each image, the left panel shows BIAS and the right panel shows RMSE. Red boxes represent a zoomed area of the RMSE plots.

The modified ocean stress accounting for the stress absorbed/released by waves is the one providing the best skill among the "single process" experiments. H4 provides the best performance with a halving of BIAS and an RMSE reduction of about 3%, due to enhanced mixing which impacts the mixed layer positioning. Figure 6b reports the metrics for salinity: in particular, the envelope of experiments exhibits quite similar BIAS values, with a progressive reduction on the BIAS for the fully-forced experiment (H4), dealing with the lowest error at the subsurface (up to 30 m).

Tables 3 and 4 recap skills for the uppermost 200 m at basin scale for both temperature and salinity, respectively: metrics are proposed as averages over every single year and over the period 2015–2019. On a yearly basis, the mean BIAS for temperature (Table 3) is close to zero in H4, but there is no significant difference in terms of RMSE.

**Table 3.** Statistics evaluated by comparing temperature vertical profile measurements and model results from circulation models from 7.5 m to 200 m deep.


**Table 4.** Statistics evaluated by comparing Salinity vertical profile measurements and model results from circulation models from 7.5 to 200 m deep.


Table 4 summarises the salinity validation results at the uppermost 200 m. Below 100 m depth, the profiles are very closely aligned. Over the given period, H4 shows a reduction of BIAS of about 20% and RMSE of about 6.5% up to 75 m. We speculate that the average error of 0.3 PSU given by the considered experiments, is related to weak representation of the boundary conditions at the Bosphorus Strait which is not properly representing the Mediterranean waters inflow at all.

The difference between H0 and H4 at the seasonal time scale from 2015 to 2019 was also investigated, as summarised in Table 5, and for the five-year average, the H4 experiment revealed consistently better performance in terms of both salinity and temperature BIAS and RMSE than H0. The best temperature performance in RMSE for the forced run occurs during Winter (Spring) with 0.56 ◦C (0.61 ◦C) compared to 0.60 ◦C (0.66 ◦C) for H0. During Summer and Autumn, the hydrodynamical model has reduced efficiency (both for H0 and H4 experiments) with almost double RMSE (1.17 ◦C and 1.18 ◦C, respectively), and no or negligible differences between H0 and H4.


**Table 5.** Seasonal evaluation of RMSE and BIAS for salinity and temperature for H0 and H4 experiments from 2015 to 2019. Bold indicates the best value between H0 and H4. DJF = Winter, MAM = Spring, JJA = Summer, SON = Autumn.
