*2.3. Waves*

The BS-WAV forecasting (BS-WAV NRT, [**?** ]) and reanalysis systems (BS-WAV MY, [**?** ]) are based on the WAM Cycle 6 Black Sea model, which replaced the former Cycle 4.6.2 (operational since April 2017) within CMEMS in December 2020. The wave model describes the ocean surface gravity waves (periods of 1.5–25 s). The regional wave model for the semi-enclosed Black Sea runs in shallow water mode with the same spatial resolution as BS-PHY. WAM calculates the two-dimensional energy density spectrum at each of the 44,699 active model grid points in the frequency and directional space. The solution of the energy balance equation is provided for 24 directional bands at 15◦, starting at 7.5◦ and measured clockwise with respect to true north, and 30 frequencies logarithmically spaced from 0.042 Hz to 0.66 Hz at intervals of Δ*f* / *f* = 0.1. Therefore, the prognostic part of the wave model covers approximately 25–1.5 s. In order to include the higher frequency waves into wave-growth/dissipation processes and the output wave characteristics, a parametric tail is fitted for frequencies above the spectral maximum [**?** ]. A detailed description is given in [**???????** ]. The WAM Cycle 6 used for the Black Sea wave hindcast is an update of the former WAM Cycle 4. The basic physics and numerics are kept in the new release. The source function integration scheme by [**?** ] and the model updates by [**?** ] are incorporated. The wave model performance is discussed in [**???** ].

The driving forces for the wave model are the 10 m wind fields provided by the ECMWF analysis and forecast fields, as for BS-PHY. The initial conditions of the BS-WAV NRT are constrained over successive cycles by including a 24 h hindcast run of the model prior to each forecast. The hindcast applies analyzed wind fields to the wave model so that the model is forced by the best available descriptions of the atmosphere and ocean. This prevents any drifts in the initial conditions of the wave model since the key response in the wave model is to the wind, and the analyzed forcing fields reduce the impact of any systematic drifts in the atmospheric model. The parameterization of the wave growth in the wind input source term is adapted to the driving wind fields.

The WAM model estimates the sea-state-dependent momentum and energy fluxes, and the Stokes Coriolis forcing diagnostics needed in order for it to be coupled to the ocean model [**? ?** ]. The wave-induced processes have a significant impact on the skills of drifter estimations, e.g., [**? ?** ]. WAM cycle 6.0 considers the new extreme wave diagnostics (maximal wave height and wave crest, [**? ?** ]) that are included in the new BS-WAV products. In the new BS-WAV NRT wave-breaking parameterization is considered as well as the time-dependent depth and current fields from the BS-PHY NRT. A novel feature of BS-WAV MY data is that radar altimeter data are assimilated. Besides a significant wave height, the assimilation includes wind speed data. The data measured are assimilated into the wave model fields using an optimal interpolation (OI) scheme [**?** ]. Given the lack of

available systematic in situ observations in the Black Sea, the satellite data add value to the wave simulations.

#### **3. Evaluating the Quality of the BS Products in the Operational Framework and for Monitoring**

Evaluating the quality of BS-MFC products is key to providing users with a reliable service. Quality assurance of the product is coordinated within the CMEMS through a dedicated working group, while BS-MFC implements state-of-the-art metrics and supports the centralized product quality dashboard (https://pqd.mercator-ocean.fr/ accessed on 11 October 2021). BS-MFC product quality is based on GODAE/OceanPredict and MERSEA/MyOcean standards for the evaluation of product accuracy [**?** ].

The following sections provide some examples of pre-qualification and operational capacities for the BS-MFC components.

#### *3.1. Validation of BS-PHY Systems*

Estimated accuracy numbers (EANs) for BS-PHY NRT were computed using the daily means analysis fields and compared with available observations (Table **??**). Root mean square difference (RMSD) and bias were estimated over the period of 2019–2020. **????** show EANs computed in the period of 2019–2020 for temperature and salinity, respectively, using ARGO in situ profiles. The evaluation was performed in specific layers. Table **??** shows the averaged RMSD for temperature, which ranges from a maximum of 2.1 ◦C at 20–30 m to below 0.55 ◦C for depths greater than 75 m. The temperature bias is around zero over the water column. Regarding salinity (Table **??**), the error is approximatively 0.26 PSU in the 5–50 m layer, and then increases to about 0.4 PSU in the halocline. Below 200 m, it is less than 0.1 PSU. Salinity bias is generally slightly negative on the subsurface, but never greater than 0.2 PSU. Regarding sea level, the system has an overall error of about 2.9 cm, while for sea surface temperature, the error is about 0.54 ◦C [**?** ]. In order to highlight how BS-PHY NRT reproduces the main features of the circulation at the basin scale, **????** show surface currents as the mean in 2019 and 2020, respectively.

Both figures show the persistency of the Rim current and submesoscale eddies during the selected period—major eddies in the Western basin, the central cyclonic gyre particularly evident in 2020 and less pronounced in 2019, small coastal eddies, and bifurcation of the Rim along the Crimean peninsula. The full overview of the skills is also operationally available through the CMEMS Product Quality Dashboard https://pqd.mercator-ocean.fr/ (accessed on 11 October 2021). Additionally, BS-PHY uses a regional validation website (https://bsfs.cmcc.it/, accessed on 11 October 2021) to evaluate and monitor the quality of the near-real-time physical products.


**Table 2.** EANs for temperature (◦C) in terms of bias and RMSD for 2019–2020, including number of observations used in the comparison.


**Table 3.** EANs for salinity (PSU) in terms of bias and RMSD for 2019–2020, including number of observations used in the comparison.

**Figure 3.** Surface currents in 2019.

**Figure 4.** Surface currents in 2020.

In the BS-PHY MY, the highest temperature RMSDs are at the seasonal thermocline depths of around 20 m as the values exceed 1.5 ◦C (Figure **??**a). The temperature bias is satisfactory, with a maximum absolute value of about 0.07 ◦C at the surface, and around zero from 60 down to 500 m, then it increases up to 0.12 ◦C below 500 m. Salinity has the largest RMSD on the surface, with a maximum of 0.7 PSU (Figure **??**b). Similar to the temperature, salinity RMSDs are relatively high at the halocline depths. The Hovmöller diagram of temperature RMSD reveals a clear seasonal pattern: the values are low (high) in winter (summer) so that the highest errors are in the seasonal thermocline (Figure **??**a). The Hovmöller diagram of salinity RMSD shows that there is a seasonal signal in the error on the surface: RMSD exceeds 1.5 PSU near the surface before 2008. At depths of 50–150 m, the error decreases up to around 0.8 PSU; however, it is almost zero below 150 m (Figure **??**b).

**Figure 5.** Vertical profiles of the root mean square deviation (left panel), bias (middle panel), and number of observations (right panel) for (**a**) temperature (◦C) and (**b**) salinity (PSU) by comparing the reanalysis results with in situ profilers in the Black Sea domain from 1 January 1993 to 31 December 2019.

**Figure 6.** Hovmöller diagrams of root mean square difference for temperature in ◦C (**a**) and salinity in PSU (**b**), by comparing the BS-PHY MY results with in situ profilers in the Black Sea domain from 1 January 1993 to 31 December 2019.

The Hovmöller diagrams reveal the lack of in situ observations between 1997 and 2003 in the Black Sea, which means that for many years the reanalysis system was only constrained by altimeter observations. In comparison with along-track sea level anomaly observations, BS-PHY MY product had a mean sea level anomaly RMSD of 2.24 cm from 1993 to 2019. Finally, the geostrophic currents estimated from the time-averaged sea level (1993–2019) reveal the Rim current, which surrounds the entire Black Sea and forms a large-scale cyclonic gyre (Figure **??**). Major details are provided in [**?** ].

**Figure 7.** Mean sea-surface height and corresponding mean sea-surface circulation from the BS-PHY MY product evaluated for the period of 1993–2019.

#### *3.2. Validation of BS-BIO Systems*

The quality of the BS-BIO NRT and BS-BIO MY systems was assessed using all the observations available for the Black Sea in existing databases (e.g., CMEMS INS and OC TACs, World Ocean database [**?** ], EMODnet at https://portal.emodnet-physics.eu/ (accessed on 11 October 2021), R/V KNORR at https://www.whoi.edu/multimedia/r-vknorr/, accessed on 11 October 2021) using EANs in the form of BIAS\_log10 between the model (M) and the observation (O), as follows (Table **??**):

$$\text{BIAS}\_{\log 10} = \log 10(\text{M}) - \log 10(\text{O}),\tag{1}$$

**Figure 8.** Sub-regions schematization from [**?** ] used for computing the error statistics for chlorophyll. Regions 1–3 are the "open sea", regions 4–5 are the North-Western Shelf, and region 7 is located at the shelf break and is considered a transition region between eutrophic and oligotrophic conditions.


**Table 4.** BIASlog10 statistics for chlorophyll obtained for BS-BIO NRT, considering the observation minus model prediction pairs for the various regions (1–11) shown in Figure **??** Data are from the L3 CMEMS satellite chlorophyll observations (Table **??**).

In the BS-BIO NRT, these error statistics can only be computed for oxygen, chlorophyll, Kd, and PAR, which are derived from satellites and BG-ARGO. Much of the validation has entailed the assessment of the model's capacity to simulate oxygen. Oxygen integrates the balance of physical and biogeochemical processes such as the formation of the cold intermediate layer, the long-term stability of the main pycnocline, dissolution, photosynthesis, and pelagic and benthic respiration. BS-BIO NRT simulates the main open sea oxycline but overestimates its vertical extension compared to BGC-ARGO data profiles that have a deeper injection of rich oxygen filaments (Figure **??**).

**Figure 9.** BGC-ARGO (**top**) and modelled (**bottom**) vertical profiles of oxygen from 2018 to 2020.

In April 2019, the dissolution of oxygen was overestimated, which is probably due to an underestimation of the SST. As shown in Figure **??**, on the north-western shelf, the model (in orange) seems to capture the main characteristics of the observed seasonal cycle (in blue), with an oxygenation of the water column in winter due to ventilation and surface photosynthesis, and a reduction later in the year in response to decreased dissolution and the intensification of respiration. The model tends to underestimate the oxygen concentration on the surface in winter time (Jan–Mar) (Figure **??** left), and does the reverse in the bottom layer in summer–autumn time (Jun–Dec), except in winter when the water column is not stratified (Figure **??** right). The model underestimates the oxygenation of the water column in winter, while later in the year, the underestimation of the surface values and the overestimation of the bottom oxygen could indicate an underestimation of the level of primary production. The lower difference between surface and bottom values in the model as compared to the observations suggests excessive mixing in the model. From July to September, the model simulates hypoxia (O2 < 63 μmol kg−1) in the northernmost part of the shelf; however, the number of hypoxic records is underestimated.

**Figure 10.** Seasonal cycle between 1992 and 2019 of the modelled (orange) and observed (blue) oxygen concentrations over the north-western shelf (at depths of 15–80 m) on the surface (**left**) and at the bottom (**right**). The number of observations for each month is reported along the *X*-axis.

A comparison of the chlorophyll-a (Chla) simulated by BS-BIO NRT and retrieved from the satellite for characteristic optical regions in the Black Sea, as in Figure **??** [**?** ], is shown in Figure **??** (subplots from a to f are related to some crucial regions as shown in Figure **??**). Subplots in Figure **??** show that the model reproduces the winter–early spring bloom typical of the deep Black Sea (Figure **??**a,c,e) and the higher level of production on the north-western shelf with the presence of several peaks (Figure **??**b,d,f). Table **??** shows the log10 BIAS statistics for chlorophyll as the estimated accuracy number (EAN). Considering regions as in Figure **??**, EAN is around zero, except in the North-West coastal region (region 4) and at the Danube Delta (region 5), strongly impacted by water discharge (e.g., from the Danube, the Dniepr and the Dniester). Dynamics of the bloom are strongly conditioned by the amount of nutrients discharged on the shelf and the model has the tendency to slightly underestimate the process, as shown by the positive EAN. This is explained using climatological averaged values for the nutrients discharged by the rivers instead of NRT values. The quality of the boundary conditions at the interface with the rivers can hamper model performance.

**Figure 11.** *Cont*.

**Figure 11.** Time series of satellite (in blue) and model (in orange) chlorophyll spatial averages computed for optical characteristic regions as highlighted by [**?** ].

#### *3.3. Validation of BS-WAV Systems*

This section describes the quality of the BS-WAV NRT and MY products. Validation is performed by comparing model results with satellite observations, as provided by Sentinel-3a, Sentinel-3b, Cryosat-2, Jason-3, and SARAL/Altika for the period between 01 July 2018 and 30 June 2020, and available in situ observations. The assessment of the corresponding wave hindcast is the best way to understand the validity of the WAM domain model underpinning these products, since the wave analysis-forecast system provided to CMEMS considers surface currents and water-level deviations from BS-PHY NRT. The data are incorporated into the WAM model grid. The growth of errors in the wave forecasts is dominated by growth errors in the forcing fields, which are the 10 m wind fields from the ECMWF IFS (Integrated Forecasting System).

In the new MY product (released to users in December 2020), wave breaking and the assimilation of measured satellite data were taken into account. The required radar altimeter data are available on the public server of AVISO at ftp-access.aviso.altimetry.fr (accessed on 11 October 2021) and CMEMS WAVE TAC, and includes significant wave height and also wind speed. The measured data are assimilated into the wave model fields using an optimal interpolation (OI) scheme. Satellites cross the Black Sea once or twice a day for less than two minutes, so very few measured values are available for assimilation into the wave and wind fields. The currents and water level deviations were not considered for the MYP.

Figure **??** depicts the scatter plots for a comparison between the modelled significant wave height (SWH) (provided by BS-WAV MY) and the one taken from satellites (Jason-1, Jason-2, Jason-3, and combinations thereof). The computed and measured mean values of SWH are generally all located between 0.9 and 1.0 m. The bias is nearly zero for all the merged satellite values, although the values for Jason-2 and Jason-3 show a small underestimation of the wave heights in WAM. The calculated biases of WAM for the different satellites confirm the good agreemen<sup>t</sup> between measurements and model results; in fact, they are in the range between −8 cm (Jason-3) and 3 cm (Jason-1). The RMSD varies from 17 cm (Jason-1) to 26 cm (Jason-3). Another deviation is in the difference in the

simulated and modelled variability, given as the standard deviations in this article. The differences range from about 5 cm (Jason-2) to 6 cm (Jason-3). Note that the measurement errors and noise that our initial quality control has not filtered out can also impact metrics, thus potentially degrading the model's skill level. The wind fields also influence the variability and skills of the model as they force the wind-wave model WAM.

**Figure 12.** Scatter plots showing satellite measurements versus modelled significant wave heights (BS-WAV MY) for 2002–2013 (Jason-1), 2008–2017 (Jason-2), 2016–2018 (Jason-3), and 2002–2018 (all satellites merged). Additionally, the following are presented: the estimated bivariate probability density (colored area), the linear slope-fit regression of modelled and observed wave heights (red line), specific quantiles taken from the empirical cumulative density function (black line), and the diagonal (blue line). Summary statistics and skill scores are also included. R: reference (satellite) data, M: model data.

Figure **??** shows two examples of a comparison between BS-WAV NRT significant wave height (SWH) and Jason-3 satellite measurements. The figure includes the distribution of SWH combined with a track of the Jason-3 satellite (upper panel), and the corresponding time series of measured and modelled SWH along the satellite track (lower panel).

On the left is the ascending path on 30 December 2019 16:53:13–16:54:56 UTC, which touches the area of the maximum wave height of about 6 m (Figure **??** left). The second example describes the conditions for the path on 18 December 2018 21:49:15–21:50:58 UTC, which crosses the area of maximum wave height of around 4 m (Figure **??** left). Both comparisons show a good agreemen<sup>t</sup> between satellite measurements and model results, with a small overestimation by the model for the first example near the maximum SWH. Regarding the quarterly validation procedure for the BS-WAV system, the comparisons between modelled and measured satellite data were analyzed for the two quarters of 2018 and all four quarters of 2019 (Table **??**). Representative results are shown in Figure **??** for the last quarter of 2018 for the four different satellites SARAL/Altika, Sentinel-3a, Jason-3, and Cryosat-2. The statistics in Table **??** show an underestimation of the measured data by the wave model simulations, particularly revealing a negative bias in the period of 2019–2020.

**Figure 13. Left**: distribution of BS-WAV NRT SWH on 30 December 2019 (17:00 UTC) and a Jason-3 satellite track also on 30 December 2019 (16:53:13–16:54:56 UTC). **Right**: distribution of SWH on 18 December 2018 (22:00 UTC) and a Jason-3 satellite track on 18 December 2018 (21:49:15–21:50:58 UTC).

**Table 5.** EANs for the current BS-WAV NRT system (all values in centimeters).


**Figure 14.** Scatter plots as in Figure **??**, but for Q4 2018 and four different satellites.

Figure **??** includes a time series of statistical parameters of the BS-WAV MY product for the complete validation period between 2002 and 2019. On this annual scale, the metrics show that the skill level in all of the statistical parameters varies for the individual satellites. While the values for the bias are always around zero for Jason-1 and Jason-2, the bias for Jason-3 is slightly higher, as shown by the green curve in Figure **??** for the Annual Bias from 2016, and proves the underestimation of 7–10 cm of the measurements by the wave model.

**Figure 15.** Yearly values of significant wave height metrics (using the BS-WAV MY system) for the Black Sea derived from individual and combined satellite measurements.

The very good behavior of all statistical parameters for Jason-1 is notable, as there is a very small RMSD value of about 17 cm, a scatter index of about 17%, a correlation of more than 0.96, and a reduction in variance with values above 0.9. That is of course due to the assimilation of Jason-1 data into the wave fields during the time period that they were available.

The parameters for Jason-2 went up for the RMSD (30 cm) around 2010, and down again after 2013 with a low of 17 cm. For the correlation, this is the other way around. The values went down from 0.94 to 0.88 in 2010, and then up again to around 0.97 after 2013. This was expected since after 2013, Jason-2 data were assimilated into the wave fields. The Jason-3 satellite provided data for a short time only. It started in 2016 and showed the most pronounced deviation between model results and the measurements. The RMSD was about 25 cm, the scatter index 25%, and the correlation 0.92. This also generated less satisfactory values for the combined satellite parameters between 2016 and 2019, which were fairly good for the previous years.

Overall, the skill scores depend to some extent on the number of collocated measurements. Clearly, with a decreasing number of observations, the values for the statistical parameters such as RMSD, bias, and scatter index go up, while those for the reduction of variance and the correlation go down. For 11 days, from 4 February 2012 to 15 February 2012, ADCP data were available at the location Pasha Dere (located at 28.03◦ E, 43.08◦ N)

in the western part of the Black Sea near the Bulgarian coast. This period covers a storm event that occurred between 7 and 9 February 2012. A comparison of the significant wave height and peak period with WAM is shown in Figure **??**. The measured and modelled time series of the significant wave height and peak period show a very good agreement. The skill scores support the findings. The bias for the significant wave height is only about −2 cm, while the RMSD is 37 cm. The corresponding values for the peak period, which depend on the resolution of the model frequency, are also very good, with a bias of 0.1 s and an RMSD of about 1 s. The correlation is very high for both integrated parameters: about 0.97 for the significant wave height and 0.91 for the peak period. The peak of the storm was also simulated quite accurately and perfectly illustrates the capacity of the wave model to represent the arrival of a storm as captured by the observations.

**Figure 16.** Significant wave height and peak period from WAM (using the BS-WAV MY system) and from the ADCP station Pasha Dere from 4 to 15 February 2012.

#### **4. Summary and Conclusions**

The BS-MFC products support users and down-stream services in terms of both sustainable industry and society, by means of accurate datasets, for a number of ocean variables such as temperature, salinity, currents for the physical core, nutrients, phytoplankton, chlorophyll, oxygen for the biogeochemical core, significant wave height and period, Stokes drift, and spectra for the wave core.

Since 2016, the BS-MFC has provided high-quality analyses, 10-day forecasts (as daily and hourly means, including hourly instantaneous for wave fields), and reanalyses (as monthly and daily means from 1993, including hourly predictions for wave fields available from 1979) of the ocean variables to describe the ocean state—both for physics and waves—and biogeochemical processes in the Black Sea region. This is part of the CMEMS product catalogue, available through https://marine.copernicus.eu/ (accessed on 11 October 2021). These systems benefit from constant improvements in the modelling components (e.g., structure, formulation, parameterization) and data assimilation capabilities. BS-MFC has developed reliable interfaces with upstream data providers—ECMWF for atmospheric forcing; CMEMS TACs for access to in situ and satellite observations, including SeaDataNet; EMODnet for historical in situ ones—and the CMEMS Dissemination Unit for the operational delivery of NRT and MY products.

The availability of in situ (i.e., ARGO, mooring stations, radar stations) and satellite observations (i.e., sea surface temperature, sea level anomaly, ocean color, significant wave height, etc.) is fundamental for the evolution of the BS-MFC systems, as reported in [**?** ]. They are used for the performance of validation exercises and data assimilation. Observing capabilities in the Black Sea have significantly improved in recent years, thanks to national and international initiatives such as CMEMS by enforcing Black Sea thematic assembly centers. It is crucial that this effort is maintained in various parts of the sea; for example, in the north-western shelf, which receives vast river inflows, adequate monitoring of physical and biogeochemical variables is needed. BS-MFC models can be used to support the design phase of the observational network and to optimize it by using observing system (simulation) experiments.

In the short term, the BS-MFC implementation phase is going to further improve the overall product offer. The BS-PHY team is building new operational systems based on open boundary conditions at the Marmara Sea to improve the representation of the Bosporus Strait dynamics, including the use of historical data for the Danube River, and increased vertical resolution. BS-BIO plans to improve the data assimilation framework, moving towards ensemble approaches to better represent the uncertainty in the prediction. The BS-WAV team proposes new core models based on WAM Cycle 6.0; the NRT is now one-way offline coupled with hourly means of ocean surface currents and sea levels provided by BS-PHY NRT system. MY now assimilates SWH from SL TAC and, similarly to PHY and BIO, is forced by ECMWF ERA5, with which it will also be coupled in the future.

BS-MFC evolutions are carried out by science-driven activities based on the use of state-of-the-art modelling platforms—NEMO, BAMHBI, and WAM, and data assimilation schemes—devoted to solving physical and biogeochemical processes in an accurate way. For the next generation of BS-MFC systems, research and development efforts will focus on ensemble forecasting, refining the coupling strategy of Ocean-Atmosphere-Waves-Biogeochemistry, developing data assimilation capabilities according to newly-available high-resolution data, developing bio-optical models, and enforcing the optimal interface with the Mediterranean Sea using a high-resolution model for the Marmara Sea, which provides high-quality open boundary conditions for both basins.

**Author Contributions:** S.A.C., M.G. (Marilaure Grégoire) and J.S. designed the work and coordinated the preparation of the contributions in collaboration with A.P., G.C., and the following teams: E.J., E.P., L.L., A.A., S.C. (Sergio Creti'), L.S., D.A., S.C. (Salvatore Causio) for the BS-PHY; L.V., A.C., C.M., E.I. for the BS-BIO; A.B., M.R., G.G. for BS-WAV. Service managemen<sup>t</sup> contribution was provided by R.L., V.M., F.P., N.V. Product quality activities were coordinated by E.P for the overall BS-MFC. M.I., M.G. (Murat Gunduz), M.M., S.M., G.C., N.P. contributed to the scientific discussions for the evolutions of the systems. N.V. and P.A. contributed to the overall project management. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Copernicus Marine Environment Monitoring Service for the Black Sea Monitoring and Forecasting Center (Contract No. 72-CMEMS-MFS-BS).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** This work has been used CMEMS data from the BS-MFC, in particular Black Sea physical, biogeochemical and wave products (analysis, forecast, reanalysis) (https:// marine.copernicus.eu/ (accessed on 11 October 2021)).

**Conflicts of Interest:** The authors declare no conflict of interest.
