The Black Sea Physics Analysis and Forecasting System within the Framework of the Copernicus Marine Service

: This work describes the design, implementation and validation of the Black Sea physics analysis and forecasting system, developed by the Black Sea Physics production unit within the Black Sea Monitoring and Forecasting Center as part of the Copernicus Marine Environment and Monitoring Service. The system provides analyses and forecasts of the temperature, salinity, sea surface height, mixed layer depth and currents for the whole Black Sea basin, excluding the Azov Sea, and has been operational since 2016. The system is composed of the NEMO (v 3.4) numerical model and an OceanVar scheme, which brings together real time observations (in-situ temperature and salinity proﬁles, sea level anomaly and sea surface temperature satellite data). An operational quality assessment framework is used to evaluate the accuracy of the products which set the basic standards for the future upgrades, highlighting the strengths and weaknesses of the model and the observing system in the Black Sea.


Introduction
Operational forecasting is now a reality in most ocean areas around the world. It follows the example of weather forecasting [1] and provides products that are crucial for the sustainable development of activities at sea and along the coasts.
The operational quality of the products has steadily improved since the late 1990s, when only a few centers around the world were engaged in ocean analysis and forecasting [2]. The operational system assessment allows Earth System Science to estimate the quality and fitness of the numerical ocean model for purpose of the observing system.
Operational forecasting in the Black Sea region is part of the data production architecture of the Copernicus Marine Environment and Monitoring Service (CMEMS, https://marine.copernicus.eu/ accessed on 29 December 2021, [3]). The Black Sea Monitoring and Forecasting Center (BS-MFC) has been operative since 2016. The center provides capacity in the Southern Regional European Seas [6]. These operational products also reconstruct the past state of the Black Sea, thus providing the optimal data set to study the ocean climate variability and in general the Black Sea general circulation. We are currently at the third generation of forecasting systems, and this paper analyzes the accuracy of the first ocean forecasting system for the Black Sea physics.
This study presents the numerical setup, operational implementation and product quality assessment for the period January 2018 to December 2020, using observations from in-situ (temperature and salinity profiles) and satellite (sea surface temperature and sea level anomaly) platforms, provided by CMEMS Thematic Assembly Centers (TAC) [3].
The paper is organized as follows: Section 2 presents the system, including the ocean model, the data assimilation method used for the forecasting cycle and the processing chain; Section 3 describes the operational products and discusses the product quality; Section 4 presents the conclusions and future evolutions of the BSFS.

Ocean Numerical Model
BSFS is based on a free surface implementation of the NEMO hydrodynamical model (v3.4, [7]) for the Black Sea region. The horizontal resolution is approximately 3 km (1/27° in zonal and 1/36° in meridional directions), which conforms to the mesoscale eddyresolving scale (Rossby radius of deformation in the Black Sea is ~20 km, [8]). On the other hand, 31 unevenly distributed z-levels are used along the vertical direction, with an initial layer at about 2.5 m in depth. BSFS bathymetry is derived from GEBCO at 1′ resolution [9]; a local refinement of the coastline, using a high-resolution NOAA dataset [10], is performed to represent the main coastal features of the basin.
The model is forced by momentum, water and heat fluxes interactively computed by bulk formulae implemented first for the Mediterranean Sea [11] and adapted for the Black Sea, as described in Appendix A. The analysis and forecast atmospheric fields are The system is based on three major components: collection of upstream data, including atmospheric forcing and observations, the numerical ocean model, and a variational data assimilation scheme. Together with the Mediterranean Sea forecasting system, which began operation in 1998 [4,5], the Black Sea one completed the operational capacity in the Southern Regional European Seas [6]. These operational products also reconstruct the past state of the Black Sea, thus providing the optimal data set to study the ocean climate variability and in general the Black Sea general circulation. We are currently at the third generation of forecasting systems, and this paper analyzes the accuracy of the first ocean forecasting system for the Black Sea physics.
This study presents the numerical setup, operational implementation and product quality assessment for the period January 2018 to December 2020, using observations from in-situ (temperature and salinity profiles) and satellite (sea surface temperature and sea level anomaly) platforms, provided by CMEMS Thematic Assembly Centers (TAC) [3].
The paper is organized as follows: Section 2 presents the system, including the ocean model, the data assimilation method used for the forecasting cycle and the processing chain; Section 3 describes the operational products and discusses the product quality; Section 4 presents the conclusions and future evolutions of the BSFS.

Ocean Numerical Model
BSFS is based on a free surface implementation of the NEMO hydrodynamical model (v3.4, [7]) for the Black Sea region. The horizontal resolution is approximately 3 km (1/27 • in zonal and 1/36 • in meridional directions), which conforms to the mesoscale eddy-resolving scale (Rossby radius of deformation in the Black Sea is~20 km, [8]). On the other hand, 31 unevenly distributed z-levels are used along the vertical direction, with an initial layer at about 2.5 m in depth. BSFS bathymetry is derived from GEBCO at 1 resolution [9]; a local refinement of the coastline, using a high-resolution NOAA dataset [10], is performed to represent the main coastal features of the basin.
The model is forced by momentum, water and heat fluxes interactively computed by bulk formulae implemented first for the Mediterranean Sea [11] and adapted for the Black Sea, as described in Appendix A. The analysis and forecast atmospheric fields are provided by the European Centre for Medium Range Weather Forecast (ECMWF) at 0.125 • horizontal resolution, at 3-6 h frequency.
The momentum equation is written in flux form and is solved with a leapfrog time stepping scheme. The free surface equation uses a linearized form (i.e., the barotropic velocity field is defined with an integral from the flat surface to the bottom) and is integrated implicitly [12] with the time step of the total velocity field equal to 150 s. The advection scheme for the temperature and salinity tracers is the total variance dissipation (TVD) scheme. Horizontal eddy diffusivity is applied for tracers using a Laplacian operator, with a coefficient of 50 m 2 /s, while a bi-Laplacian viscosity is applied for momentum, with a coefficient of 1.0 × 10 8 m 4 /s. The vertical diffusivity and viscosity are computed from a turbulent kinetic energy (TKE) closure model, with the parameters set as in [7]. Vertical eddy viscosity and diffusivity coefficients are selected with values of 1.2 × 10 −5 m 2 /s and 1.2 × 10 −6 m 2 /s, respectively. The model also implements free-slip lateral boundary conditions and a classical quadratic bottom friction scheme, with a drag coefficient of The Bosporus Strait connects the Black Sea with the Sea of Marmara: to model it, a twolayer, narrow strait water mass exchange which maintains a relatively steady state water and salt balance in the Black Sea [13][14][15][16][17][18] is used, with a closed boundary condition. The excess precipitation and river runoff over evaporation is removed using the outflow from the strait, thus leading to a zero balance in the Black Sea [19]. Considering the horizontal average of the free surface equation [20], the net transport at the Bosporus Strait is given by: where Tr B is the net transport at the Bosporus Strait, and A is the basin surface area, η is the sea surface height, E, P, R are, respectively, evaporation, precipitation and runoff, the Dirac δ is different from zero at 72 river mouths and the triangular brackets mean horizonal basin average. First, we redefine the transport at the Bosporus Strait in terms of a "discharge" R B : Thus, the discharge at the Bosporus Strait is given by: where H i,B and L i,B are depth and width in one-grid-sea-point, respectively. We calculated the Bosporus discharge in (3) from a 10-year simulation by computing the mean free surface tendency and the mean water flux. The values of the Bosporus discharge are stored as monthly mean values and set as vertical velocity boundary conditions as done for the rivers, except with the negative sign, indicating a discharge out of the basin or a "negative river". This parametrization is quite robust for decadal long simulations that do not consider climate change trends in sea level and water fluxes.
With regard to the 72 real river runoff contributions, we use monthly mean discharge data from the SESAME dataset [21] for all rivers, including the Danube, the Dniepr, the Dniester, the Rioni, the Kizil Irmak, and the Sakarya.

Data Assimilation Scheme
The ocean model is coupled with a data assimilation system in order to produce analyses for optimal initial conditions of the forecasts. The DA ocean state vector x contains the temperature T, salinity S and sea level values at each model grid location. For a model setup with n vertical levels, this vector is defined as: x = (T 0 , . . . , T n , S 0 , . . . , S n , η) Given an observation vector y 0 and the background model state x b it is possible to define the innovation, i.e., the difference between the observations and their respective model predictions, denoted d: with H(x b ) being the observation operator that projects the model state onto the observation space. The aim of the DA system is to find a correction to the model state: This correction needs to minimize the analysis error while taking into account the error covariances of the model state vector and the observations. There are several methods to calculate such corrections; however, they can be broadly categorized as variational methods or Kalman filters. In BSFS, a 3D-variational method is used [22,23]. In this method, the corrections are derived by iteratively minimizing a cost function J, defined as: Here the matrices B and R are the error covariance matrices of, respectively, the background state and the observations. For a number of state variables n b , the model covariance matrix B is of size n 2 b . In addition, as the cost function contains B −1 , the covariance matrix needs to be inverted. It is clear that for large numbers of state variables, this is a computationally costly calculation, which can be avoided if additional constraints are imposed on B, in particular if B is of the form: with V an arbitrary matrix. In this case, Equation (7) can be written as: introducing a coordinate transformation: The minimization of the cost function can now be performed in terms of the new control variable v, without the need to calculate B −1 . It is sufficient to perform the minimization in terms of v before transforming back to the model state space increment δx with (10).
As the covariance is by definition a positive definite matrix, the matrix V exists and could be found by performing a Cholesky decomposition on B. However, since the model error covariance B is usually derived from model state variable anomalies in a long model run, V is simply the anomaly matrix and B is already calculated according to (8).
One of the main features of OceanVar is the covariance decomposition. The matrix V of (8) is expressed as a product of different components: where V H is the horizontal and V V is the vertical component of the covariance. The horizontal component provides the correlation between neighboring grid points at each model level. This is implemented by means of a recursive filter with a radius of approximately 25 km in open sea and a gradual falloff near the coast. The vertical component is calculated from a long model simulation and decomposed using empirical orthogonal functions (EOF). This EOF decomposition significantly reduces the number of variables, n v n x , and re-duces spurious correlations due to the finite dataset used to estimate the covariance. V η is the dynamic height operator, which uses the local hydrostatic adjustment scheme [23] to transform the sea level anomaly innovations into increments for T and S in regions deeper than 1000 m. For the BSFS system with 31 vertical levels, the vertical covariance of T and S is represented by 15 EOF. The EOF are calculated separately for each month and each grid location to adequately capture the variability of the covariance in time and in space. The DA system is used to assimilate in-situ observations of temperature and salinity from available ARGO profilers in the Black Sea, and satellite observations of sea level anomaly (SLA) and of sea surface temperature (SST). For the latter, the foundation temperature is used and assimilated only at nighttime (as in e.g., [24]). This approach limits possible biases due to the diurnal cycle of SST.

Operational Chain
BSFS version here considered is EAS3 (European Analysis and forecasting System, version 3, referred also in the next sections) for the Black Sea. It implements technical interfaces with upstream data, used for model runs-including ECMWF atmospheric forcing, CMEMS in-situ and satellite observations from corresponding TACs-as well as with the CMEMS Dissemination Unit for the delivery of analysis and forecast products to users. Observations used by BSFS for assimilation and validation are summarized in Table 1. The atmospheric forcing (previous day analyses and 10 days forecast starting at 00:00 UTC) is downloaded from ECMWF through Aeronautica Militare Italiana. As soon as it is available, the BSFS processing system starts (at around 7:00 UTC). The atmospheric forcing availability is a major source of delay for the oceanic forecast.
Every day the BSFS (Figure 1) produces ten days of forecast (J to J+9), one day of simulation (J−), and three days of analysis (J−4 to J−1). Once a week (on Tuesday), 14 days of analysis are produced, from J−15 to J−1, with the assimilation of all available satellite and in-situ data collected over the past weeks. Analysis produced at the weekly cycle represents the best estimation of the Black Sea state because all in-situ and satellite altimetry high quality processing data is used up to J−3. Analysis and simulation runs are forced by ECMWF atmospheric analysis fields at 6 h frequency; the forecast cycle is forced with ECMWF atmospheric forecast fields at 3 h frequency for the first 3 days and 6 h frequency for the remaining 7 days.
At the end of the forecasting cycle, post-processing involves the preparation of all the BSFS files in a format that is compliant with the CMEMS and CF standards, delivering to the CMEMS DU, and archiving of the BSFS native products at CMCC's supercomputing facilities. Every day, the forecast is released to CMEMS within less than 3 h (target delivery time: 12 UTC).
The BSFS product catalogue includes daily and hourly mean datasets, centered at midday of each J, for the Black Sea essential variables: 3D temperature (T), salinity (S), zonal and meridional velocity components (U and V, respectively) and 2D sea surface height (SSH), bottom temperature (BottomT) and mixed layer depth (MLD).
The BSFS operational chain is very similar to that of the Mediterranean Forecasting System (MFS, [25]) in terms of general setup (technical interface, analysis-simulationforecast chain). Unlike MFS, it implements a more frequent analysis cycle to provide higher quality Black Sea products using daily observations.

Quality Assessment of the Operational System
The aim of the product quality assessment is to monitor the analysis quality and forecast accuracy of BSFS products using quasi-independent validation assessment. Operationally, a regional website provides daily bulletins (http://bsfs.cmcc.it/, accessed on 29 December 2021) and skill scores (http://oceanlab.cmcc.it/bsfs-evaluation/, accessed on 29 December 2021). Daily bulletins consist of a collection of interactive 2D maps for visualizing the essential variables of the Black Sea (as provided in the CMEMS catalogue). Weekly skill scores are provided through the evaluation section: during the assimilation weekly cycle, as described in Section 2.3, the difference between the model analysis and the observations at the time and location of the observations (i.e., misfits) is stored to compute statistics. Two main statistical metrics for the assessment of the analysis quality are computed using the misfits: the first is the bias, given by: and the Root Mean Square Error (RMSE) given by: where N is the number of data used in the evaluation, φ M is the model analysis field, and φ O the observation.
On the other hand, to assess the quality of the forecast, we evaluated the differences of the forecast fields with respect to the analysis considered to be the best estimate of the truth. Murphy (1993, [26]), revisited by [2], defined the "forecast goodness" methodology and here we use his "quality measure" methodology that compares the forecast with the analysis and the "persistence" or the best estimate at initial time. If the difference between forecast and analysis is better than the difference between forecast and persistence, then the forecast is valuable. For a modern oceanographic forecasting system, [27] defined a metrics that we will partially follow here by computing:

•
The difference between the analysis and the forecast (AF): where φ AN and φ FC is the temperature (salinity) daily mean forecast and analysis, respectively, at each forecast day t = day1, day2, . . . , day10 (ID day as in Figure 2); T corresponds to the time period covered by the evaluation (here we consider a time period of 1 month). The metrics is normalized horizontally by averaging over the area and at a specific selected depth.
lated initial condition at day J−2.

Sea Surface Temperature
The BSFS analysis sea surface temperature is assessed by comparing analysis model fields against SST satellite data remapped over the Black Sea basin at 1/16° spatial resolution and representative of night SST values and delivered by CMEMS (Table 1). Figure 3 shows the time series from 2018 to 2019 of the difference between BSFS and satellite SST (BIAS, top) and RMSE (bottom). The numerical SST is slightly warmer than the measured one (+0.1 °C), and the error is about 0.5 °C.

•
The difference between the forecast and the persistence (PF): Here, φ AN (t = day0) is the persistence (ID day as in Figure 2), e.g., the initial condition for the forecast. The latter is the last time step of a simulation started from an assimilated initial condition at day J−2.
In the following subsections, we present BSFS analysis statistics for the period 2018-2019 for SST and SLA, as well as temperature and salinity at given layers

Sea Surface Temperature
The BSFS analysis sea surface temperature is assessed by comparing analysis model fields against SST satellite data remapped over the Black Sea basin at 1/16 • spatial resolution and representative of night SST values and delivered by CMEMS (Table 1). Figure 3 shows the time series from 2018 to 2019 of the difference between BSFS and satellite SST (BIAS, top) and RMSE (bottom). The numerical SST is slightly warmer than the measured one (+0.1 • C), and the error is about 0.5 • C.

Analysis Quality
Water column properties given by BSFS are assessed after 3D temperat ity are validated against all available observations. Table 2 shows the basin scale RMSE misfits at specific layers (m) for the 2-5 m surface layer, the average RMSE for temperature at the whole bas °C ( Figure 5(a1)) and 0.2 PSU (Figure 5(b1)). However, the error in temper increase in the sub-surface, from 5 m to 30 m, from 0.5 °C to 2 °C in the summ 0.5 °C in the winter ( Figure 5(a2-a4)).
The maximum error for temperature occurs in the 10-20 m layer with 2.5 °C computed for August-October in both 2018 and 2019 ( Figure 5(a3)).
The increased error in the thermocline (10-30 m layer) and in the halo m layer) in the summer is likely due to a deficiency in the vertical discretiz only 31 z-levels) and in the vertical mixing parameterization.
The errors are due to (a) the Bosporus Strait being represented as clos (b) mixing processes that were not completely resolved since effects indu and tides are not accounted for, and (c) simplified representation of the hea and water fluxes through bulk formulation, whose effects are particularly pecially during the summer. These are significant challenges to be conside velopment of the next generation of the Black Sea physical system.

Analysis Quality
Water column properties given by BSFS are assessed after 3D temperature and salinity are validated against all available observations. Table 2 shows the basin scale RMSE misfits at specific layers (m) for 2018-2019. In the 2-5 m surface layer, the average RMSE for temperature at the whole basin is about 0.4 • C ( Figure 5(a1)) and 0.2 PSU (Figure 5(b1)). However, the error in temperature starts to increase in the sub-surface, from 5 m to 30 m, from 0.5 • C to 2 • C in the summer and below 0.5 • C in the winter ( Figure 5(a2-a4)). The maximum error for temperature occurs in the 10-20 m layer with a value up to 2.5 • C computed for August-October in both 2018 and 2019 ( Figure 5(a3)).
Regarding salinity, the error in the subsurface, on average, is around 0.2 PSU ( Figure 5(b2-b4)), which increases from 30 m up to 100 m up to a maximum value of 0.4 PSU ( Figure 5(b5-b7)). In the intermediate layers up to 1000 m, the error decreases to quasi-zero for both temperature ( Figure 5(a8)) and salinity ( Figure 5(b8)).
The increased error in the thermocline (10-30 m layer) and in the halocline (50-100 m layer) in the summer is likely due to a deficiency in the vertical discretization (we use only 31 z-levels) and in the vertical mixing parameterization.   The errors are due to (a) the Bosporus Strait being represented as closed boundary, (b) mixing processes that were not completely resolved since effects induced by waves and tides are not accounted for, and (c) simplified representation of the heat, momentum and water fluxes through bulk formulation, whose effects are particularly important especially during the summer. These are significant challenges to be considered in the development of the next generation of the Black Sea physical system.

Forecast Assessment
In Figures 6 and 7 we present the BSFS forecast quality assessment for temperature and salinity considering 2 reference months, February 2020 and August 2020, respectively. Metrics are given at different depths-2.5, 30 and 150 m. If the forecast is skillful, AF score should be better than PF score. and salinity considering 2 reference months, February 2020 and August 2020, respectively. Metrics are given at different depths-2.5, 30 and 150 m. If the forecast is skillful, score should be better than score. In February 2020, is bigger than after 2 (Figures 6 and 7): this means that after the first day the advantage of doing the forecast is clear over persisting the initial condition. It is interesting to note that for salinity this advantage grows more rapidly, i.e., for salinity the forecast overcomes the persistence even before the first day. Over the water column, temperature and salinity analysis and forecast exhibit a similar variability: higher error reaching a maximum around 0.6° Celsius and 0.5 PSU at 10 is shown for temperature ( Figure 6(a1,a2)) and salinity ( Figure 6(b1,b2)) in the subsurface up to 30 m; at 150 m, both variables exhibit a lower error (never higher than 0.1° Celsius and 0.3 PSU over the reference 10-day period). In August 2020, errors at 30 m are higher than the ones computed in the previous analyzed month: starting from 2, in particular for temperature (Figure 7(a1,a2)), shows higher values than -almost 50% higher at 30 m, with maximum error of about 2.0° Celsius (Figure 7(a2)). This means that during the summer the forecast adds to persistence more than during winter: the dynamical model is clearly responsible for the proper changes in water column stratification even for few days. At 30 m and in August (Figure 7(a1)) the difference between and is the largest: we assume that this is an indication of the complex dynamics involved in the mixing and the CIL dynamics. Regarding salinity (Figure 7(b1-b3)

Temporal Evolution of the CIL and Stratification
The BSFS is able to represent the formation and evolution of the cold intermediate layer (CIL), which is a typical water mass property of the Black Sea basin, and well documented in the literature ( [16] for a very comprehensive literature review and updated evaluation). It is identified by the 8.0 °C isotherm, arising at the surface during winter, penetrating the subsurface (typically at a depth range of 50-100 m), intruding the warmer zone, and persisting over time. The time versus depth diagrams at the observation location are shown for temperature and for salinity in Figure 8, between 2015-2019. Figure 8a shows the cold event of 2017, which is also captured by observations as demonstrated by the quite low difference between the model and observations (Figure 8b). It also shows the reduction in the CIL width because of the warming period from 2018 and its perforation and partial disappearance in 2019 [16]. The major differences of about 2 °C occur in the upper layer and are shown in Figure 8b. Salinity stratification is well represented in Figure 9a with small differences with respect to observations, as shown in Figure 9b  In February 2020, PF is bigger than AF after day2 (Figures 6 and 7): this means that after the first day the advantage of doing the forecast is clear over persisting the initial condition. It is interesting to note that for salinity this advantage grows more rapidly, i.e., for salinity the forecast overcomes the persistence even before the first day. Over the water column, temperature and salinity analysis and forecast exhibit a similar variability: higher error reaching a maximum around 0.6 • Celsius and 0.5 PSU at day10 is shown for temperature ( Figure 6(a1,a2)) and salinity ( Figure 6(b1,b2)) in the subsurface up to 30 m; at 150 m, both variables exhibit a lower error (never higher than 0.1 • Celsius and 0.3 PSU over the reference 10-day period).
In August 2020, errors at 30 m are higher than the ones computed in the previous analyzed month: starting from day2, in particular for temperature (Figure 7(a1,a2)), PF shows higher values than AF-almost 50% higher at 30 m, with maximum error of about 2.0 • Celsius (Figure 7(a2)). This means that during the summer the forecast adds to persistence more than during winter: the dynamical model is clearly responsible for the proper changes in water column stratification even for few days. At 30 m and in August (Figure 7(a1)) the difference between PF and AF is the largest: we assume that this is an indication of the complex dynamics involved in the mixing and the CIL dynamics. Regarding salinity (Figure 7(b1-b3)), PF and AF errors are quite similar comparing to winter and of the order of 0.3-0.4 PSU.

Temporal Evolution of the CIL and Stratification
The BSFS is able to represent the formation and evolution of the cold intermediate layer (CIL), which is a typical water mass property of the Black Sea basin, and well documented in the literature ( [16] for a very comprehensive literature review and updated evaluation). It is identified by the 8.0 • C isotherm, arising at the surface during winter, penetrating the subsurface (typically at a depth range of 50-100 m), intruding the warmer zone, and persisting over time. The time versus depth diagrams at the observation location are shown for temperature and for salinity in Figure 8, between 2015-2019. Figure 8a shows the cold event of 2017, which is also captured by observations as demonstrated by the quite low difference between the model and observations (Figure 8b). It also shows the reduction in the CIL width because of the warming period from 2018 and its perforation and partial disappearance in 2019 [16]. The major differences of about 2 • C occur in the upper layer and are shown in Figure 8b. Salinity stratification is well represented in Figure 9a with small differences with respect to observations, as shown in Figure 9b

Temporal Evolution of the CIL and Stratification
The BSFS is able to represent the formation and evolution of the cold intermediate layer (CIL), which is a typical water mass property of the Black Sea basin, and well documented in the literature ( [16] for a very comprehensive literature review and updated evaluation). It is identified by the 8.0 °C isotherm, arising at the surface during winter, penetrating the subsurface (typically at a depth range of 50-100 m), intruding the warmer zone, and persisting over time. The time versus depth diagrams at the observation location are shown for temperature and for salinity in Figure 8, between 2015-2019. Figure 8a shows the cold event of 2017, which is also captured by observations as demonstrated by the quite low difference between the model and observations (Figure 8b). It also shows the reduction in the CIL width because of the warming period from 2018 and its perforation and partial disappearance in 2019 [16]. The major differences of about 2 °C occur in the upper layer and are shown in Figure 8b. Salinity stratification is well represented in Figure 9a with small differences with respect to observations, as shown in Figure 9b

Circulation
The surface mean circulation is shown in Figure 10 as monthly means for 2019. Considering the poor availability of observations, a robust validation exercise is extremely difficult, thus the main features emerging from the annual monthly means are described

Circulation
The surface mean circulation is shown in Figure 10 as monthly means for 2019. Considering the poor availability of observations, a robust validation exercise is extremely difficult, thus the main features emerging from the annual monthly means are described in detail with respect to those extensively described in [28][29][30][31].  The Black Sea surface dynamics is characterized by a main cyclonic gyre, the Rim current, encompassing the basin and a variety of mesoscale eddies along the coast, some of them quasi-stationary. BSFS captures most of the particular dynamical structures in the basin, such as the Rim current, which persists over 2019, and small-scale structures such as coastal anticyclonic eddies, which appear along the Russian-Georgian coastline (one quite intense in February and progressively weaker in March) and coastal cyclonic eddies, which are much more intermittent over the year and weak). The Batumi Gyre, located near the Georgian coast, forms in April and lasts until June, when the circulation in that region is fragmented in more unstable and weak eddies until its regeneration in August-September. Progressively towards the southwestern region, small coastal eddies become part of the Rim current. On the Turkish coastline, cyclonic eddies like the one in the area of Trabzon (January-February; November-December) coexist with anticyclonic eddies, as in the Synop area (June-October). In the area between Istanbul and Burgas-the Bosporus region-an anticyclonic eddy appears from October to December. Close to the Danube, coastal currents are quite well captured by the model, then going further towards the continental slope (30 • E-32 • E and 44 • N-45 • N) the Sevastopol anticyclonic gyre starts in May (very weak), becoming stronger in November and December.

Conclusions
Within the framework of the CMEMS and the BS-MFC, the Black Sea physical analysis and forecasting system provides essential variables for understanding the physical processes and dynamics of the Black Sea basin. The BSFS has been operational since the end of 2016 and has been developed and maintained at CMCC in collaboration with the USOF (University of Sofia, Bulgaria, scientific partner in the BS-MFC consortium).
We have presented the BSFS ocean model based on NEMO v3.4 able to assimilate near real time, in-situ and satellite observational products using the OceanVar scheme. The Black Sea hydrodynamic model has about 3 km horizontal resolution and uses 31 levels with partial steps. It implements a closed boundary condition at the Bosporus Strait and is forced by ECMWF analysis and forecasting atmospheric fields.
The BSFS implements two production cycles, one daily (which includes 3 days analysis) and one weekly (based on 14 days analysis), the latter to assimilate a higher number of collected observations to provide the best quality initial conditions for the forecasting cycle: the processing system is completed every day by 1-day simulation and 10-day forecasts.
We have also described an operational dashboard for product quality monitoring which assesses the skill scores of the analysis products for 2018-2019 and gives an accuracy of around 1 • C in the sea surface temperature. Errors increase at the sub-surface (10-30 m layer): in particular, the thermocline experiences a maximum error of 2 • C during the summer period, while salinity reaches an error of about 0.5 PSU. Considering the sea level, the assimilation of the along-track satellite SLA guarantees an average error of about 2.3 cm. Such skills put the BSFS on the same track of quality and robustness as state-of-the-art regional configurations in the CMEMS framework. A regional website is operationally maintained at CMCC to provide daily bulletins and metrics for monitoring the lifecycle and performances of the system. The future forecasting system for the Black Sea will include at least four new main components which will significantly improve the quality of the BSFS analysis fields and forecasting skill scores by means of improved hydrodynamical core model: (1) increased vertical resolution for a better representation of open ocean dynamics and mixing processes, combined with a revised bathymetry and coastline and data assimilation upgrades; (2) the Bosporus Strait will work as an open boundary, in order to improve the connection with the Mediterranean Sea. This will be achieved thanks to a novel implementation of the Marmara Sea model, based on an unstructured grid method that provides open boundary conditions to the Black Sea through the Bosporus and the Mediterranean Sea through the Dardanelles; (3) improved representation of the land forcing: in particular, the representation of the Danube River using historical discharge datasets provided by the NIHWM (National Institute of Hydrology and Water Management, Romania-scientific partner in the BS-MFC consortium); (4) an online wave-current model to improve the small scales dynamics at the surface.
Author Contributions: S.A.C. coordinated the work in collaboration with G.C. and N.P. The ocean circulation model has been designed and improved by S.A.C., E.P., N.P., G.C., M.I.; the data assimilation scheme has been designed by E.J. with the collaboration of L.L.; R.L. contributed on upstream data access; Funding acquisition, A.P.; S.C. (Sergio Creti') and L.S. contributed on operational chain design and operational maintenance of the system. Validation has been developed in collaboration with E.J., D.A. and S.C. (Salvatore Causio). All authors have read and agreed to the published version of the manuscript. where ρ A is the density of the moist air, C P is the specific heat capacity, C h is the turbulent exchange coefficient for humidity set as a constant and equal to 1.3·10 −3 .
The Q e is computed as follows: where L e is the latent heat of vaporization, C e 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 q A is the specific humidity of air. The two constants, C h and C e , 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: where V x , V y are the wind components, while C D is the drag coefficient given as a function of wind speed and temperature difference T A − 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: