The Mediterranean Forecasting System – Part 1: Evolution and performance

The Mediterranean Forecasting System produces operational analyses and reanalyses and 10 d forecasts for many essential ocean variables (EOVs), from currents, temperature, salinity, and sea level to wind waves and pelagic biogeochemistry. The products are available at a horizontal resolution of $\mathrm{1}/\mathrm{24}$^∘ (approximately 4 km) and with 141 unevenly spaced vertical levels.

The core of the Mediterranean Forecasting System is constituted by the physical (PHY), the biogeochemical (BIO), and the wave (WAV) components, consisting of both numerical models and data assimilation modules. The three components together constitute the so-called Mediterranean Monitoring and Forecasting Center (Med-MFC) of the Copernicus Marine Service.

Daily 10 d forecasts and analyses are produced by the PHY, BIO, and WAV operational systems, while reanalyses are produced every ∼ 3 years for the past 30 years and are extended (yearly). The modelling systems, their coupling strategy, and their evolutions are illustrated in detail. For the first time, the quality of the products is documented in terms of skill metrics evaluated over a common 3-year period (2018–2020), giving the first complete assessment of uncertainties for all the Mediterranean environmental variable analyses.

1 Introduction

Ocean analysis and forecasting systems are now available for the global ocean and regional seas at different spatial scales and with different numbers of essential ocean variables (EOVs) considered (Tonani et al., 2015). The societal drivers for the operational products that stem from the ocean analysis and forecasting products are the safety of maritime transport, multiple coastal hazards, and climate anomalies. Moreover, the operational products are at the basis of a new understanding of the dynamics of the ocean circulation (Pinardi et al., 2015); its linked biogeochemical cycles, namely carbon uptake and eutrophication, among others (Canu et al., 2015; von Schuckmann et al., 2020); and extreme storm surge events (Giesen et al., 2021).

The ocean analysis and forecasting system for the entire Mediterranean Sea was set up over the past 15 years (Pinardi and Coppini, 2010; Pinardi et al., 2017; Lazzari et al., 2010; Salon et al., 2019; Ravdas et al., 2018; Katsafados et al., 2016), and in 2015, it became operational in the framework of the Copernicus Marine Service, which is the marine component of the Copernicus Programme European Union service for sustainable use of the ocean, providing free, regular, and systematic information on the state of the blue (physical), white (sea ice), and green (biogeochemical) ocean on the global and regional scales. The Copernicus Marine Service in Europe has shown the strength of a state-of-the-art operational service implemented by hundreds of experts and teams who are distributed throughout Europe and come from public and private sectors, from operational and research organizations, from different countries, and from diverse cultures and relations to the ocean (Le Traon et al., 2017; Alvarez Fanjul et al., 2022). In this paper, we give an overview of the core components of the system, i.e. the numerical models and the data assimilation modules that represent the eddy-resolving ocean general circulation, the biogeochemical tracers, and the wind waves. Furthermore, we will document the quality of EOV products using goodness indices (Brassington, 2017). The core components constitute the so-called Mediterranean Monitoring and Forecasting Center (Med-MFC) of the Copernicus Marine Service (Le Traon et al., 2019). The integrated approach of the Med-MFC system represents a unique opportunity for the users to access state-of-the-art information provided in a uniform manner (e.g. same grid, unique format, and unique point of access). This ocean analysis and forecasting system, hereafter Med-MFC, produces analyses, 10 d forecasts, and reanalyses (Adani et al., 2011; Pinardi et al., 2015; von Schuckmann et al., 2016, 2018, 2019; Terzic et al., 2021; Simoncelli et al., 2016, 2019; Ravdas et al., 2018; Escudier et al., 2020, 2021; Cossarini et al., 2021).

An essential task of the production activities concerns the continuous assessment of the quality of the products (Sotillo et al., 2021; Alvarez Fanjul et al., 2022), which is achieved at two levels: (i) the pre-qualification of the systems before delivering a new release, including an extensive scientific validation of the products, published in the QUality Information Documents (QUIDs) available on the Copernicus Marine Product Catalogue; (ii) the operational evaluation of the skill metrics during operations, made available through the Copernicus Marine Product Quality Dashboard website (https://pqd.mercator-ocean.fr, last access: 30 June 2023), as well as through the Mediterranean regional validation websites implemented at the level of the Med-MFC production units (PHY: https://medfs.cmcc.it/, last access: 30 June 2023; WAV: http://Med-MFC-wav.hcmr.gr/, last access: 30 June 2023; BIO: https://medeaf.ogs.it/nrt-validation, last access: 30 June 2023). All the delivered variables are thus validated with respect to satellite and in situ observations using Copernicus Marine observational datasets, as well as additional datasets, climatologies, or literature information when needed.

The Mediterranean Sea is a semi-enclosed basin with an anti-estuarine circulation corresponding to a 0.9 and 0.8 ± 0.06 Sv baroclinic inflow and outflow at the Strait of Gibraltar, positive energy inputs by the winds, and net buoyancy losses inducing a vigorous overturning circulation (Cessi et al., 2014; Pinardi et al., 2019). The basin-scale circulation is dominated by mesoscale and sub-mesoscale variability (Pinardi et al., 2015; Bergamasco and Malanotte-Rizzoli, 2010; Pinardi et al., 2006; Robinson et al., 2001; Ayoub et al., 1998), with the former being subdivided into semi-permanent and synoptic mesoscales with spatial scales larger than 4–6 times the local Rossby radius of deformation. The stratification is large during summer in the first 50 m, and during winter the water column is practically unstratified. The Mediterranean Sea is an oligotrophic basin (Siokou-Frangou et al., 2010) with a west-to-east decreasing productivity gradient (Lazzari et al., 2012) and relatively high primary productivity in open-ocean areas where winter mixing increases surface nutrients (Cossarini et al., 2019). The wave conditions are driven by the winter storminess, while summer is characterized by low significant wave height values and higher value scatter (Ravdas et al., 2018). The yearly mean wave period, as estimated from available wave buoys over the Mediterranean Sea, amounts to 3.82 s with typical deviations of 0.92 s, while the mean significant wave height is 0.82 m (1.28 m as estimated by satellite observations) with typical deviations of 0.67 m (0.77 m for satellite data).

In this paper, we describe the final setup of the Med-MFC core components for the period 2017–2020. The Med-MFC modelling systems share the same grid resolution ($\mathrm{1}/\mathrm{24}$^∘) and bathymetry and use the same atmospheric and river forcing fields. Moreover, daily mean fields evaluated by the physical model are used to force the wave component (surface currents) and the transport–biogeochemical model (temperature, salinity, horizontal and vertical velocities, sea level, and diffusivity). This allows several model parameterizations to be calibrated to obtain the best result in terms of the specific environmental variable considered by each component. In the Copernicus Marine Service, the approach of forcing waves and biogeochemistry models with information from the hydrodynamic models is used and represents a standard which is also applied for the other MFCs. Several MFCs also use the online coupling between physics and wave models and between physics and biogeochemical models. Furthermore, this weakly coupled system ensures an efficient development of the data assimilation modules connected to each numerical model module and specific input datasets. It is a distributed system that shares information with efficiency and effectiveness when and how it is required by relevant processes. Due to its rather unique structure and the quality of its products, the system described could be used as a basic standard for new systems to be developed.

The paper is organized as follows. Section 2 provides an overview of the technical specifications of the Med-MFC components, Sect. 3 describes the quality of the system for a reference period from 2018 to 2020 and the quality of the forcing, and Sect. 4 concludes the paper and presents future perspectives.

2 Description of the Med-MFC core components

The structure of the Med-MFC core components is shown in Fig. 1: the physical component (PHY) is composed of the NEMO general circulation model (Madec and the NEMO system Team, 2019) coupled to the WaveWatch-III (WW3) wave model (Clementi et al., 2017a) and the ocean data assimilation OceanVar 3DVAR (Dobricic and Pinardi, 2008 and Storto et al., 2015); the biogeochemical component (BIO) is composed of the Biogeochemical Flux Model (BFM), the tracer transport OGS Transport Model (OGSTM), and a data assimilation scheme (Lazzari et al., 2012, 2016; Cossarini et al., 2015; Vichi et al., 2020), forced daily by the daily mean of the PHY component fields; the wave component (WAV) is composed of the wave model WAM (WAMDI Group, 1988) and its assimilation scheme, forced daily by the daily mean of the PHY component fields. Daily 10 d forecasts, as well as analyses and reanalyses, are produced with all PHY, BIO, and WAVE components, as described below.

Figure 1 The Med-MFC core components and the offline coupling scheme. The blue arrows are the exchanged fields at daily frequency between the three components.

Each component of the Med-MFC has its own data assimilation system; as such, an important effort was made to extract the most relevant information from satellite and in situ observations to produce the analysis and correct initial conditions for the forecast in order to benefit the forecasting skills. The main goal of the paper is to present the current quality of the operational system components by comparing the analysis and – for specific variables, such as significant wave height – the background (simulation) with in situ and/or satellite observations. The skill of the wave and biogeochemical models is assessed by considering intercomparisons of the model solution during the 24 h analysis phase with in situ and remotely sensed observations. As the latter are ingested into the model through data assimilation, the first-guess model fields (i.e. model background) are used instead of analyses.

2.1 The general circulation model component

2.1.1 Numerical model description

The PHY numerical model component comprises a two-way coupled current–wave model based on NEMO and WW3 implemented over the whole Mediterranean basin and extended into the Atlantic Sea in order to better resolve the exchanges with the Atlantic Ocean (Fig. 2). The model horizontal grid resolution is $\mathrm{1}/\mathrm{24}$^∘ (ca. 4 km) and is resolved along 141 unevenly spaced vertical levels (Clementi et al., 2017b, 2019). The topography is an interpolation of the GEBCO (General Bathymetric Chart of the Oceans) 30 arcsec grid (Weatherall et al., 2015) that is filtered and specifically modified in critical areas such as the eastern Adriatic coastal areas (to avoid instabilities in circulation due to the presence of a large number of small islands), the Gibraltar and Messina straits (to better represent the transports), and the Atlantic edges external border (to avoid large bathymetric inconsistencies with respect to the Copernicus Global Analysis and Forecast product in which the model is nested). All the numerical model choices are documented in Table A1 in the Appendix.

Figure 2 The solid red box presents the domain of the PHY and WAV Mediterranean components. For BIO, the domain extends into the Atlantic as far as the dashed red line. The blue box presents one of the WAM domains, producing boundary conditions for the Mediterranean WAV component which extends only into the solid red box.

The general circulation model considers the nonlinear free-surface formulation and vertical z-star coordinates. The numerical scheme uses the time-splitting formulation to solve the free-surface and the barotropic equations with a (100 times) smaller time step with respect to the one used to evaluate the prognostic 3D variables (240 s). The active-tracer (temperature and salinity) advection scheme is a mixed upstream–MUSCL (monotonic upwind scheme for conservation laws; Levy et al., 2001) scheme, as modified in Oddo et al. (2009). The vertical diffusion and viscosity terms are defined as a function of the Richardson number following Pacanowski and Philander (1981). The air–sea surface fluxes of momentum, mass, and heat are computed using bulk formulae described in Pettenuzzo et al. (2010), and the Copernicus satellite gridded sea surface temperature (SST) data (Buongiorno Nardelli et al., 2013) are used to correct the non-solar heat flux using a relaxation constant of 110 W m⁻² K⁻¹ centred at midnight. A detailed description of other specific features of the model implementations can be found in Tonani et al. (2008) and Oddo et al. (2009, 2014).

The wave model WW3 is discretized by means of 24 directional bins (15^∘ resolution) and 30 frequency bins (ranging between 0.05 and 0.7931 Hz) to represent the wave spectral distribution. The wave model is implemented using the same bathymetry and grid of the hydrodynamic model and uses the surface currents to evaluate the wave refraction but assumes no interactions with the ocean bottom. The Mediterranean implementation of WW3 follows WAM Cycle 4 model physics (Guenther et al., 1992); the wind input and dissipation terms are based on Janssen's quasi-linear theory for wind wave generation (Janssen, 1989, 1991), the wave dissipation term is based on Hasselmann's (1974) whitecapping theory according to Komen et al. (1984), and the nonlinear wave–wave interaction is modelled using the Discrete Interaction Approximation (DIA, Hasselmann et al., 1985). The exchanges between the circulation and wave models are performed using an online two-way coupling between NEMO and WW3. The models are forced by the same atmospheric fields (high-resolution ECMWF analysis and forecast winds) and are two-way coupled at hourly intervals, exchanging the following fields: NEMO sends to WW3 the sea surface currents and temperature, which are then used to evaluate the wave refraction and the wind speed stability parameter, respectively. The neutral drag coefficient computed by WW3 is passed to NEMO to compute the surface wind stress.

The NEMO–WW3 coupled system is intended to provide the representation of current–wave interaction processes in the ocean general circulation. At the moment, the feedback is considered only for the surface wind stress drag coefficient, and more details on this wave–current model coupling can be found in Clementi et al. (2017a).

2.1.2 Model initialization, external forcing, and boundary conditions

The PHY component was initialized in January 2015 using temperature and salinity winter climatological fields from WOA13 V2 (World Ocean Atlas 2013 V2, https://www.nodc.noaa.gov/OC5/woa13/woa13data.html, last access: 30 June 2023). The atmospheric forcing fields for both NEMO and WW3 models are from the $\mathrm{1}/\mathrm{8}$^∘ horizontal resolution and 6 h temporal frequency (a 3 h frequency is used to force the first 3 d of forecasting) operational analysis fields from the European Centre for Medium-Range Weather Forecast (ECMWF) Integrated Forecasting System (IFS), and a higher spatial resolution of $\mathrm{1}/\mathrm{10}$^∘ (with a higher forecast temporal frequency of 1–3–6 h according to the forecast leading time) is used starting from year 2020.

The circulation model's lateral open boundary conditions (LOBCs) in the Atlantic Ocean are provided by the Copernicus Global Analysis and Forecast product (Lellouche et al., 2018) at $\mathrm{1}/\mathrm{12}$^∘ horizontal resolution and 50 vertical levels. Daily mean fields are used, and the numerical schemes applied at the open boundaries are the Flather (1976) radiation scheme for the barotropic velocity and the Orlanski (1976) radiation condition (normal projection of oblique radiation case) with adaptive nudging (Marchesiello et al., 2001) for the baroclinic velocity and the tracers. The nesting technique is detailed in Oddo et al. (2009), who also show a marked improvement in the salinity characteristics of the Modified Atlantic Water and in the Mediterranean sea level seasonal variability. The Dardanelles Strait boundary conditions (Delrosso, 2020) consist of a merge between the Copernicus global ocean products and daily climatology derived from a Marmara Sea box model (Maderich et al., 2015). The WW3 model implementation considers closed boundaries in both the Atlantic Ocean and the Dardanelles Strait.

The river runoff inputs consist of monthly climatological data for 39 major rivers (characterized by an average discharge larger than 50 m³ s⁻¹) with a prescribed constant salinity at the river mouth (Delrosso, 2020) that is evaluated by means of sensitivity experiments and is listed in Table A4. More realistic and time-varying river salinity values (at least for major rivers) will be evaluated in future modelling evolutions using an estuary box model, such as the one presented in Verri et al. (2020), coupled to the hydrodynamic model.

2.1.3 The data assimilation component

A 3D variational data assimilation scheme, called OceanVar, initially developed by Dobricic and Pinardi (2008) and further improved for a wide range of ocean data assimilation applications (Storto et al., 2015), is coupled to NEMO.

The OceanVar scheme aims to minimize the cost function, as described in the following equation:

$$\begin{array}{}\text{(1)}& J={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\delta }{x}^T{B}^{-\mathrm{1}}\mathit{\delta }x+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(H\mathit{\delta }x-d{)}^T{R}^{-\mathrm{1}}(H\mathit{\delta }x-d),\end{array}$$

where $\mathit{\delta }x=x-x_{\mathrm{b}}$; x is the unknown ocean state, equal to the analysis x_a at the minimum of J; x_b is the background state; and $d=y-H\left(x_{\mathrm{b}}\right)$ is the misfit between an observation y and its modelled correspondent mapped onto the observation space to the observation location by the observation operator, H.

In OceanVar, the background error covariance matrix is considered to be B=VV^T, where V is a sequence of linear operators: V=V_ηV_HV_V. Multivariate EOFs (empirical orthogonal functions, described in Dobricic et al., 2008 and Pistoia et al., 2017) compose the vertical component operator, V_V. EOFs are computed in every grid point for the sea surface height, temperature, and salinity using a 3-year simulation in order to capture the mesoscale eddy variability that is assumed to represent the unbalanced component of the background error covariance. The horizontal covariances, V_H, are modelled by an iterative recursive filter (Dobric and Pinardi, 2008; Storto et al., 2014). In order to assimilate altimeter observations, the dynamic height operator, V_η, developed in Storto et al. (2011), is used. A reference level of 1000 m is used for this operator so sea level anomaly (SLA) along-track observations over water shallower than this depth are not assimilated.

The observational error covariance matrix, R, is estimated following the relationship of Desroziers et al. (2005). The assimilated observations include along-track altimeter sea level anomaly from six satellites and in situ vertical temperature and salinity profiles from Argo floats. The SLA tracks provided by nadir altimeters are assimilated by subsampling every second observation to reduce the spatial correlation between consecutive measurements. A special quality control procedure is applied to real-time Argo data before they are assimilated. It consists of removing real-time profiles with negative temperature and/or salinity, temperatures higher than 45 ^∘C, and salinities higher than 45 PSU; profiles with gaps in the observations of more than 40 m in the first 300 m depth (to avoid possible inconsistencies in the thermocline); and profiles with observations provided only below 35 m depth and observations in the first model layer (0–2 m). Moreover, a background quality check is implemented to reject observations whose square departure exceeds the sum of the observational and background error variances 64 times in the case of SLA and 25 times in the case of in situ temperature and salinity. The quality checks are applied to each individual observation of each Argo vertical profile and for each altimeter track. The misfits are computed at the observation time by applying the FGAT (First Guess at the Appropriate Time) procedure, and the corrections to the background are applied once a day to the restart file using observations within a 1 d time window.

2.2 The wind wave component

2.2.1 Numerical model description

The WAV component consists of two nested wave model implementations: the first grid covers the whole Mediterranean Sea at $\mathrm{1}/\mathrm{24}$^∘ horizontal resolution, and it is nested within a coarser resolution wave model grid at $\mathrm{1}/\mathrm{6}$^∘ horizontal resolution implemented over the Atlantic Ocean (Fig. 2).

The wave model is based on the state-of-the-art third-generation WAM Cycle 4.6.2, which is a modernized and improved version of the well-known and extensively used WAM Cycle 4 wave model (WAMDI Group, 1988; Komen et al., 1994). WAM solves the wave transport equation explicitly without any presumption regarding the shape of the wave spectrum. Its source terms include the wind input, whitecapping dissipation, nonlinear transfer, and bottom friction. The wind input term is adopted from Snyder et al. (1981). The whitecapping dissipation term is based on Hasselman's (1974) whitecapping theory. The wind input and whitecapping-dissipation source terms of the present cycle of the wave model are a further development based on Janssen's quasi-linear theory of wind wave generation (Janssen, 1989, 1991). The nonlinear transfer term is a parameterization of the exact nonlinear interactions (Komen et al., 1984; Hasselman et al., 1985). Lastly, the bottom-friction term is based on the empirical JONSWAP model of Hasselman et al. (1973).

The bathymetric map has been constructed using the GEBCO 30 arcsec bathymetric dataset for the Mediterranean Sea model and the ETOPO 2 dataset (US Department of Commerce, National Oceanic and Atmospheric Administration, NOAA National Centers for Environmental Information, 2022) 2 min gridded global relief data for the North Atlantic model. In both cases, mapping on the model grid was done using bi-linear interpolation accompanied by some degree of isotropic Laplacian smoothing. This bathymetry is different from the one used for the PHY component, optimized for the specific quality of the wave products.

The wave spectrum is discretized using 32 frequencies, which cover a logarithmically scaled frequency band from 0.04177 to 0.8018 Hz (covering wave periods ranging from approximately 1 to 24 s) at intervals of $\mathrm{d}f/f$ = 0.1 and 24 equally spaced directions (15^∘ bin). The WAV model component runs in shallow-water mode considering wave refraction due to depth and currents in addition to depth-induced wave breaking. Modifications from default values of WAM 4.6.2 have been performed in the input source functions as a result of a tuning procedure. Specifically, the value of the wave age shift parameter (ZALP¹ in the wind input source function was set to 0.011 (0.008 is the default) for the Mediterranean model, and the tunable whitecapping-dissipation coefficients C_DS and δ were altered from their default values to become C_DS=1.33 (2.1 default) and δ=0.5 (default value was 0.6). Finally, a limitation to the high-frequency part of the wave spectrum corresponding to the Cy43r1 ECMWF wave-forecasting system (ECMWF, 2016) was also implemented and tested in order to reduce the wave steepness at very high wind speeds.

2.2.2 Model initialization, external forcing, and boundary conditions

The WAV component is forced with 10 m above sea surface analysis and forecasted ECMWF winds at $\mathrm{1}/\mathrm{8}$^∘ dissemination resolution. The temporal resolution is 6 h for the analysis, 3 h for the first 3 d of the forecast, and 6 h for the rest of the forecast cycle. From year 2021, a higher spatial ($\mathrm{1}/\mathrm{10}$^∘ for both analysis and forecast) and temporal (hourly for forecast days 1–3, 3-hourly for days 4–6, and 6-hourly for days 7–10) resolution dataset is used to force the WAV component. The wind is bi-linearly interpolated onto the model grids. Sea ice coverage fields used by the North Atlantic wave model are also obtained from ECMWF. With respect to current forcing, the WAV model is forced by daily averaged surface currents obtained from Copernicus Marine Service Med-MFC at $\mathrm{1}/\mathrm{24}$^∘ resolution, and the North Atlantic model is forced by daily averaged surface currents obtained from the Copernicus global physical model at $\mathrm{1}/\mathrm{12}$^∘ resolution. The WAV component runs one cycle per day, operating in analysis (for 24 h in the past (previous day)) and forecast (for 10 d in the future) modes. During the analysis phase, the model background is blended through data assimilation with available significant wave height (SWH) satellite observations at 3-hourly intervals and forced with ECMWF analysis 6-hourly winds and daily averaged surface currents.

The Mediterranean Sea model receives a full wave spectrum at 5 min intervals at its Atlantic Ocean open boundary from the WAM implementation in the North Atlantic. The latter model is considered to have all four of its boundaries closed, assuming no wave energy propagation from the adjacent seas. This assumption is readily justified for the north and west boundaries of the North Atlantic model considering the adjacent topography which restricts the development and propagation of swell into the model domain.

2.2.3 The wave data assimilation component

The assimilation module of the WAV component is based on the data assimilation scheme of WAM Cycle 4.6.2, which consists of an optimal interpolation (OI) of the along-track significant wave height (SWH) observations retrieved by altimetry and then a re-adjustment of the wave spectrum at each grid point accordingly. This assimilation approach was initially developed by Lionello et al. (1992) and consists of two steps. First, a best-guess (analysed) field of significant wave height is determined by OI with appropriate assumptions regarding the error covariance matrix. One of the key issues is the specification of the background error covariance matrix for the waves called P and the observation error covariance matrix, R. The first is defined as in the following equation:

$$\begin{array}{}\text{(2)}& \mathbf{P}=\exp \left({\displaystyle \frac{d_{ij}}{l_{\mathrm{c}}}}\right),\end{array}$$

while the second is defined according to Eq. (3) as follows:

$$\begin{array}{}\text{(3)}& \mathbf{R}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{2}}}},\end{array}$$

where i and j are the model grid points in the longitudinal and latitudinal directions respectively, d is the distance of the observation location to the grid point, l_c is the field correlation length, and ${\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{2}}$ stand for the observation and model errors respectively. In the above expressions, the error is considered to be homogeneous and isotropic. We use R=1 and the correlation length l_c equal to 3^∘ (∼ 300 km).

Finally, the weights assigned to the observations are the elements of the gain matrix K, as presented in Eq. (4):

$$\begin{array}{}\text{(4)}& \mathbf{K}=\mathbf{P}{H}^T[H\mathbf{P}{H}^T+\mathbf{R}{]}^{-\mathrm{1}},\end{array}$$

where H is the observation operator that projects the model solution to the observation location. For the current version of Med-waves, the OI analysis procedure is applied only to altimeter along-track SWH measurements, although wind at 10 m measurements can be assimilated as well. Prior to the OI procedure, quality-checked SWH observations which are available in a ±1.5 h time window are collocated with the closest model grid point and are averaged.

During the second step, the analysed significant wave height field is used to retrieve the full dimensional wave spectrum from a first-guess spectrum provided by the model itself, introducing additional assumptions to transform the information of a single wave height spectrum into separate corrections for the wind sea and swell components of the spectrum. Two-dimensional wave spectra are regarded either as wind sea spectra if the wind sea energy is larger than 3–4 times the total energy, or, if this condition is not satisfied, as swell. If the first-guess spectrum is mainly wind–sea, the spectrum is updated using empirical energy growth curves from the model. In the case of swell, the spectrum is updated assuming the average wave steepness provided by the first-guess spectrum is correct, but the wind is not updated.

Prior to assimilation, all altimeter SWH observations are subject to a quality control procedure. Every day, the system is scheduled to simulate 264 h: 24 h in the past (analysis) blending through data assimilation model results with all satellite SWH observations available followed by a 240 h forecast. The assimilation step adopted for the current version of the Med-waves system is equal to 3 h.

2.3 Mediterranean biogeochemical component

2.3.1 Numerical model description

The BIO component consists of the Biogeochemical Flux Model (BFM, Vichi et al., 2007) coupled with the transport (OGSTM) module (Salon et al., 2019). Advection, vertical and horizontal diffusion, and the sinking term for the biogeochemical tracers (Foujols et al., 2000) are solved by the OGSTM module that uses daily 3D velocity, diffusivities, and 2D atmospheric fields provided by the PHY component through the offline coupling scheme (Fig. 1). A source-splitting numerical time integration is used to couple advection and diffusion to the biochemical tracer rates.

BFM describes the biogeochemical cycles of carbon, nitrogen, phosphorus, silicon, and oxygen through the dissolved inorganic compartment and the particulate living and non-living organic compartments (Lazzari et al., 2012, 2016). The model includes four phytoplankton functional groups (i.e. diatoms, flagellates, picophytoplankton, and dinoflagellates), four zooplankton groups (i.e. carnivorous and omnivorous mesozooplankton, heterotrophic nanoflagellates, and microzooplankton) and heterotrophic bacteria. Among the nutrients, dissolved inorganic nitrogen is simulated in terms of nitrate and ammonia. The non-living dissolved organic compartment includes labile, semi-labile, and refractory organic matter. A carbonate system component (Cossarini et al., 2015) includes alkalinity (ALK), dissolved inorganic carbon (DIC), and particulate inorganic carbon (PIC) as prognostic variables; computes CO₂ air–sea gas exchange according to Wanninkhof (2014); and provides diagnostics variables such as pH, CO₂ concentration, and calcite saturation horizon.

2.3.2 Model initialization, external forcing, and boundary conditions

Initial conditions of nutrients (nitrate, ammonia, silicate, and phosphate), oxygen, and carbonate variables (DIC and alkalinity) consist of 16 climatological profiles homogeneously applied in each of the subregions represented in Fig. 3. Climatological profiles are computed from the European Marine Observation and Data Network (EMODnet) dataset (Buga et at., 2018). The other biogeochemical state variables (phytoplankton, zooplankton, and bacteria biomasses) are initialized in the photic layer (0–200 m) according to the standard BFM values. A 5-year hindcast is run using the first year (i.e. 2017) in perpetual mode. The model has two open lateral conditions: in the Atlantic Ocean and in the Dardanelles Strait. Nutrients, oxygen, DIC, and alkalinity in the Atlantic (i.e. boundary at long = 9^∘ W) are provided through seasonally varying climatological profiles derived from the World Ocean Atlas (Boyer et al., 2018) and the literature (Alvarez et al., 2014), and a Newtonian dumping is applied. The Newtonian dumping is set between the longitudes 9^∘ and 6.5^∘ W with a timescale relaxation term linearly varying from 1/24 per day at 9^∘ W to 90 per day at 6.5^∘ W. A Dirichlet-type scheme with constant concentration values of nutrients, DIC, and alkalinity derived from the literature (Yalcin et al., 2017; Tugrul et al., 2002; Souvermezoglou et al., 2014; Copin-Montegut, 1993; Schneider et al., 2007; Krasakopoulou et al., 2017) is applied at the Dardanelles Strait. The concentrations are also tuned to provide input fluxes from the Black Sea to the Mediterranean Sea that are consistent with published estimates (Deliverable of Perseus, 2012; Yalcin et al., 2017; Tugrul et al., 2002; Copin-Montegut, 1993). A radiative condition is set for the other BFM tracers.

Figure 3 The Mediterranean Sea domain and subregion subdivisions for analysis of the skill scores: Alboran (1), southwest Mediterranean-1 (2), northwest Mediterranean (3), southwest Mediterranean-2 (4), south Tyrrenian (5), north Tyrrenian (6), west Ionian (7), east Ionian (8), northeast Ionian (9), south Adriatic (10), north Adriatic (11), west Levantine (12), Aegean (13), north central Levantine (14), south central Levantine (15), east Levantine (16).

Terrestrial inputs include 39 rivers, consistent with the PHY component. Annual nutrients input are about 46 500 × 10⁶ molN yr⁻¹ and 881 × 10⁶ molP yr⁻¹ (Salon et al., 2019). Carbon and alkalinity inputs are 9300 × 10⁹ gC yr⁻¹ and 800 × 10⁹ mol yr⁻¹ respectively. Estimates are derived considering typical concentrations per freshwater mass in macro-coastal areas of the Mediterranean Sea (Copin-Montegut, 1993; Meybeck and Ragu, 1995; Kempe et al., 1991) and the river water discharges from the Perseus dataset (Deliverable of Perseus, 2012, as before). Annual atmospheric nutrient depositions are 81 300 × 10⁶ molN yr⁻¹ and 1194 × 10⁶ molP yr⁻¹ for nitrogen and phosphorus respectively (Ribera d'Alcalà et al., 2003). Spatially constant values of atmospheric pCO₂ are derived from the 1992–2018 time series of the ENEA Lampedusa station (Trisolino et al., 2021), with the 2019 and 2020 values extrapolated by linear trend.

2.3.3 The biogeochemical data assimilation component

The BIO component features a variational data assimilation scheme (3DVarBio) which is based on the minimization of the cost function (Eq. 1) (Teruzzi et al., 2014). Minimization is computed iteratively in a reduced space using an efficient parallel PETSc/TAO (Portable, Extensible Toolkit for Scientific Computation and the Toolkit for Advanced Observation) solver (Teruzzi et al., 2019), and the background error covariance matrix, B, is factored as B=VV^T, where V is a sequence of linear operators: V=V_BV_HV_V. The horizontal error covariance operator (V_H) is a Gaussian filter and includes a non-uniform and direction-dependent length-scale correlation radius to account for anisotropic coastal assimilation (Terruzzi et al., 2018) and vertical profile assimilation (Cossarini et al., 2019). The vertical error covariance operator (V_V) is based on a set of 0–200 m vertical error profiles obtained using an empirical orthogonal function (EOF) decomposition of a 20-year-long pre-existing biogeochemical simulation. EOFs are computed monthly for the 16 subregions with the actual vertical resolution and are rescaled at each grid point considering the ratio between observation and model variances (Teruzzi et al., 2018). The biogeochemical error covariance operator (V_B) is designed to preserve the ratios among phytoplankton functional types and their internal carbon-to-nutrient quotas (Teruzzi et al., 2014) and supports monthly and spatially varying covariances between dissolved inorganic nutrients (Teruzzi et al., 2021). In the most recent BIO model configuration (Teruzzi et al., 2021; Cossarini et al., 2019), the assimilated biogeochemical observations are satellite multi-sensor (MODIS, VIIRS (Visible Infrared Imaging Radiometer Suite), and OLCI – Ocean and Land Colour Imager) surface chlorophyll data (Volpe et al., 2019) and quality-controlled BGC-Argo nitrate and chlorophyll profiles (Schmechtig et al., 2018; Johnson et al., 2018). Ocean colour data are interpolated from original 1 km resolution to the $\mathrm{1}/\mathrm{24}$^∘ model resolution.

2.4 System evolutions

The Mediterranean has been the site of major forecasting research activities since the late nineties (Pinardi et al., 2003; Pinardi and Coppini, 2010). The changes and evolution of the MED-MFS components are presented in Table 1. Before 2008, only the PHY and BIO components were present. The PHY component was based on the Ocean Parallelise (OPA) code (Madec et al., 1998) with the highest available horizontal and vertical resolutions of $\mathrm{1}/\mathrm{16}$^∘ (approx. 6.5 km) in the horizontal and 72 vertical levels, with closed lateral boundaries; only seven major rivers; and the implementation of a weekly 3D-VAR assimilation scheme (Dobricic et al., 2007), assimilating temperature and salinity vertical profiles, sea level anomaly (SLA), and track altimeter data. Moreover, a non-solar heat flux correction was imposed through a nudging for the whole day with sea surface temperature (SST) satellite gridded data.

A major upgrade of the PHY component was achieved in 2009 by implementing a version of the numerical model NEMOv3.1 that included LOBCs in the Atlantic Ocean (Oddo et al., 2009) and by moving to a daily assimilation cycle. The first exchanges with a wave model were implemented in 2010 when the PHY component was coupled hourly with WAM, receiving the surface drag coefficient to better represent the wind stress. In 2013, the whole operational modelling system was updated by implementing an upgraded two-way online coupled system based on NEMOv3.4 and WW3 (Clementi et al., 2017a), allowing for a more consistent exchange between the two models. The following year, the PHY general circulation module was improved by accounting for the effect of atmospheric pressure (in addition to wind and buoyancy fluxes) and an explicit linear free-surface formulation using a time-splitting scheme (Oddo et al., 2014), while the assimilation scheme was enhanced thanks to the assimilation of Tailored Altimetry Products for Assimilation Systems (TAPAS) SLA data, allowing for the application of specific corrections of the altimetric original signal (Dobricic et al., 2012).

The PHY component delivered in 2015 included the nesting in the Atlantic Ocean through daily analysis and forecast fields from the global system, while 1 year later, the assimilation scheme was enhanced by including the computation of monthly and grid point EOFs and vertical observational errors varying with depth.

Another major PHY component evolution was achieved in 2017 when the resolution of the operational system was increased to $\mathrm{1}/\mathrm{24}$^∘ (approx. 4 km) horizontal resolution and 141 vertical levels using the z-star vertical coordinate system, a nonlinear free-surface formulation, and when the NEMOv3.6 version and 39 rivers were introduced. From the year 2019, the Dardanelles Strait inflow was set as a lateral open boundary condition (instead of as a river runoff climatological input), allowing for a daily update of the fluxes, and an improved nudging with the satellite sea surface temperature was included by correcting the heat fluxes only close to midnight.

The WAV component was developed and released for the first time in 2017 based on WAM Cycle 4.5.4, providing on a daily basis 5 d wave forecasts and simulations for the Mediterranean Sea at $\mathrm{1}/\mathrm{24}$^∘ horizontal resolution (Ravdas et al., 2018) nested within a North Atlantic model at $\mathrm{1}/\mathrm{6}$^∘ resolution and forced with ECMWF 10 m winds and PHY component surface currents. In March 2018, the system was upgraded by incorporating the data assimilation component to utilize available track SWH satellite observations from Sentinel-3A and Jason-3. In 2019, the wave model was upgraded to Cycle 4.6.2, and the durations of the forecasts were extended to 10 d. Additionally, a limitation to the high-frequency part of the wave spectrum was applied, while modifications from default values were introduced in the input source and dissipation functions: ZALP was set to 0.011, and C_DS and δ became 1.33 and 0.5 respectively.

In 2009, the first pre-operational version of the BIO component featured early versions of the OGSTM transport model and the BFM model (Lazzari et al., 2010). The spatial resolution was $\mathrm{1}/\mathrm{8}$^∘, which required a subsampling of the PHY component fields from the $\mathrm{1}/\mathrm{16}$^∘ resolution. The Atlantic boundary was closed with a nudging term for nutrients, and the land nutrients input included the three major Mediterranean rivers (i.e. Po, Rhone, and Nile), and the Dardanelles was treated as a river. BFM used constant daily averaged irradiance to force photosynthesis (Lazzari et al., 2010).

Horizontal resolution aligned with the physical model in 2013 and was refined to $\mathrm{1}/\mathrm{24}$^∘ in 2017. Full alignment between the PHY and BIO components in terms of the same horizontal and vertical resolutions, bathymetry, and boundaries (number and position of rivers) was introduced in 2018 and remained a standard that mitigated possible approximation errors related to the use of daily outputs of the eddy-resolving ocean general circulation model to force the transport of tracers (Salon et al., 2019). Additionally, nutrient and carbon land inputs from 39 rivers were introduced in 2017, and open boundary conditions were introduced in Dardanelles Strait in 2019 and in the Atlantic Ocean in 2020 (Salon et al., 2019).

Since 2008, three major improvements of the BFM model have been integrated: (i) the addition of the carbonate system to predict alkalinity, ocean acidity, and CO₂ air–sea exchanges in 2016 (Cossarini et al., 2015); (ii) the revision of nutrient formulation of phytoplankton in 2018 (Lazzari et al., 2016); and (iii) in 2020, the introduction of the day–night cycle in the light-dependent formulation of phytoplankton (Salon et al., 2019) and of the novel light extinction coefficient (Terzic et al., 2021).

A major system evolution and quality improvement was achieved in 2013 with the inclusion of the assimilation of satellite chlorophyll through a variational scheme with prescribed background error covariance (Teruzzi et al., 2014). The assimilation method was improved in 2018 to include coastal components (i.e. non-uniform and direction-dependent horizontal covariance; Teruzzi et al., 2018) and in 2019 to integrate new observations (i.e. BGC-Argo float profiles), including a new parameterization for the vertical and biogeochemical background error covariance (Cossarini et al., 2019).

In terms of operational product delivery, the BIO component has produced daily 10 d forecasts and weekly 7 d analyses since 2020, fully aligned with the PHY component (Salon et al., 2019). Before that, the system produced 7 d analyses and a 7 d forecast once per week since 2013, while a second cycle of 7 d forecasts was added each week in 2015.

3 Quality assessment

The evaluation of the quality of the Med-MFC is given here only for the analysis products, leaving the assessment of the forecast skill for future work. One overarching driver for the Med-MFC evolution is the continuous improvement of the numerical model and data assimilation modules with respect to a well-defined set of goodness indices established for all the European regional seas (Hernandez et al., 2009). Ocean model uncertainties emerge from sources of error relevant to the ocean state, including physics, biogeochemistry, and sea ice, as well as errors in the initial state and boundary conditions (i.e. atmospheric forcing and lateral open boundary conditions). Model uncertainties in ocean physics have a significant impact on all other system components, as in, for example, biogeochemistry and sea ice (Alvarez Fanjul et al., 2022). Our results describe the quality of the Med-MFC products, presenting the statistics and accuracy numbers based on a reference simulation produced to calibrate and validate the operational forecasting systems, whereas the analysis of model uncertainty sources is outlined in the discussion part with reference to previous specific publications.

3.1 PHY component skill

The skill of the physical component is assessed over a 3-year period from 2018 to 2020 (Clementi et al., 2019). The evaluation is done by means of estimated accuracy numbers (EANs), which consist of the root mean square difference (RMSD) and bias (model minus observations) of the daily mean of model outputs against satellite and in situ observations. EANs are evaluated using the daily mean of model estimates interpolated on the available observations for that day; this goodness score is somewhat approximated, especially at the surface where daily variability is large, but this is a score used by many forecasting systems (Ciliberti et al., 2022; Toledano et al., 2022; Sotillo et al., 2021; Nagy et al., 2020), and we will show it for reference purposes. We also use misfits, which are the differences between the model solutions and the observations at the observational time during the forward model integration, for this assessment. The misfits provide quasi-independent and more accurate skill assessment since they are calculated before the variational analysis and at the observational time.

Table 1 Changes in the Mediterranean forecasting components since 2008.

^* The sea level anomaly (SLA) TAPAS product is produced to give information about the different corrections of the altimetric original signal.

Table 2 EAN estimates with in situ observations. The differences (bias) and their squared values (RMSD) are then averaged over the whole Mediterranean Sea region and over nine vertical layers for the years 2018–2020.

Table 2 summarizes the EAN of 3D model temperature and salinity daily mean values compared to in situ observations, in particular Argo floats and CTD profiles averaged over the 3 reference years. Model temperature shows small positive and negative biases depending on the depth, with the largest error (maximum value of the period is 0.85 ^∘C) in the sub-surface layers between 10 and 60 m, decreasing with depth. Salinity is characterized by an almost generally negative small bias, meaning generally lower salinities than measured, along the whole water column except for the first layer. The salinity RMSD mean value is generally lower than 0.2 PSU, and the error is larger in the first layers and decreases significantly below 150 m. The comparison with other Copernicus Marine Service forecasting-system EAN values presented in the Quality Information Document (QUID), considering the fact that the validation periods are different, shows that the Mediterranean temperature and salinity qualities in terms of RMSD are aligned with all the other Copernicus forecasting systems. In particular, the average sea surface temperature RMSD with respect to satellite data ranges from 0.48 ^∘C in the Northwest Shelf (derived from the QUID of the product NORTHWESTSHELF_ANALYSIS_FORECAST_PHY_004_013 https://doi.org/10.48670/moi-00054) to 0.8 ^∘C in the Baltic Sea (derived from the QUID of the product BALTICSEA_ANALYSISFORECAST_PHY_003_006 https://doi.org/10.48670/moi-00010), while the 3D mean temperature RMSD with respect to in situ data ranges from 0.4 ^∘C in the Mediterranean and Northwest Shelf to 0.7 ^∘C in the Black Sea (derived from the QUID of the product BLKSEA ANALYSISFORECAST_PHY_007_001 https://doi.org/10.25423/cmcc/blksea_analysisforecast_phy_007_001_eas4), and the salinity mean RMSD varies from 0.1 PSU in the Mediterranean and Northwest Shelf to 0.3 PSU in the Iberia–Biscay–Ireland area (derived from the QUID of the product IBI_ANALYSISFORECAST_PHY_005_001 https://doi.org/10.48670/moi-00027). The sea level anomaly skill is also aligned with those of other operational systems within the Copernicus Marine Service when compared with satellite altimeter observations (from 2.2 cm in the Black Sea to 9 cm in the Northwest Shelf area).

The other goodness index is computed as weekly mean root mean square error and bias using temperature and salinity misfits that are computed at FGAT. The misfits are more precise in accounting for surface errors since the observations are compared with the model at the exact time of the day when observations are taken. This index is represented as a depth–time Hovmöller diagram in Fig. 4. The temperature error is seasonal (Fig. 4a), with maximum values of ∼ 1.8 ^∘C in the range of 30–60 m depth, corresponding to the depth of the mixed layer and the seasonal thermocline during the stratified season from June to November. The error is reduced to an average value of around 0.4 ^∘C during the vertically mixed season from December to May. The temperature misfits (Fig. 4c) indicate an overall overestimation of the temperature, except for in the subsurface layers, during winter and spring.

Figure 4 Hovmöller (depth–time) diagrams: (a) weekly root mean square (rms) of temperature misfits, (b) weekly rms salinity misfits, (c) weekly bias of temperature, and (d) weekly bias of salinity evaluated along the water column and averaged over the whole Mediterranean Sea.

The salinity error (Fig. 4b) is defined by two main structures: one that is constant throughout the year down to about 150 m and the seasonal amplification during summer, as for the temperature errors. The maximum errors reach values of 0.35 PSU in the summer period and decrease to 0.025 PSU below ∼ 150 m. We argue that the background error, uniform throughout the year, could be due to inaccurate advection of salinity in different sub-areas of the Mediterranean Sea. Moreover, the model salinity bias is generally negative; i.e. the model salinity is lower than the observations (Fig. 4d). This could be related to the larger Atlantic water inflow with respect to the literature (Soto-Navaro et al., 2010) at Gibraltar, as reported in Table 3, and to inaccurate mixing at Gibraltar due to the lack of tides.

Table 3 Gibraltar mean and standard deviation volume transports [Sv] from the Med-PHY numerical system averaged over the period 2018–2020 compared to literature values (current metre observations from October 2004 to January 2009).

Sea surface temperature (SST) and sea level anomaly (SLA) skills are evaluated by comparing them with satellite observations: model daily mean SST is compared to SST satellite L4 gridded data at $\mathrm{1}/\mathrm{16}$^∘ resolution (Buongiorno Nardelli et al., 2013), while SLA is compared to along-track satellite altimeter observations (Taburet et al., 2019) in terms of model misfits. Table 4 presents the RMSD and bias values computed for SST, as well as SLA RMSD averaged in the Mediterranean Sea and over the 16 subregions (see Fig. 3). Considering SST, the RMSD values range between 0.47 and 0.69 ^∘C (mean Mediterranean Sea error is 0.54 ^∘C), and the bias is generally positive, possibly caused by an overestimation of the downward shortwave radiation flux, which is estimated according to Reed's (1977) formula, as already discussed in Byun and Pinardi (2007) and Pettenuzzo et al. (2010). The SLA error ranges between 2.3 and 5.3 cm (mean error is 3.8 cm). The SLA skill scores vary in different regions; this could be related to the spatial coverage of the observations (not homogeneous in the basin) and to the limit of the 1000 m assimilation depth (due to the dynamic height operator, which assumes a level of no motion to compute the sea level increments from temperature and salinity increments; see Sect. 2.1.3).

Table 4 EAN RMSD and bias of SST and SLA RMSD averaged over the whole of the Mediterranean Sea and the 16 subregions (see Fig. 3) for the period 2018–2020.

The time variability of the model SLA accuracy is also provided by means of weekly model misfits evaluated for each available satellite altimeter and averaged over the whole Mediterranean Sea, as shown in Fig. 5. The error ranges between 2.5 and 5.5 cm (maximum error with respect to CryoSat) with a large variability among the different satellites and with a generalized increase of error during autumn and winter seasons.

Figure 5 Time Series of weekly mean rms misfit error for SLA evaluated with respect to available satellite altimeters and averaged over the whole Mediterranean Sea. The bold black line represents the mean error with respect to the whole set of satellites, which are separately shown with different colours. The grey area indicates the number of observations used for the validation.

3.2 WAV component skill

The quality of the wave analysis and forecast product is assessed over a 3-year period from January 2018 to December 2020. The skill of the Mediterranean wave model is assessed by considering intercomparisons of the model solution during the 24 h analysis phase with available in situ (SWH and mean wave period from wave buoys) and remotely sensed (SWH) observations. As the latter are ingested into the model through data assimilation, the model first-guess SWH (i.e. model background) is used instead of model analysis.

Significant wave height (SWH) and mean wave period (MWP) measurements are used for data validation from 28 wave buoys in the Mediterranean Sea (lower panel of Fig. 7). Data quality control procedures have been applied to the in situ observations (Copernicus Marine In Situ Team and Copernicus In Situ TAC, 2020), and measurements associated with a bad quality flag are not taken into consideration.

Figure 6 Scatter plots of (a) significant wave height (H_s) and (b) mean wave period (T_m) versus wave buoy observations for the 28 stations of the Mediterranean Sea (Fig. 7c) for a 3-year period (2018–2020). The graphs also include quantile–quantile plots (grey crosses), 45^∘ reference lines (dashed orange line), and least-squares best-fit lines (green line). At the top left of each picture, statistical scores are given: entries refer to the amount of data available for computing the statistics, R and M refer to the observed and modelled values respectively, SI is the scatter index (defined as the standard deviation of model–observation differences relative to the observed mean), and CORR is the Pearson correlation coefficient.

Figure 7 Significant wave height difference between model and observations (a, b) at the 28 buoy locations (c) for a 3-year period (2018–2020). For all locations, the performance of the model is evaluated against buoy data by means of bias, root mean square difference (RMSD), scatter index (SI), and deviations of the Pearson correlation coefficient from unity (CORR deviation).

Figure 6 depicts scatter plots of the evaluation of the observed SWH and MWP against measurements obtained from the 28 buoys. For the immense amount of match-up data (within the range 0–1.25 m), the model overestimates SWH with respect to the buoy measurements (left side panel). Additionally, the model underestimates SWH during more energetic events (>1.25 m), except in the range 5.5–6.2 m. For large wave heights, model results underestimate SWH compared to the buoys, which agrees with past findings for the Mediterranean Sea (Ardhuin et al., 2007; Korres et al., 2011). Negative SWH bias can be attributed to errors in the forcing or inaccurate wave growth and dissipation at high wind speeds (Pineau-Guillou et al., 2018). The dashed orange line (i.e. the 45^∘ reference line) in the quantile–quantile (Q–Q) plot stands for the unit gradient line. We observe that model results follow the dashed orange line very closely, meaning the model produces well the distribution of SWH observations. Although for higher waves (>1.25 m) the model tends to underestimate SWH (except in the range 5.5–6.2 m), it overproduces very large wave heights (100th, 99.97th, 99.96th, and 99.95th percentiles); hence, a deviation from the dashed orange reference line in the Q–Q plot becomes prominent for very high waves. Concerning MWP, the model systematically underestimates it (right side panel). Despite the overall modelled MWP underestimation (bias = −0.314 s), the system tends to overestimate MWP for high percentiles and/or very long waves (hence, we observe the deviation of the Q–Q plot from the unit gradient line for very high periods). Seasonal results (not shown) for both variables SWH and MWP indicated that the model adequately captures the seasonal variability. For SWH, RMSD values vary from 0.154 m in summer to 0.231 m in winter. Nevertheless, the scatter index (SI) is higher in summer (0.26) than during the other seasons. Additionally, the highest Pearson correlation coefficient (CORR) is observed in winter (0.963), while the lower one is equal to 0.932 and is observed in summer. The metrics reveal that the model better follows the observations in winter than during the other months since the former is associated with more well-defined weather patterns and higher waves. A similar conclusion has also been reached in other studies (e.g. Ardhuin et al., 2007) for the Mediterranean Sea. Summer and autumn are characterized by higher SI values (0.244 and 0.260 respectively), while lower values are obtained for winter and spring (0.231 and 0.227 respectively). Finally, small positive bias values are met for all seasons, with the highest values found in summer (0.012 m). Regarding mean wave period, RMSD varies from 0.610 s in summer to 0.66 s in winter, and bias is negative for all seasons. SI does not present significant seasonal variability, with the highest value encountered in summer. Finally, CORR for MWP is higher than 0.8 in all seasons (values are within the range 0.859–0.878, while during summer, CORR equals 0.792). These metrics demonstrate that the model wave period (similarly to the wave height) correctly follows the observations in well-defined weather conditions characterized by higher waves and longer periods, agreeing with past studies (Cavaleri and Sclavo, 2006; Ravdas et al., 2018).

The qualification metrics for the different buoy locations in Fig. 6 are plotted in Fig. 7 (upper panel). RMSD at the different buoy locations varies from 0.13 to 0.31 m. SI varies from 0.17 at buoy 3732621 to 0.35 at the buoys of Malaga and SARON (Aegean Sea). In general, SI values above the mean value for the whole Mediterranean Sea (0.24) are obtained at wave buoys located near the coast, particularly if these are sheltered by land masses to the north–northwest (e.g. western French coastline) and/or within enclosed basins characterized by a complex topography such as the Aegean Sea. As explained in several studies (Ravdas et al., 2018), in these cases, the spatial resolution of the wave model is often not adequate to resolve the fine bathymetric features, whilst the spatial resolution of the wind forcing is incapable of reproducing the fine orographic effects, introducing errors to the wave analysis. The Pearson correlation coefficient (CORR) mostly follows the pattern of variation of SI (in this figure, we present the CORR deviation from unity). CORR ranges from 0.87 at SARON in the Aegean Sea to 0.97 at the deep-water buoy 6100196 offshore Spain, which is well exposed to the prevailing northwesterly winds in the region. The bias varies from −0.13 m at buoy 3732621 (located north of Crete) to 0.13 m at buoy 6100021 located near the French coast. Its sign varies, with positive and negative values computed at almost the same number of locations respectively. Figure 8 (right) shows the scatter plot between the first-guess SWH and satellite observations. Here, the initial-guess SWH refers to the model SWH before data assimilation, thus meaning semi-independent model data. In addition, a scatter plot resulting from the comparison of the ECMWF forcing wind speeds (U₁₀) and satellite measurements of U₁₀ is shown in Fig. 8 (left). It is seen that ECMWF forcing overestimates U₁₀ with respect to observations throughout most of U₁₀ range, while some underestimation is observed for high wind speeds (14–19 m s⁻¹). An overall ECMWF overestimation of 3 % is computed. On the other hand, the SWH model underestimation is about 6 %. Compared to the equivalent results obtained from the model–buoy comparison, a smaller scatter (by about 7 %) with a larger overall bias is associated with the model–satellite comparison, i.e. open-ocean waves. SI values compare well at the more exposed wave buoys in the Mediterranean Sea.

Figure 8 Scatter plots of ECMWF forcing wind speed U₁₀ versus satellite U₁₀ observations (a) and model significant wave height (H_s) versus satellite observations (b) over the entire Mediterranean basin for the 3-year period (2018–2020).

Figure 9 SWH evaluation against satellite data: maps of scatter index (SI) (a) and relative bias (b) over the Mediterranean Sea subregions (shown in Fig. 3) for the 3-year period (2018–2020).

Figure 9 maps statistics of the comparison of model first guesses and satellite observations of SWH for the different subregions of the Mediterranean Sea. The Aegean and Alboran seas have relatively high SI values (0.21). The highest value of SI is obtained for the north Adriatic Sea (0.26) followed by the south Adriatic (0.23). The lowest values (0.13–0.15) are found in the Levantine Basin, the Ionian Sea, and the southwest Mediterranean Sea. Relatively low values (0.16) are also found west of the islands of Sardinia and Corsica. As discussed above, the error is due to inaccuracies associated with orographic winds and/or local sea breezes and the missing representation of the complicated bathymetry in the fetch-limited, enclosed regions. SWH negative bias is present in all subregions.

Finally, intercompared to ECMWF, UK Met Office, and DMI (Danish Meteorological Institute) wave forecasting systems for a different year (2014), Med-waves shows a better skill in terms of SWH, with RMSEs for the western Med buoys equal to 0.227 m (0.234 m for ECMWF and 0.281 m for UK Met Office) and 0.201 m for the central and eastern Mediterranean (0.227 m for ECMWF and 0.268 m for DMI).

3.3 BIO component skill

The BIO component state variables can be validated at three different uncertainty levels, providing degrees of confirmation (Oreskes et al., 1994) of different scales of variability based on the availability of reference data.

Near-real-time satellite and BGC-Argo float data provide a rigorous skill performance validation dataset down to the scales of the week and mesoscale dynamics for a limited set of variables: chlorophyll, nitrate, and oxygen. Datasets of historical oceanographic data (SOCAT dataset, Baker et al., 2016; EMODnet data collection, Buga et al., 2018; Cossarini et al., 2017; Lazzari et al., 2016) are used to build a reference framework of subregions and annual and seasonal climatological profiles to validate model performance to simulate the basin-wide gradients, the mean vertical profiles, and the seasonal cycle. For these datasets, it is possible to have nutrients such as nitrate, phosphate, ammonia, and silicate, as well as dissolved oxygen, dissolved inorganic carbon, alkalinity, and surface pCO₂.

Lastly, a third level of validation regards those variables whose observability level is very scarce (e.g. phytoplankton biomass) or based on indirect estimations (e.g. primary production, air–sea CO₂ fluxes). Only a confirmation of the range of variability and a general uncertainty estimation can be provided for those variables (see, for example, the validation of model primary production in von Schuckmann et al., 2020; Cossarini et al., 2020).

Considering the 2018–2020 reference period, the chlorophyll is very well reproduced by the BIO component, both in terms of seasonal cycle and spatial gradient at the surface (Fig. 10) and in terms of vertical profiles at the BGC-Argo float positions (Table 5). The uncertainty of the surface chlorophyll is lower than 0.03 mg m⁻³ with larger values registered in winter and western subregions where the variability and the chlorophyll values are higher (Fig. 10a and b). Regarding profiles, chlorophyll values and vertical shapes driven by mesoscale dynamics are simulated with a high level of accuracy by the model (Salon et al., 2019; Cossarini et al., 2019, 2021). Daily values of RMSD and of Pearson correlation are computed between satellite and model output maps and are then averaged over the two periods (Fig. 10c and d). The plot of RMDS (Fig. 10c) shows that higher errors are registered in the western subregions and in winter when chlorophyll levels and variability are higher. On the other hand, spatial correlation values are moderate and high in all subregions (i.e. values always above 0.5, except for a few subregions), with summer values being better than winter values. Considering the number of grid points in each subregion, all values in Fig. 10d should be considered to be significantly non-zero at the 0.05 level. Indeed, Salon et al. (2019) show how, using novel metrics, the BIO component reproduces with a high level of accuracy not only the concentrations in the euphotic layer but also the seasonal evolution of the shape of the profiles. The depth of the deep chlorophyll maximum during summer and of the surface bloom during winter, as well as the depth of the nitracline and the depth of the maximum oxygen layer, which result from the interaction of physical and biogeochemical processes, are reproduced with an uncertainty of O(10¹) m (Table 5). However, the conclusions about the mesoscale accuracy of the BIO component should be taken with caution since the BGC-Argo observations are still relatively few in number (about one-eighth of the Argo floats have biochemical sensors and are unevenly spaced (e.g. the southern Mediterranean Sea is less observed than northern areas).

Table 5 The rms of the difference between MedBFM and Argo-BGC profiles for ecosystem metrics. RMSDs of the metrics are computed for each profile, then averaged over time and space considering the 2017–2020 period. Subregions: swm (reg2 + reg4), nwm (reg3), tyr (reg5 + reg6), adr (reg10 + reg11), ion (reg7 + reg8 + reg9), and lev (reg13 + reg14 + reg15 + reg16).

Figure 10 Time series of surface chlorophyll for centred composite 7 d satellite analysis (green) and the model analysis (black) in two selected subregions (a, b). The rms of differences (c) and Pearson correlation (d) between maps of satellite and model forecast for the day before the assimilation in the 16 subregions of Fig. 5c. Metrics are averaged over the winter (from October to April) and summer (from May to September) periods.

As explained above, an additional verification of biogeochemical variables can be achieved for an additional seven variables (not considering chlorophyll) and two other derived variables with climatological data. An example of such a comparison is shown in Fig. 11 for the carbonate system variables. Average maps and profiles of alkalinity and DIC in selected subregions in the zonal directions (coloured lines) are superimposed well onto the range of variability in terms of the historical in situ data (grey-shaded areas), demonstrating the capability of the BIO component to reproduce both horizontal basin-wide gradients and vertical profiles in the different areas. A slight overestimation of DIC and alkalinity (underestimation of alkalinity) is simulated in the Alboran subregion in the upper 0–100 m layer.

Figure 11 Spatial distribution of modelled alkalinity (a) and DIC (b) and comparison of vertical profiles of alkalinity (c) and DIC (d) for model (average and range of variability represented by solid and dashed coloured lines respectively) and EMODnet climatology (average and range of variability represented black dots and lines and grey-shaded areas, respectively) for selected macro areas. Climatological data are computed using historical data (Buga et al., 2018; Bakker et al., 2014). The range of variability is the average ± standard deviation.

As a summary of the skill performance analysis, statistics based on RMSD for all the considered model variables (Table 6) report the model uncertainty in reproducing the basin-wide values and gradients for the selected layers. Generally, larger errors are computed for the upper layers where the variability (both spatial and temporal) is higher. Ammonia also indicates high errors in subsurface layers, which is due to a possible incorrect initialization of deep layers as a result of the lack of data in 9 out of 16 subregions. These numbers, which were found in response to the request for a synthetic measurement of the Copernicus Marine Service product accuracy (Hernandez et al., 2018), are consolidated by means of deep skill performance analysis of the BFM model in terms of reproducing chlorophyll (Lazzari et al., 2012; Teruzzi et al., 2018), nutrients (Lazzari et al., 2016; Salon et al., 2019), and carbonate system variables (Cossarini et al., 2015).

Table 6 RMSD of the difference between model and climatological profiles at different depths evaluated in the 2017–2020 reference period. Statistics are computed using the 16 subregions in Fig. 3. Reference datasets for validation are: (1) EMODnet data collections (Buga et al., 2018) integrated with additional oceanographic cruises (Cossarini et al., 2015), and (2) Socat dataset (Bakker et al., 2014).

Chlorophyll from ocean colour is the most common variable used for validation and near-real-time assessment of operational biogeochemical models and allows for a comparison of the forecast skill performance among the Marine Copernicus systems. Results of surface chlorophyll skill scores show that the quality of the first day of the forecast of the BIO component is in line with that of other Copernicus models² (Spruch et al., 2020; Vandenbulcke et al., 2021; McEwan et al., 2021; McGovern et al., 2020). In particular, the two proposed accuracy indexes (i.e. one minus scatter index and one minus the root mean square error normalized on variability) of the MED model equal to 34 % and 47 %, which are within the ranges of the other Copernicus systems: 11 %–38 % and 13 %–73 % for the two skill scores respectively (Spruch et al., 2020; Vandenbulcke et al., 2021; McEwan et al., 2021; McGovern et al., 2020).

For other biogeochemical variables, a direct comparison of the accuracy among Copernicus models is not straightforward, given the different protocols for metrics computation, the representativeness of the available observations, and the large range of variability of observed values of biogeochemical variables among the European seas. Nevertheless, a rough comparative assessment of the quality of Marine Copernicus biogeochemical models can be provided using published estimated EANs normalized by the typical values of the variables (McEwan et al., 2021; Feudale et al., 2021; Spruch et al., 2020; Melsom and Yumruktepe, 2021; McGovern et al., 2020; Vandenbulcke et al., 2021) to derive a common index of relative uncertainty. As for examples, the relative uncertainty of oxygen of the MED system is on the order of 2 %, which is in line with the other Copernicus systems, except for the Baltic and Black Sea systems, which show slightly higher relative errors. For nutrients, nitrate and phosphate uncertainties of the MED are about 50 % and 35 %, which are similar to or slightly better than most of the other Copernicus marine biogeochemical systems (i.e. ranges of 30 %–75 % and 30 %–50 % for nitrate and phosphate respectively). Finally, the relative uncertainty of pH simulated by the MED system is less than 0.5 %, while other Copernicus systems report relative errors on the order of 1 %–2 %.

Beside the aforementioned comparison, the MED biogeochemical system exhibits some distinguishable features: the continuous monitoring of the forecast skill of surface chlorophyll since the beginning of the operational biogeochemical system dating back to 2010 (Salon et al., 2019), a large number of validated variables with in situ data (i.e. up to 10 variables, Table 6), and the thorough use of BGC-Argo observations for near-real-time forecast validation (Salon et al., 2019; Cossarini et al., 2021; https://medeaf.ogs.it/nrt-validation, last access: August 2022).

3.4 ECMWF forcing skill

A calibration and validation system of the ECMWF forcing fields used by the Med-MFC operational systems has been developed using in situ ground meteorological observations from METeorological Aerodrome Report (METAR) stations and numerical model data from ECMWF (see Fig. 12). Four well-established statistical indices for validating 2 m temperature, dew point temperature, air pressure, and wind speed have been defined: (a) bias, (b) RMSE, (c) Nash–Sutcliffe model efficiency coefficient, and (d) correlation coefficient.

Figure 12 Example of ECMWF wind speed validation with respect to METAR ground observations in 2019 in the area of the Gulf of Lion. (a) Time series of daily mean wind speed time series from METAR station (blue line) and from ECMWF (red line). (b) Time series of main skill metrics (bias, root mean square error (RMSE), and Nash–Sutcliffe model efficiency coefficient (nse)).

The atmospheric forcing calibration and validation system will become publicly available, and an example of this validation is provided in Fig. 12, showing daily mean wind speed time series from a METAR station (blue line) and ECMWF (red line) in the area of the Gulf of Lion during the year 2019, as well as time series of the main skill metrics.

4 Conclusions and future perspectives

In this paper, the Med-MFC components (PHY, BIO, and WAV) have been described, providing an overview of their technical specifications. The PHY component provides 3D currents, temperature, and salinity with the BIO and WAV components daily, with daily mean values. This model system is flexible enough that improvements can be carried out separately on the three components, considering the different levels of maturity of the numerical modelling parameterizations, the data assimilation components, and the validation datasets. A different data assimilation system is run for each component, making the best use of all available data from satellite and in situ observations; the effort is to assimilate as much data as possible and use background or model uncertainties to account for the missing couplings. The three components' accuracy has been evaluated for a common 3-year period, from January 2018 to December 2020.

The PHY component has been validated by comparing model data to in situ and satellite observations, showing good accuracy in representing the spatial pattern and the temporal variability of the temperature, salinity, and sea level in the Mediterranean Sea. In particular, the model has a warm surface temperature bias of +0.12 ^∘C when compared to satellite SST. The water column temperature error has a clear seasonal signal with the largest errors at the depth of the surface mixed layer and the seasonal thermocline. The model error in salinity is higher in the first layers and decreases significantly below 150 m. The SLA presents a mean average error of 3.8 cm over the 3-year-averaged period for the whole basin.

The WAV component was extensively validated for the 3-year period using all available in situ and satellite observations in the Mediterranean Sea. All statistical values calculated and presented here showed a very good system performance. It is concluded that the Mediterranean SWH is accurately simulated by the WAV component. The typical SWH difference with observations (RMSE) over the whole basin is 0.21 m (0.197 m for in situ and 0.228 m for satellite observations) with a bias ranging from −0.137 to −0.005 m when the comparison is against the in situ observations and from to −0.088 to 0.131 m when the comparison is with satellites. The scatter index (SI) exhibits low values (0.13–0.17) over the majority of the basin and relatively higher values (0.18–0.21) over the Aegean, Alboran, Ligurian, and east Levantine seas, with the highest SI value encountered in the north Adriatic Sea (0.26). As explained, the occurrence of higher SI values is mainly related to the quality of ECMWF winds in fetch-limited areas of the basin where the orographic effects play an important role and to the difficulties of wave models in appropriately resolving complicated bathymetry and coastline.

Overall, the quality of the WAV component stems from the ECMWF wind forcing that drives the wave dynamics, data assimilation, forcing from Med-PHY surface currents, and improved parameterization of the wave wind source and dissipation terms of the WAM model. In particular, the WAV component assimilates satellite altimetry data with a well-calibrated stand-alone OI scheme and implements regular updates and improved parameterization independently from the other components. Given that wind-forcing quality has a substantial influence on the model response, a considerable part of the wave product uncertainty, especially under high winds or extreme conditions, is related to the wind-forcing uncertainty and can be substantially improved by undertaking the ensemble approach in wave forecasting. The lower accuracy of the wave product in semi-enclosed regions of the Mediterranean Sea (e.g. Adriatic and Aegean seas) can be related to the current spatiotemporal resolution of the wind forcing. Near the coast, unresolved topography by the wind and wave models and fetch limitations causes the wave model performance to deteriorate. In particular, the WAV component assimilates satellite altimetry data with a well-calibrated stand-alone OI scheme and implements regular updates and improved parameterization independently from the other components. Given that wind-forcing quality has a substantial influence on the model response, a considerable part of the wave product uncertainty, especially under high winds or extreme conditions, is related to the wind-forcing uncertainty and can be substantially improved by undertaking the ensemble approach in wave forecasting. The lower accuracy of the wave product in semi-enclosed regions of the Mediterranean Sea (e.g. Adriatic and Aegean seas) can be related to the current spatiotemporal resolution of the wind forcing. Near the coast, unresolved topography by the wind and wave models and fetch limitations cause the wave model performance to deteriorate.

The BIO system has defined a validation framework (Salon et al., 2019) based on multivariate (e.g. more than 10 variables) and multilevel metrics that include Global Ocean Data Assimilation Experiment (GODAE) class-1 and class-4 statistics and process-oriented metrics. Of particular interest is that the present validation framework also includes near-real-time observations (i.e. satellite and BGC-Argo) that show average errors in the 0–200 m layer of 0.04 mg m⁻³, 0.4 mmol m⁻³, and 16.8 mmol m⁻³ for chlorophyll, nitrate, and oxygen respectively. Thus, the validation framework represents a robust benchmark for the future improvements of the Mediterranean BIO model. Indeed, as detailed in Salon et al. (2019) and Cossarini et al. (2021), critical sources of the BIO model errors include unresolved Atlantic boundary conditions, as well as land–sea and atmosphere–sea forcing uncertainty in model parameterization and the inconsistency of coupled physical–biogeochemical processes.

The value and reliability of the Med-MFC systems are demonstrated by the several downscaling coastal model systems and downstream applications that use its outputs operationally. The CYCOFOS (Cyprus Coastal Ocean Forecasting and Observing System) (Zodiatis et al., 2003), which is a subregional forecasting and observing system in the eastern Mediterranean Levantine Basin, uses the Med-MFC outputs to set its boundary conditions. The Med-MFC outputs are used as initial and lateral boundary conditions by the physical and wave ocean system MIT Operational (MITO), which provides 5 d forecasts at a resolution up to $\mathrm{1}/\mathrm{48}$^∘ (Napolitano et al., 2022). The Southern Adriatic Northern Ionian coastal Forecasting System (SANIFS) is a coastal–ocean operational system providing short-term forecasts since September 2014 (Federico et al., 2017). It is built on the unstructured-grid, finite-element, three-dimensional, hydrodynamic Shallow water HYdrodynamic Finite Element Model (SHYFEM) model and is based on a downscaling approach starting from the large-scale system Med-MFC, which provides the open-sea fields.

The CADEAU (Assimilation of national water quality data in Coastal Areas for a marine Directives oriEnted downstreAm prodUct) physical–biogeochemical forecasting system of the northern Adriatic Sea (Bruschi et al., 2021) is based on a high-resolution (up to 700 m) application of the MITgcm-BFM model (Cossarini et al., 2017), targeting water quality and eutrophication and using the daily Med-MFC products for initialization and to constrain the southern boundary.

In addition, the Gutta savinG fUel and emissions from mariTime Transport in the Adriatic region (GUTTA-VISIR) system, which can be defined as a tactical, global-optimization, single-objective, deterministic model system for ship route planning (Mannarini et al., 2016a, b; Mannarini and Carelli, 2019), uses the analysis and forecast wave and current fields from the Med-MFC in conjunction with wind fields from ECMWF.

Since 2008, the Med-MFC components have been continuously upgraded and substantially improved. The system evolution will also continue in the future following the main drivers of the three components: the Copernicus Marine Service users. Considering the PHY system, the users need finer spatial scales and higher time frequencies of the products, especially for improving the representation of the coastal scale and limited area processes in nested models, thus providing a unique opportunity to model the coastal areas at a resolution of a few hundred metres using nesting schemes as demonstrated in Federico et al. (2017), Trotta et al. (2021), and others. Users also require higher accuracy in storm surge forecasting, which can be achieved by including the explicit representation of the tidal forcing to resolve nonlinear interactions between astronomical and internal tides with the baroclinic circulation. An upgrade of lateral open boundary conditions in the Atlantic and the Black Sea would provide a better evaluation of the transport at Gibraltar and improved dynamics in the north Aegean Sea. Higher-frequency river runoff data from hydrological models, as well as more accurate salinity values at river mouths, would provide better salinity skill not only along the coastal areas but in the whole basin. Another important goal for the future is to assimilate Argo and drifter trajectory (Nilsson et al., 2010) and glider (Dobricic et al., 2010) data, as well as sea level anomalies in coastal areas. Finally, future efforts should consider ensemble forecasting to recast the deterministic forecast within a probabilistic framework, assessing the modelling uncertainties (Pinardi et al., 2011; Millif et al., 2011; Thoppil et al., 2021; Barton et al., 2021).

User needs for the future evolution of the WAV component indicate the necessity of increasing the frequency of the wave analyses, making available larger datasets, such as of the wave spectra and dedicated products (like the directional spread at peak frequency and different parts of the wave spectrum). The required increased accuracy in wave height and mean period predictions can be mainly achieved by improving the quality of the wind forcing, which is the main driving force of wave models. Bias correction of ECMWF winds and further downscaling of ECMWF forecasts are expected to improve winds and consequently wave product quality, especially in semi-enclosed areas (e.g. the Adriatic and Aegean seas) and near the coast. Assimilation upgrades with the ingestion of multi-mission significant wave heights at 5 Hz and in situ wave height measurements from high-frequency (HF) radars will improve accuracy in coastal areas of the Mediterranean Sea, while the inclusion of spectral information in the near future (e.g. Chinese-French Oceanography Satellite – CFOSAT – wave spectrum) will further improve the prediction of the sea state. Finally, the development of a WAV ensemble prediction based on ECMWF operational ensemble winds is expected to improve the existing accuracy of the deterministic forecast at lead times beyond 48 h, providing, in parallel, uncertainty estimates of wave parameters.

User requirements for the BIO component developments include improved quality and products tailored to ecosystem and coastal applications. The validation results have contributed to identifying ameliorable model process representations and model parameter estimates that can be improved. These include a better representation of vertical nutrient and plankton dynamics and a greater number of phytoplankton functional types and zooplankton compartments to describe the diversity of the plankton community and the different energy and matter pathways in the ecosystem. In addition, the integration of optics and biogeochemistry, including novel hyperspectral and high-resolution radiometric data, can be used to better represent photosynthesis and light-related processes and to calibrate parameters of important ecosystem processes (Lazzari et al., 2021). Assimilation of new in situ profile sensors and variables (e.g. BGC-Argo float and glider) will help increase the reliability of BIO products, especially along the water column (Cossarini et al., 2019). Higher-quality vertical dynamics can be achieved through a better representation of vertical model error covariances by means of ensemble (Carrassi et al., 2018) or joint physical–biogeochemical data assimilation techniques. Finally, revising nutrient and carbon inputs from rivers (e.g. from monthly climatologies to daily observations or model predictions) will allow better resolution of coastal dynamics and coastal–offshore patterns in critical areas.

A companion paper (Part 2), which is currently being prepared, will demonstrate the capacities of the Med-MFC components in describing the medicane effects on the ocean. In particular, the Med-MFC physics, biogeochemistry, and wave components will be used to describe the effects of Medicane Zorbas (27–30 September 2018) on the ocean variables.

Appendix A

Table A1 List of the NEMO and WW3 numerical setups for the PHY component.

Table A2 List of the WAM model setup for the WAV component.

Table A3 List of the OGSTM-BFM model setup for the BIO component.

Table A4 River sources implemented as freshwater inputs in the physical and biogeochemical models, including river name, the annual mean runoff, and the imposed salinity at the river mouth.