**1. Introduction**

Waves constitute the interface between ocean and atmosphere and have an important role in terms of exchanges through this interface [1]. Their representation is necessary to accurately compute the different air–sea fluxes of heat and momentum [2].

There is a widespread worldwide offer of accurate and reliable wave forecast services. A variety of operational wave forecast services, ranging from global to local coastal scales, are run by different operational oceanographic centres (some of them national weather offices); the wave forecast products benefit different end-users, supporting day-to-day operations at sea and contributing to warning systems that minimize potential risks for marine safety (among others). The authors in [3], in their review of European Operational oceanographic capacities, indicated how several wave models (i.e., WAM (Wave Model), SWaN (Simulation Waves Nearshore), WaveWatch-iii, WWM-II, etc.) [4–6] are used in the forecast services delivered by the operational oceanographic centres. Some of these operational services use operational assimilation schemes to account for near real-time observational wave information, especially from satellite altimeters [7].

**Citation:** Toledano, C.; Ghantous, M.; Lorente, P.; Dalphinet, A.; Aouf, L.; Sotillo, M.G. Impacts of an Altimetric Wave Data Assimilation Scheme and Currents-Wave Coupling in an Operational Wave System: The New Copernicus Marine IBI Wave Forecast Service. *J. Mar. Sci. Eng.* **2022**, *10*, 457. https://doi.org/10.3390/ jmse10040457

 Academic Editor: Liliana Rusu

Received: 11 February 2022 Accepted: 22 March 2022 Published: 24 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Upper ocean dynamics are strongly affected by sea-state dependent processes, inducing the impact of waves on the ocean's small- and large-scale circulation.

On the one hand, waves affect the ocean surface layer through different processes [8]: Waves induce surface current via the Stokes drift, adding a term on the Coriolis effect in the momentum equation (the so-called Stokes–Coriolis force). Part of the atmospheric wind stress contributes to the wind-wave growth, thus, subtracting a quantity of energy to the ocean currents. Furthermore, during wave breaking, turbulent kinetic energy is produced and affects the upper ocean surface layer, enhancing the turbulent mixing. Recent studies have attempted to determine the impacts of wave effects on the representation of the ocean surface layer at different spatio-temporal scales. Among others: wave-induced mix-layer depths representation [9], relevant impacts on the atmospheric surface temperature, pressure, and precipitation [10,11], modifications in wind stress by the rise of roughness length and friction velocity [12]. This is especially true during storm events, when wave–current interactions might represent a leading order process of the upper ocean. In this context, ref. [13] strongly recommends using an ocean-waves-atmosphere coupled system to improve the representation of tropical cyclones' intensity, structure and motion. Indeed, ref. [14] studied the effect of sea waves on the typhoon Imodu (15–19 July 2003). Moreover, ref. [15] demonstrated how a coupled system simulates more accurate surface dynamics than uncoupled models, with larger improvement on the shelf, showing that (especially during extreme events) ocean-wave coupling improves the accuracy of the surface dynamics, with larger improvements in the simulation of ocean currents over the shelf due to the synergy between strong tidal currents and more mature decaying waves.

On the other hand, the presence of ocean currents affects the waves, changing their amplitude, frequency and direction. This is generally due to the energy bunching, accounted in the wave energy balance when the velocity of the wave energy propagates across the current, the energy transfer between waves and currents, the frequency shifting (including Doppler shifting) and current-induced refraction [16]. Ref. [17] accounts for significant wave height changes in the Baltic Sea due to the impact of ocean currents (up to 20% in specific severe storm conditions, mostly in shallower waters and when waves and surface currents propagate in opposite directions [18,19]). The Copernicus Marine Service [20,21], one of the streamlined six thematic streams of the Copernicus Services (Atmosphere, Marine, Land, Climate Change, Security and Emergency) [22,23] and internationally recognized as one of the most advanced service capabilities in terms of ocean monitoring and forecasting, provides regular systematic reference information on the physical, biogeochemical and sea-ice state for the European regional seas and the global ocean. This service recently included, in its product portfolio, essential ocean variables related to the sea state, and near-real-time wave forecasts, and multi-year wave reanalysis products were progressively incorporated (along the 2015–2018 development phase) in the Copernicus Marine Service offer. The Copernicus Marine Service high-level strategy includes a roadmap with associated Research and Development (R&D) priorities [24], which identifies some developments per thematic area that are key for the future service evolution. Among others, (i) upgrade of data assimilation schemes (to improve the analysis and reanalysis capabilities) and (ii) enhancement of the representations of coupling effects between ocean-wave-sea-ice-atmosphere-land components (to improve forecast model solutions) are seen as prioritized research lines for any Copernicus Marine Monitoring and Forecasting Centres (MFC).

Specifically, for the European Atlantic Façade, the Copernicus Marine IBI-MFC (Iberia– Biscay–Ireland Monitoring Forecasting Centre) delivers daily ocean model estimates and forecasts of different physical and biogeochemical parameters, including, since 2016, hourly wave forecasts and multi-year products [25].

The present work focuses on the research performed to develop the current operational version of the IBI-MFC wave model application. This research was mainly conducted to improve the accuracy of these IBI-MFC wave model products, by means of developing a new coupled ocean-wave modelling framework that also includes wave data assimilation. In that sense, this study has two specific objectives:


To address these questions, different wave model sensitivity tests are performed. Several wave model simulations generated with the IBI-MFC wave model set-up, and only differing from the operational version (available in 2018) for the activation of the new data assimilation scheme and in the degree of the ocean current forcing applied, are run. The assessment of these model simulations is conducted using several local available in-situ and satellite wave observations.

The paper is organised as follows: Section 2 provides a description of the Copernicus IBI-MFC wave model system and the different model sensitivity tests performed, together with the model assessment proposed. Section 3 presents the main results, with an analysis of the proposed updates. Finally, in Section 4, the impacts of the current forcing interactions and the data assimilation scheme proposed are discussed, providing a look ahead to related benefits on the IBI wave operational forecast capabilities.

#### **2. Methodology and Sensitivity Tests for Copernicus Marine IBI-MFC Wave System** *2.1.TheIBIAreaandIBI-MFCWaveModel*

The Copernicus Marine IBI-MFC (Iberia–Biscay–Ireland Monitoring Forecasting Centre) o ffers a comprehensive portfolio of regular and systematic regional information on the state of the ocean for the European Atlantic façade, supporting all kinds of marine applications. As part of this IBI-MFC service, a short term (10-days) high-resolution wave forecast is updated twice a day for the IBI area. Hourly instantaneous data for significant wave height, wave direction, wave period variables, together with wind sea and swell (primary and secondary partitioned wave spectra) parameters are delivered as part of this regional Copernicus Marine IBI wave product.

The MFWAM model configuration for IBI MFC is implemented on the IBI domain (26–56◦ N and −19–5◦ E; see geographical domain in Figure 1) with a horizontal resolution of 5 km approximately (1/20◦).

The wave model used as base of this IBI-MFC operational system is the MFWAM [26]. This MFWAM model is based on the IFS-ECWAM 41R2 code [27], with changes regarding the dissipation by wave breaking and the swell damping source terms as developed by [28]. The current version of the model includes major improvements achieved within the FP7 European research MyWave Project [2,29]. The IBI-wave model performs a partitioning technique on wave spectra over all ocean grid points of the IBI domain. The partitioning technique is based on the watershed method developed for image processing [30]. This process effectively treats the wave spectrum as a topographic map from which individual peaks in wave energy can be identified to define the separate wave components. First, wave spectrum is split in wind sea and swell wave spectra. Then, partitioning is applied for the swell wave spectrum. The peaks on the spectrum are isolated and they are considered as partitions. Afterward, classification of swell partitions in primary and secondary swell is performed depending on the mean energy of each partition.

**Figure 1.** The CMEMS IBI-MFC wave Forecast/Analysis model Application. Model and product service spatial coverages, and details on the forcing, Open Boundary conditions and external data sources used for the ocean current forcing and the Data Assimilation applied.

The IBI-MFC wave model was upgraded (January 2018 Operational release) to improve the drag coefficient variation with the wind speed, resulting in positive impacts in the surface stress characterization. This improvement of the surface stress is certainly needed for the coupling with the IBI ocean currents here tested. To this end, a new setting on the wave dissipation term, the sheltering parameter, and the use of Phillips spectrum tail for the high-frequency part of the wave spectrum was also implemented. Moreover, the minimum water depth is taken as 5 m (instead of the 1 m value used in former IBI wave model versions). Associated to these upgrades, slight improvements in terms of significant wave height were obtained (reduced scatter index, around 1.9% when comparing model results with observations from altimeters [31]).

The bathymetry used is derived from the 1 arc minute ETOPO1 ocean bathymetry by National Geophysical Data (NOAA) [32]. The wave spectrum is discretized in 24 directions and 30 frequencies ranging from 0.035 Hz to 0.56 Hz. The Copernicus Marine IBI wave model is driven by short-range forecasted and analyzed IFS-ECMWF winds [33] at 1/8º hourly winds, which are used for the first 90 h, decreasing time frequency to 3-hourly until T + 144 h, and 6-hourly forecasts until T + 240 h. It uses as boundary conditions (wave spectra) Copernicus Marine GLOBAL wave data at 1/10º spatial resolution [34].

An IBI-MFC wave forecast product (10 days hourly data updated twice a day) is delivered to end-users for the IBI service domain (see spatial coverage in Figure 1), together with a 2-year historic timeseries, composed of IBI best estimates (analysis data for the Day-1 date).

#### *2.2. Model Sensitivity Test: The Proposed IBI Wave Model Upgrades*

The Copernicus Marine IBI-MFC identified, among others, the following two major shortcomings in its wave forecast service:


Bridging these 2 identified gaps in the Copernicus IBI wave service was considered as a major goal for the IBI-MFC service evolution planned for the last Copernicus-1 Service Phase (2018–2021) and 2 different research working lines were followed to achieve the objective. These specific research lines were fully aligned with the Copernicus Marine Evolution Strategy [20], and their scientific research priorities (implementation of data assimilation schemes and enhancement of the model coupling between different earth system components) were two of the major amelioration axes proposed for the Copernicus Marine products and services.

The present work aims to quantify the potential added value of new IBI wave analyses (to be generated by means of a new data assimilation scheme implemented to assimilate altimetric significant wave height satellite observations) with respect to the IBI best model wave estimates (note that traditionally, the IBI-MFC was delivered as historic best model estimate 2 years of model hindcast wave data, wave run for day D–1 forced with analyzed winds). Likewise, the impact of including the contribution of the surface ocean currents on the IBI wave model solution is evaluated.

The impacts on the IBI wave solution of both the new IBI wave DA system and the current–wave coupling are evaluated through specific IBI-like wave model scenarios. To this aim, different model sensitivity test runs, based on the IBI-MFC operational model configuration, have been designed and run. Table 1 shows an overview of the four proposed IBI wave model scenarios: the Control run (IBI-CO) was performed using the same wave model set-up that was used in the IBI operations in December 2019. Two more IBI wave sensitivity test runs were performed: one with the current forcing activated, but without Data Assimilation (the IBI-CU run) and another with data assimilation activated, but without any coupling contribution (the IBI-DA run). Finally, a last model scenario that includes the two novelties (data assimilation and ocean current forcing) was assessed. This model scenario (hereafter named as IBI-OP) is consistent with the new IBI-MFC operational CMEMS IBI wave model system (in operations since July 2020).

**Table 1.** IBI wave model scenarios: overview of model runs performed to test the 2 proposed novelties (the current forcing and the new Data Assimilation scheme). The Control run and all the test runs using the IBI wave model set-up used in operations in Copernicus release in December 2019. The last sensitivity test, switching both the DA and the currents coupling, to be proposed as base of the new Copernicus IBI wave forecast system (in operations since July 2020).


The four different model scenarios proposed were run over the year 2018. The IBI wave model set-up novelties were assessed, and all the proposed sensitivity model runs were validated using in-situ and remote sensing observations (i.e., from coastal and deep-water mooring buoys, and altimetric SWH). A complete description of the different sensitivity model tests performed to improve the operational wave model application that comprise the IBI-MFC operational wave forecast service is provided below. Results from the different model test runs are provided in the following section.

#### 2.2.1. The Altimetric Wave DA Scheme Proposed for IBI

The data assimilation scheme proposed to be applied in the IBI wave service and tested through the IBI-DA test run presented here, is based on the optimal interpolation scheme described by [35] and it is the same scheme used in the Copernicus Marine GLOBAL wave system. The variable to be assimilated is the significant wave height, Hs. Because Hs is not a state variable of the system, this introduces an extra complication in that the energy must be repartitioned from a frequency and direction integrated parameter (the Hs) to the full directional frequency spectrum. This involves making several assumptions and is by nature inexact, but in practice performs well [7]. What follows is a brief description of the method, which has been adjusted for the ST4 physics used in the IBI-wave model.

For a state vector x, optimal interpolation seeks a weighted combination of the background (or mode forecast), denoted by *xf*, with observations, *y<sup>o</sup>*, in order to produce an analysis *<sup>x</sup>a*. The fundamental equation is:

$$\mathbf{x}^{a} = \mathbf{x}^{f} + \mathbf{K} \left(\mathbf{y}^{\rho} - \mathbf{H} \mathbf{x}^{f}\right), \tag{1}$$

where **H** is the observation matrix. **K** is a weighting matrix.

$$\mathbf{K} = \mathbf{P} \mathbf{H}^{\mathrm{T}} \left( \mathbf{H} \mathbf{P} \mathbf{H}^{\mathrm{T}} + \mathbf{R} \right)^{-1},\tag{2}$$

where the matrices **P** and **R** are respectively the model and observation error covariance matrices. In MFWAM these matrices are expressed as correlation matrices:

$$\mathbf{K} = \mathbf{C} \mathbf{H}^{\mathrm{T}} \left( \mathbf{H} \mathbf{C} \mathbf{H}^{\mathrm{T}} + \mathbf{I} \right)^{-1},\tag{3}$$

The ratio of background and observation errors is kept constant over the IBI domain and set to 1 (i.e., model and observation errors are assumed to be equal everywhere). Although a different ratio may be theoretically justifiable, we have found that this value works best for this model in this domain; it is the same as that used in the global configuration. Observational errors are additionally assumed to be uncorrelated. With these assumptions, R is none other than the identity matrix, I. P becomes the correlation matrix *C* defined in terms of the correlation length *λc*:

$$x\_{i\bar{j}} = e^{-\left(\frac{d\_{i\bar{j}}}{\lambda\_{\bar{\varepsilon}}}\right)^{\mu}},\tag{4}$$

where *dij* is the grea<sup>t</sup> circle distance between points *i* and *j* and a is a tuning parameter.

With this simplified OI scheme the only tunable parameters are in the correlation function. For the IBI configuration used here the correlation length *λc* is set to 170 km, significantly less than the 300 km length used in the global configuration, and the tuning parameter a is set to 1. We performed some experiments with *λc* and *a*, in particular testing values determined from a correlation study of the global model divided by basin [36]. These did not result in an improvement in performance, however, so the original values were kept. The cutoff distance, beyond which observations are not included in the analysis, is set to 650 km. The analysis, Equation (1), produces a corrected estimate for the significant wave height. Since *Hs* is not a state vector of the wave model, but rather an integrated parameter, in order to correct the model itself MFWAM redistributes the energy in the wave spectrum using the method of [35], which is based on empirical wave growth laws. The analyzed spectrum is expressed as:

$$F^a = ABF^f(Bf, \theta)\_\prime \tag{5}$$

where *F* denotes the wave spectrum, *f* the frequency and *θ* the direction, and the superscripts *a* and *f* refer to analysis and background respectively. *A* is the ratio of analysis to background energy, which can be expressed as (*Hsa*/*Hsf*)2, where *A* determines the overall energy correction to the spectrum, while the effect of *B* is to rescale the frequencies.

Two different methods for computing *B* are used, depending on whether the spectrum is determined to be primarily a swell spectrum (the energy of the swell accounts for more than 1/4 of the total energy of the spectrum) or a wind–sea spectrum. If the spectrum is predominantly wind–sea, *B* is computed from the mean frequencies of the background and analysis as:

> *B*

$$
\overline{f} = \overline{f}^f / \overline{f}^a \,. \tag{6}
$$

The choice of mean frequency was for computation reasons; the peak frequency would be just as valid a choice (if not more so), as both are in any case approximations. If the spectrum is predominantly swell, the average steepness of the waves is assumed constant. Therefore:

$$B = \overline{f}^f / \overline{f}^a = \sqrt{\mathcal{H}\_s^a / \mathcal{H}\_s^f} \tag{7}$$

The calculation of *f f* and *Hsf* is taken care of by the model, but fa and *Hsa* have to be estimated from the *Hs* analysis. This is done by exploiting the empirically derived duration-limited growth laws. By defining the non-dimensional energy, mean frequency and duration, respectively:

$$
\varepsilon \, = \mathfrak{u}\_\*^4 \varepsilon / \mathfrak{g}^2,\tag{8}
$$

$$f' = \mathfrak{u}\_\* f / \mathfrak{g}\_\* \tag{9}$$

$$\mathbf{f}' = \mathbf{t} \mathbf{g} / \mathbf{u}\_{\*},\tag{10}$$

where *g* is gravitational acceleration, and *= Hs2/4*. The growth relations are:

$$
\epsilon^{\prime}(t) = 1877 \left( t^{\prime} \left[ t^{\prime} + 5.440 \ast 10^{5} \right] \right)^{1.9}, \tag{11}
$$

$$
\epsilon \left( \overline{f} \right) = 5.054 \ast 10^{-4} \overline{f}^{-2.959} \text{ \AA} \tag{12}
$$

From the background friction velocity and *Hs* we can use Equation (9) and the growth law equations to estimate an updated ' and *f*'. These in turn give us *f a* and *Hsa*, and with these we can calculate *B* and *A* and produce the updated spectrum (5).

The data are combined to generate super-observations. Within a set time window around the analysis time, the data are assigned to model grid points closest to them, and any data sharing a model point are averaged together and treated as one observation. This reduces the number of data to assimilate, lightening the computational load, and it helps to smooth out potential errors in the observations. Outlying data are rejected in this step.

This data assimilation method results in a final *Hs* corrected mostly in the wind sea part of the frequency spectrum [7].

#### 2.2.2. Wave–Current Coupling Proposed for IBI

To incorporate surface ocean currents forcing in the IBI wave model system, surface current data from the CMEMS-IBI analysis and forecast ocean model system were used as inputs in the IBI wave model system for wave refraction.

The presence of current may change the amplitude, frequency, and direction of waves. This is generally due to the energy bunching that is readily accounted for the energy balance equation of waves using the velocity of the wave energy propagating across the current, the energy transfer between waves and currents, the frequency shifting (including Doppler shifting) and current-induced refraction [16].

MFWAM model equations include the depth and current refraction. The propagation velocity in the relative frequency space should be computed according to [37]:

$$\mathcal{L}\_{\sigma} = \frac{\partial \sigma}{\partial d} \left[ \frac{\partial d}{\partial t} + \stackrel{\rightarrow}{\boldsymbol{\mu}} \cdot \stackrel{\rightarrow}{\nabla} d \right] - \mathcal{c}\_{\mathfrak{F}} \stackrel{\rightarrow}{k} \cdot \frac{\partial \stackrel{\rightarrow}{\boldsymbol{u}}}{\partial t} \,, \tag{13}$$

where σ is the relative frequency, *cσ* the propagation velocity in the relative frequency space, *t* is time coordinate, *s* is the space coordinate in the direction of propagation, → ∇ the gradient operator in the geographical space, *d* is the total depth, →*u* is the current, *cg* is the group velocity and → *k* is the wave number vector. As in MFWAM model current and water depth are time-independent, the term *∂d*/*∂t* in Equation (13) is not present.

The offline method of coupling with surface currents takes the files needed for the whole forecast period from the IBI-MFC ocean forecast system. The files used are not exactly the ones delivered through the Copernicus Marine catalogue, but rather the native IBI NEMO model outputs (at the 1/36◦ ORCA grid), which the IBI wave system interpolates from the IBI ocean model grid into the 1/20◦ regular grid used for the MFWAM model.

#### 2.2.3. Assessment of Model Runs: Evaluation Criteria against In-Situ and Altimeter Data

To assess the performance of the numerical model applying both novelties and to identify the main sources of uncertainty linked to the Wave–Current coupling and the application of the data assimilation in the IBI wave model simulations, the four different model scenarios were performed over the year 2018. The significant wave height and mean period fields, resulting from the different IBI sensitivity runs, are validated by means of comparison with different in-situ and satellite remoted sensed observational data sources.

In-situ measurements of significant wave height, Hs, and mean wave period, *Tm*02, were extracted from mooring buoys available in the IBI region, compiled in the product delivered by the Copernicus Marine IBI INSITU-TAC [38]. The mean wave period (*Tm*02) is defined as follows:

$$Tm\_{02} = 2\pi \left( \frac{\iint \omega^2 E(\omega, \theta) d\omega d\theta}{\iint E(\omega, \theta) d\omega d\theta} \right)^{\frac{-1}{2}},\tag{14}$$

where *E*(*<sup>ω</sup>*,*θ*) is the variance density and w the absolute radian frequency.

Measurements for *Hs* and *Tm*02 variables are available for the examined year in the area at 49 and 45 buoys, respectively (see list and locations in Figure 2).

To validate the model through model output and buoy data collocation, the time series were taken at the model grid point nearest to the buoy location. For regional validation purposes, the IBI domain is split into different sub-regions of interest, being validation metrics gathered for the whole IBI service domain and for each sub-region (see spatial domains in Figure 2). Likewise, different metrics are computed separately using coastal and deep-water mooring buoys.

The different model sensitivity test runs were validated with satellite observations of significant wave height, Hs. However, since data from Jason-2, Jason-3, Saral, Cryosat-2 and Sentinel3 altimeters are now assimilated into the model (information from this mission included in the L3 CMEMS altimeter data product used for assimilation), an independent data source is needed. Thus, the diagnostic after data assimilation is performed by comparing the model to the HY-2A satellite altimeter processed by the French Space Agency CNES. The validation procedure with altimetric observations begins with pre-processing the *Hs* data, rejecting *Hs* data lower than 0.5 m or higher than 12 m and eliminating big jumps in terms of *Hs* value from one observation to the next one (the biggest value of steps higher than 1 m in case of positive steps and −2 m, in case of negative steps, is rejected). To validate the model runs with altimetric observations we used the HY-2A satellite data for both DA and reference runs, ensuring that a data source independent from the assimilated data, HY-2A *Hs* is biased [39], so a calibration to reference mission such as Jason-3 has been implemented on crossover locations. This leads to a linear for HY-2A Hs, which is expressed as follows:

$$H\_s = 0.9476 \times H\_s^{biased} - 0.0230,\tag{15}$$

The modelled *Hs* is also post-processed with an upscaling to a 0.1 degree resolution (the nearest grid point for altimeters), in order to closer match the observations, and it is limited to values above 0.5 m. The final step to prepare for validation is then to find the points of modelled *Hs* that correspond to the observation points. The validation statistics are then computed using these two values for each point where a valid observation exists.

Apart from bias, root mean square difference (RMSD) and correlation (CCOR), validation against altimeters includes the statistical quantity scatter index (*SI*2) used for the wave model statistics defined as:

$$SI2 = \frac{\sqrt{\sum\_{i=1}^{N} \left[ (\mathbf{x}\_i - \overline{\mathbf{x}}) - (\mathbf{y}\_i - \overline{\mathbf{y}}) \right]^2}}{\sum\_{i=1}^{N} \mathbf{y}\_i},\tag{16}$$

where *yi* is the observation, *xi* the model value corresponding to the *i*th observation, *N* is the total number of observations, and the overbars refer to the population mean. This definition of scatter index differs from others in that the observations are not squared before taking the mean, so is only valid for quantities such as *Hs* which are always positive.

**Figure 2.** Locations of all the mooring buoys used in the model validation performed for the year 2018. In the accompanying table, the seven-digit WMO identifier (or buoy name provided by CMEMS) is followed by the information on the location of the mooring buoy (coastal or deep water). The 8 different sub-regions of interest for validation purposes in which the IBI service domain is divided are shown in the map (see polygons of different colors). From the list of 60 buoys in the IBI region, only the 49 available for year 2018 have been used.

A lower scatter index is not always a reliable gauge for model performance [40]. The symmetrically normalized root square error (*HH*), also used for validation against altimeters, provides more accurate information on the accuracy of simulation.

The error indicator *HH* proposed by Hanna and Henold, ref. [41], is defined as the RMSE divided by the absolute value of the mean of the product of the observations and modelled values:

$$HH = \sqrt{\frac{\sum\_{i=1}^{N} (x\_i - y\_i)^2}{\sum\_{i=1}^{N} y\_i x\_i}},\tag{17}$$

Then, in the next section, the impact for the two proposed novelties (current forcing and data assimilation) is assessed with tests performed in 2018, following the evaluation criteria against in-situ and altimeter data described above.
