**1. Introduction**

In the global water cycle, ocean-atmosphere surface net freshwater fluxes are balanced by ocean transport and land discharges within the oceans. The return of terrestrial runoff into oceans is mostly concentrated at the mouths of the world's major rivers, and it locally provides significant freshwater inflows, forcing changes in sea water densities. Consequently, realistic estimates of land freshwater discharges into oceans are needed to globally study water budgets [1]. At regional and local scales, and very specifically on coastal and shelf areas, the river runoff contribution is a major forcing of the ocean dynamics, having a strong influence on water stratification [2]. Besides, it also introduces significant fluctuations in local circulation patterns, generating important buoyant freshwater river

**Citation:** Sotillo, M.G.; Campuzano, F.; Guihou, K.; Lorente, P.; Olmedo, E.; Matulka, A.; Santos, F.; Amo-Baladrón, M.A.; Novellino, A. River Freshwater Contribution in Operational Ocean Models along the European Atlantic Façade: Impact of a New River Discharge Forcing Data on the CMEMS IBI Regional Model Solution. *J. Mar. Sci. Eng.* **2021**, *9*, 401. https://doi.org/10.3390/jmse9040401

Academic Editor: Alberto Ribotti

Received: 9 March 2021 Accepted: 1 April 2021 Published: 9 April 2021

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

**Copyright:** © 2021 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/).

plume-like flow structures [3]. Thus, noticeable on-shelf baroclinic frontal areas are caused and maintained in time by the freshwater contributions of one or various major rivers discharging in nearby coastal areas. These areas, generally located in coastal shelf seas or into estuaries where currents patterns are governed by density differences between salt sea water and fresh river water are known as Region of Fresh Water Influence (ROFI, [4]). The discharged freshwater volume, and its temporal variability, impact on the local circulation of such ROFI areas, by means of different interactions between the linked buoyancy plume structure and any existing regional ambient flow features (i.e., major current flowing along the slope, currents linked to open sea dynamics, barotropic tidal shelf currents, or ocean transports responding to prevailing wind regimes [5,6]). The coastal freshwater inputs even modulate typical coastal effects, such as upwelling events [7], having major impacts on the productivity and the whole ecological environment [8].

Due to this key role played at different scales on the ocean system, there is a growing interest in including realistic freshwater signals coming from land in operational ocean model systems.

In global ocean and basin models, which are extensively used in climate prediction and reanalysis production, this freshwater input is crucial to maintain water balances and to reproduce salinity fields with realistic latitudinal gradients that affect climate signals such as the Atlantic Meridional Overturning transports [9]. It is important to note that these global ocean model applications are not fit-for-purpose for performing on coastal areas as many processes of paramount importance to reproduce shelf dynamics are parameterized, filtered out, or simply not included. Furthermore, due to the typical global model coarse resolution, coastlines are not realistically represented. This limitation of global ocean models to deal with coastal areas has conditioned the way to force such models with punctual freshwater river discharges. Since most global models cannot accurately represent river mouths or estuaries, the freshwater river contribution has been traditionally included by means of some modification in surface ocean-atmosphere water fluxes along coastal areas (mostly modifying the evaporation-precipitation (E-P) ratio along coastal model grid points, as proposed by Bourdallé-Badie and Tréguier [10]).

The use of coastal and river freshwater forcing derived from climatological data has been a primary, and still a very commonly accepted, approach in global ocean modeling. However, available databases and services that progressively provide better estimates of coastal runoff and river discharges have started to improve ocean model freshwater coastal forcing input. Several research initiatives and established services are focused on generating at a global scale up-to-date runoff and freshwater river discharge historical datasets. Dai et al. [11] compiled for community use a global dataset with monthly streamflow at farthest downstream stations for the world's 925 largest ocean-reaching rivers; likewise, the service currently provided by the Global Runoff Data Centre [12], under the auspices of the World Meteorological Organization (WMO), delivers quality controlled "historical" databases, comprising daily and monthly mean river discharges from more than 9900 in-situ stations from 159 countries. The GloFAS-ERA5 (Global Flood Awareness System) operational global river discharge reanalysis [13] produces, consistently with the ECMWF ERA5 (European Centre for Medium-Range Weather Forecasts global reanalysis, version 5th) atmospheric reanalysis and using a river network, a global gridded data product (with a horizontal resolution of 0.1 degree and at a daily time step) that covers from 1979 to near-real-time (within a delay of 7 days).

However, when we move into coastal ocean modeling, approaches and needs are slightly different: in their review of coastal modeling science foundation, Kourafalou et al. [14] pointed out the need of applying very high-resolution models at coastal scales, able to accurately capture coastal features, and forced with realistic punctual river inflows to adequately reproduce the related dynamics occurring in ROFI areas. This demand in terms of high-resolution models is crucial when forecasting estuarine circulation. Model applications focused on estuaries show horizontal resolutions that range from kilometric scales (at coastal boundaries) to meters (within the estuary zones). These extremely high coastal model

resolutions should be accompanied by a combination of ad hoc riverine and atmospheric forcing with appropriate temporal variability (being needed daily, and even higher temporal frequency data). Therefore, neither river discharge monthly climatological products nor coarse resolution watershed model estimations are fully satisfactory when used to force operational ocean coastal models [15].

Lately, the utility of dynamical model downscaling from the deep ocean, across the continental shelf and into the coastal areas, has been extensively demonstrated [16–18]. In the framework of the Copernicus Marine Environment Monitoring Service (CMEMS), a global ocean model together with a wealth of regional ocean model systems, are operationally running to deliver short-term forecast and multi-year reanalysis information on the physical and biogeochemical ocean for the global ocean and the European regional seas (CMEMS, 2020 [19]). Specifically, along the European Atlantic façade, CMEMS provides a regional service (besides its global one) for the so-called IBI area (Iberia-Biscay-Ireland regional seas [20]). These CMEMS IBI regional products, down-scaled from the CMEMS global ones, are aimed to characterize upper-ocean dynamics at finer scales and to provide better boundary conditions for very-high-resolution coastal models, embedded into IBI solution, used by other non-CMEMS downstream services (as pointed out by Le Traon et al. [21], among others, one of the goals of the CMEMS mission).

The latest review of the current European capacity in terms of operational marine and coastal modeling systems [22] maps the organizations and the operational model systems at regional and coastal scales across Europe and points out the sustained availability of CMEMS global and regional scale core products as a positive driver to favor proliferation of "downstream" services devoted to coastal monitoring and forecasting. Likewise, this review highlights the availability of different operational ocean models, covering in some cases overlapped areas.

The availability of different operational regional and coastal model solutions in some IBI areas is seen as an opportunity for on-going initiatives, such as the MyCoast Project [23]. This EU Interreg Atlantic Program project is aimed at enhancing the regional cooperation, taking benefit of synergies among existing operational oceanographic tools and services to develop improved coastal risk services to support marine safety, fight against pollution, and generally, to increase preparedness to extreme met-ocean coastal events. MyCoast develops and upgrades different risk tools, making them suitable for using different operational model solutions as inputs. In parallel to this interoperability enhancement, a multi-model intercomparison exercise was conducted to deepen the knowledge on the different model products used as forcing by the MyCoast risk tools, identifying differences and analogies between 9 operational regional and coastal model solutions on overlapping areas. This multi-model comparison exercise, operationally updated since 2018, shows general consistency among the models, especially for sea surface temperature (SST). However, noticeable differences in sea surface salinity (SSS), mostly on shelf areas, are found (see Figure 1).

These modelled salinity differences are particularly evident in ROFI zones, such as the northwestern Iberian shelf and the French shelf in the Gulf of Biscay, and seem mostly related to the diverse treatment that each operational ocean model system gives to the coastal and river freshwater discharge forcing. Indeed, Matulka et al. [24] inform that the approaches followed to include freshwater inputs in each operational forecast system, shown in Figure 1, vary from one model set-up to another, using: (1) data derived from daily observations taken in river flow stations (i.e., the MeteoGalicia MTG-ROMS set-up and the CMEMS IBI one, but this one using observations only for some of the major rivers), (2) daily freshwater discharges at river mouths produced by hydrological models (this is the case for IBI, especially for its five-days-ahead forecast run), (3) data derived from river discharge monthly climatology (a very common approach, fully or partially used in IBI, IMI-ROMS and SHOM-HyCOM operational ocean model set-ups), and finally, (4) the no river forcing option (in this case, used only by the IST-MOHID operational version, available when the study was performed). Independently of the approach used to include

this freshwater signal, what seems clear is that having a non-optimal representation of this coastal and river freshwater forcing in ocean model set-ups can lead to biased simulated ocean states, jeopardizing the reliability of ocean forecasts, and diminishing forecast skill on ROFI areas.

**Figure 1.** 2018 Spring (March-April-May) Sea Surface Salinity (SSS) field (psu) simulated by different operational ocean models in the IBI region. The following ocean model solutions, included in the MyCoast Multi-Model Exercise, are shown (each application based on a different model and covering a different geographical domain): The CMEMS IBI regional solution, based on NEMO model and covering the whole IBI region (CMEMS-IBI-NEMO), the *Puertos del Estado* model (PdE-MITGCM) I for the Strait of Gibraltar (dark blue), the *Instituto Superior Tecnico de Lisboa* (IST-MOHID) system for the western Iberian domain (purple), the *MeteoGalicia* (MTG-ROMS) forecast service for the Northwestern Iberian domain (cyan), the Irish Marine Institute (IMI-ROMS) application for Ireland (red), the *Service hydrographique et océanographique de la Marine* application (SHOM-HyCOM) for the Gulf of Biscay (yellow).

The MyCoast multi-model intercomparison exercise is further proof of the necessity to improve coastal freshwater signals in operational ocean models. This need, already identified within the operational oceanographic community [22], includes improved standardized freshwater river inputs (considering not only flow rates but also particulate and dissolved matter key for biogeochemical processes) and foresees, as a longer-term objective, to include connection and coupling with land hydrology models.

Likewise, the CMEMS Service, as part of its Evolution Strategy and scientific research priorities (CMEMS Scientific and Technical Advisory Committee, 2017 [25]), identified the enhancement of land forcing within their medium- and long-term objectives.

CMEMS has supported in the last years different research projects on cutting-edge science and technological developments, needed to ensure future state-of-the-art service evolutions. Within the last call for tenders (2018–2020), two CMEMS service evolution projects were dedicated to make the required research to improve the freshwater river discharge inputs for global and regional CMEMS models: the BRONCO (Benefits of dynamically modelled river discharge input for ocean an coupled atmosphere-land-ocena systems) Project [26] was aimed at generating a global river discharge reanalysis, using the GloFAS system driven by the ECMWF ERA5 atmospheric reanalysis. On the other hand, the LAMBDA (Land Marine Boundary Development and Analysis) Project [27], regionally focused on the European Atlantic Façade and the North Sea, had the objective of achieving better estimations, by means of high-resolution hydrological model simulations, for characterizing major, and also some minor, river contributions and their day-to-day variability (highly needed for operational ocean coastal models). The resulting new freshwater model estimates, more realistic and adjusted using observations along historical periods, reproduce not only more reliable long-term changes in continental freshwater discharges into the oceans, but also a more detailed spatial distribution. These model products and in-situ observational data are operationally updated, made available and downloadable into the LAMBDA product viewer web interface (http://www.cmems-lambda.eu/mapviewer/ (accessed on 8 April 2021)).

The present work aims to quantify the potential added value of using this new LAMBDA freshwater discharge product as part of the fresh water forcing of a regional ocean model simulation. To this aim, different model scenarios, based on the CMEMS IBI model configuration, have been designed and run. The modelled salinity fields, resulting from the different sensitivity IBI test runs, are validated by means of comparison with different in-situ observational data sources. Despite that test runs are performed over the whole IBI model domain, this study mostly focuses on ROFI areas of the Atlantic Iberia and on the French shelf in the Gulf of Biscay.

#### **2. Materials and Methods**

#### *2.1. Description of Study Domain: In-Situ Observed Salinity Products*

The region of study, hereafter referred as IBBIS, is a subset of the IBI operational configuration (Figure 2) and spans from the Strait of Gibraltar (35◦ N) to Bretagne (48.6◦ N). The warm and saline Iberian Poleward Current (IPC, [28]) flows along the steep slope of the narrow shelf of the Western Iberian Peninsula (50 km wide at its maximum, west of 8◦ W), then around the Galician shelf, and finally, enters the Bay of Biscay, separating the continental shelf waters from the open-ocean dynamics. This semi-enclosed Bay of Biscay features a narrow shelf on its southern part that widens towards the North (~150 km at 47◦ N), separated from a deep abyssal plain (more than 4000 m depth) by a steep and rugged slope, with many canyons. The main drivers on the shelf of this IBBIS domain are wind, tides and freshwater coastal runoff. The ocean dynamic varies locally, depending on the orientation of the coast, width of the shelf, freshwater discharge rates and other factors.

In this study, the IBBIS domain has been divided in three subregions, each one featuring its own characteristics: the CADIZ area, that is the Northern shelf of the Gulf of Cádiz, the Western Iberian Shelf (WISHE area), that includes the western shelf of the Iberian Peninsula, and the Gulf of Biscay (BISCA area).

The CADIZ area, at the southern tip of the Iberian Peninsula, is the place where the exchanges between the Atlantic and Mediterranean take place, across the Strait of Gibraltar. The saline MOW (Mediterranean Outflow Waters) come as a subsurface outflow from the Strait of Gibraltar, and cascade along the slope of the narrow shelf of the bay [29,30]. The circulation on the shelf is wind-driven, with a clear seasonal variability [31]: upwelling occurs in summer under the influence of northerlies, while downwelling takes place in

winter under the influence of southerlies [32]. The variability on the shelf is explained by the front between the cold upwelled shelf-waters and the warmer southern waters (e.g., Reference [33]). Two main rivers, Guadalquivir and Guadiana, feed the shelf on its wider portion: in winter, the increased river runoff cools down the shelf waters [34] and contributes to the generation of a westward inshore flow [35], although this is not the main driver.

**Figure 2.** Study Domain. Location of the in-situ measurements used in this study: fixed moorings (purple), surface transects with thermo-salinometer (blue), coastal profiles with CTD (Conductivity, Temperature, and Depth) instrumentation (orange) and profiles with XBT (eXpendable BathyThermograph) instrumentation (red). Black rectangles show subregions in which the study is focused (IBBIS, CADIZ, WISHE, BISCA). The names of the main rivers included in the IBI model set-up and the geographical features mentioned in the paper are denoted in grey. The shaded field shows the IBI averaged surface salinity over the year 2018, from the reference simulation; the 200 m isobath, delimiting the shelf, is depicted in dashed line.

The WISHE area features a narrow shelf (30 to 90 km at its widest), and in comparison to the other proposed study subregions, this is the most sensitive to open-ocean general circulation patterns occurring in the North-eastern Atlantic. The geostrophic flow of the IPC intensifies in winter and isolates the narrow shelf's waters from the open-ocean southward circulation. Along the Galician coast, the currents are mainly wind-driven [36]. The WISHE area marks the limit of the Eastern North Atlantic Upwelling region, with a frequent development of equatorward filaments at the coast in spring and autumn [37]. The interactions and front dynamics between the IPC and the shelf waters are complex [38]. The Western Iberian Buoyant Plume extends from the Mondego river northwards and is characterized by a salinity below 35.8. It is fed by the Mondego, Douro and Minho rivers and the numerous Rias of Galicia. Under non-upwelling conditions, the plume stays confined to the inner shelf and contributes to the stratification and cooling of the shelf waters. On the other hand, under upwelling conditions, the plume extends offshore, and masks the thermal signature of the IPC.

Finally, in the BISCA area, the IPC develops instabilities and eddies along its path which propagate offshore and take part of the general anticyclonic circulation of the openocean section [28,39]. On the French shelf, the main drivers are the coastal density-driven jets which interact with the tidal signal [40,41], and the surface wind-driven currents that advect shelf waters southward, e.g., Solabarrieta et al. [42]. In winter, the activity is weak over the shelf, but intense mesoscale activity takes place in the open ocean. The Southwestern winds induce northeastward drift over the shelf. From spring, the wind weakens and the dynamic activity on the shelf reaches the same levels of energy as offshore [43]. The major source of freshwater on the BISCA shelf is the river outflow, with a seasonal cycle featuring a maximum of discharge in winter and spring. The main rivers are the Loire and Vilaine system (47◦ N), the Garonne (45◦ N) and the Adour, in the South-West corner of the bay. Freshwater filaments extend from the coast and participate in the cross-shelf exchange [44]. Patches of low-salinity surface waters can be pushed by wind-induced coastal upwellings and reach the shelf-break [45]. Along the northern Iberian coast, the numerous rivers have a small averaged runoff, but can have peaks of torrential discharge, principally in spring and autumn. The accumulated outflow from these coastal rivers in winter and spring induces a low-salinity front on the narrow shelf [46], that can reach the shelf-break and extend down to 50 m depth [47].

All the in-situ salinity observations available for the year 2018 were gathered and explored to assess the capacity of the IBI model configurations to reproduce salinity fronts and river plume extension on the shelf, depending on each freshwater forcing. Table 1 shows an overview of the available observations used in each of the study areas.


**Table 1.** Overview of the in-situ observational salinity data sources in the study domain and subregions.

> Hourly timeseries of salinity from mooring stations for the year 2018 were obtained from the CMEMS in-situ observational product (Copernicus Marine In-Situ Tac Data Management Team (2019) [48]). The buoys are mostly moored on the shelf or at the shelfbreak (see their locations, represented by purple dots, in Figure 2) and measurements are collected at surface and subsurface (mostly between 1 and 3 m deep). After processing, 13 stations are available at the whole study area: 7 in the BISCA zone (numbered M1 to M7,

following the coast southwards from Bretagne to Galicia), 5 in the WISHE one (numbered M8 to M12, from North to South) and 1 in the CADIZ subregion (M13). Data from all of them have been used for the model validation, however some stations provide a limited exploitable input; in some cases, there is no more than a few months of data available (e.g., M4), whereas in other cases, the buoy is moored too close to the coast, not being fully useful for a meaningful validation of a regional model such as IBI (e.g., M3 or M10). However, as salinity observational data are sparse, it has been decided not to discard them, but use them with caution.

Weekly salinity profiles from 5 conductivity-temperature-depth (CTD) stations over the year 2018 located at Rias Baixas within the WISHE area were provided by the *Instituto Tecnolóxico para o Control do medio Mariño de Galicia* (INTECMAR). These CTD measures have been running since 1992 and the current INTECMAR oceanographic network is formed by 43 stations distributed along the Galician Coast [49]. These weekly campaigns are very coastal, with most of the stations in Rias and estuaries. The depth values are not uniform over the year, depending on the devices used to measure the data [50]. In this study, only samples from the 5 most offshore INTECMAR CTD stations (identified by orange stars in Figure 2) over the year 2018 are used, and these measurements are compared with daily modelled salinity profiles extracted, at the closet model grid point, from the different IBI model runs.

Surface Transects measured every 10 min with thermo-salinometers (TSG) installed on board of the research vessel Noruega were performed during the year 2018. Data from this specific campaign were provided by the Instituto Português do Mar e da Atmosfera (IPMA). Blue transects depicted in Figure 2 show the on-shelf measurements taken from 42◦ N southwards. This IPMA campaign covers part of the WISHE and CADIZ zones for a short time (from 28 April until 30 May 2018). The salinity is calculated as the salinity of the water inside the TSG and it represents salinity of the oceanic surface layer. Therefore, these measurements are used to validate the surface salinity extracted from the different IBI model simulations performed.

Daily salinity expendable bathythermograph (XBT) profiles are obtained from the Coriolis Ocean Dataset for Reanalysis for the Ireland-Biscay-Iberia region [51] real-time website (www.coriolis.eu.org (accessed on 8 April 2021)). These data are located on the shelf (red dots in Figure 2) in the BISCA area and are a part of the RECOPESCA ( *Réseau de mesure de l'activité de pêche spatialisé et de données environnementales à usage scientifique*) Project campaigns (IFREMER; *L'Institut Français de Recherche pour l'Exploitation de la Mer*), which collect in-situ data in the Bay of Biscay onboard fishing vessels [52]. The salinity XBT profiles are compared to the daily modelled salinity profiles extracted from IBI simulations at the closest grid point.

#### *2.2. Description of the New River Freshwater LAMBDA Database*

The CMEMS Service Evolution LAMBDA project [27] aimed to improve land boundary conditions by coupling watershed models to ocean regional models. The project study area simulated the territories draining to the European Atlantic front and the North Sea, to provide boundary conditions to the CMEMS IBI and NWS (North West Shelves) Monitoring and Forecasting Centers, using the MOHID (Water Modelling System) Land model [53] with regular grids with a 0.05◦ horizontal resolution. The MOHID Land model, part of the MOHID Modeling System (http://www.mohid.com (accessed on 8 April 2021); [54]), is a physically based, spatially distributed, continuous, variable time step model for the water and property cycles in inland waters that allows simulation of single and multi-catchment simulations.

Each drainage basin domain was simulated for a 11-year period (1 January 2008–1 January 2019) using ERA5 reanalysis meteorology from the Copernicus Climate Change Service (https://climate.copernicus.eu/ (accessed on 8 April 2021)). These simulations resulted in the first LAMBDA data product, used in this study, consisting of daily flow and temperature values near the river mouth for a total of 54 main European rivers. In a second phase, 70 extra discharges for minor rivers in Western Iberia and 364 extra discharges for Ireland and UK domain were extracted to obtain a more comprehensive freshwater budget in these areas. However, these extra rivers included in a second release of the LAMBDA database are not fully validated ye<sup>t</sup> and the model resolution (too coarse) used may be a limitation to accurately represent flows from such minor rivers. Therefore, in this work, only used freshwater data for rivers (the ones used currently as forcing of the IBI model application) included in the first release of the LAMBDA database was used. More information about the LAMBDA project products and results can be seen in Campuzano et al. [55].

#### *2.3. The CMEMS IBI Salinity Model Product*

The CMEMS-IBI Near-Real-Time forecast and analysis operational product [56] is generated daily by means of a NEMO model application [57], run over a regional grid (see geographical coverage in Figure 1) with a horizontal resolution of 1/36◦ (around 3 km) and 50 unevenly distributed z vertical levels. A best estimate of the ocean state (IBI analysis generated through a model run with data assimilation on a weekly basis) and a five-dayahead forecast (daily updated free model run) are routinely produced for a number of hydrodynamic variables, among others: temperature, salinity, mixed layer depth, zonal and meridional velocity currents and sea surface height. While hourly-averaged estimations are provided for the sea surface, daily-averaged values are computed for the rest of the three-dimensional water column in the entire IBI regional domain. Hourly 3D physical forecast data, covering exclusively the continental shelf and coastal areas, are additionally delivered to foster the dynamical downscaling and subsequent implementation of coastal downstream services.

The CMEMS IBI system uses a SAM2-based (second generation of Mercator Assimilation System) data assimilation scheme [58] to enhance its predictive skills by constraining the model in a multivariate way with a wealth of observations: Altimeter data (i.e., along-track sea-level anomalies), in-situ temperature and salinity vertical profiles (from the CMEMS CORA—Coriolis Ocean Dataset for Reanalysis- 4.1 database) and satellitederived sea surface temperature (OSTIA—Operational Sea Surface Temperature and Sea Ice Analysis- product) are regularly assimilated to estimate periodically initial conditions for numerical forecast simulation. The resulting product from this latest IBI analysis run is what is delivered as historical best estimates.

At the time of the experiment, the IBI ocean model operational simulation is forced by 3-h daily updated high-resolution (1/8◦) meteorological conditions provided by the European Center for Medium-Range Weather Forecast (ECMWF). Lateral open boundary conditions for temperature, salinity and currents, imposed from daily outputs from the parent system (the CMEMS GLOBAL ocean forecast [59]), are complemented with 11 tidal harmonics built from the FES2014 [60] tidal model solution. Land freshwater forcing is included in the IBI operational model set-up. Discharged flows for 33 of the main rivers presented in the IBI area are considered. In-depth insights into this IBI coastal and riverine freshwater forcing are provided in the next section devoted to describing the IBI model sensitivity tests performed in this study. For further technical specifications of the CMEMS-IBI model solution, the reader is referred to Sotillo et al. [20] and Aznar et al. [61].

The multi-parameter skill assessment of this CMEMS IBI ocean forecast system is operationally conducted through the NARVAL (North Atlantic Regional VALidation) toolbox, which routinely monitors its performance and prognostic capabilities by computing a variety of quality indicators [62]. Both real-time ('online mode') and regular-scheduled ('delayed-mode') comparisons are performed using a wealth of independent observational sources as a benchmark, encompassing in-situ observations (collected by moorings, tidegauges, drifters, gliders and Argo floats networks) and remote-sensed estimations (from satellites and high-frequency radars). For salinity field, the CMEMS-IBI model performance is evaluated using a multiple spatio-temporal scale adopting a multi-platform approach (see list of observational data sources used in the IBI validation in Table 2).


**Table 2.** Summary of observational (Obs) platforms used in NARVAL (North Atlantic Regional VALidation) toolbox to validate IBI model salinity (\* the climatology used is the World Ocean Atlas (WOA) 2013).

> With respect to in-situ observations, the NARVAL software tool makes routine comparisons of CMEMS IBI salinity model products, extracting sea surface salinity timeseries on the grid points closest to the mooring buoy locations, and salinity profiles where Argo or glider observations are available.

> Additionally, the IBI SSS field is validated using daily remotely sensed satellite estimations derived from the European Space Agency's Soil Moisture and Ocean Salinity (SMOS) mission. In recent years, the SMOS products have steadily evolved to mitigate systematic biases detected in coastal areas [63], and today, they better portray plume patterns in major river basins [64,65]. However, in the IBI area, there is still room for improvement in these SMOS-derived products, especially in the operational ones and at very coastal areas where systematic biases are identified. The operational SMOS products' limitations (not able to fully capture variability associated to minor coastal freshwater contributions) have conditioned their use in NARVAL to operationally assess IBI salinity product quality, focusing the validation only on deep open waters.

> It is pointed out that recent updates of SMOS products, taking benefit of new LAMBDA Project developments, has resulted in a new 9-year (2011–2019) reprocessing of global SMOS SSS daily maps. This new dataset introduces significant improvements with respect to the operational and previous SMOS versions, tending to improve salinity gradients, preserving the small-scale ones close to the coast [66], and minimizing latitudinal and seasonal biases [67]. These improvements support the possibility of using this new SMOS satellite data as a salinity reference in model assessments, and the IBI-MFC has introduced it in NARVAL as an extra salinity data source for IBI delayed mode validations.

> Highlights from the operational validation of IBI salinity products with in-situ observations and comparisons with the new SMOS satellite product are shown in Section 3.1.

#### *2.4. Description of the IBI Model Sensitivity Tests Performed Using LAMBDA Product as Forcing*

In the CMEMS IBI operational NEMO model set-up, inputs from 33 rivers are prescribed as a daily flux of freshwater (0.1 psu. The main rivers of the European Atlantic Shelf are considered (e.g., Garonne or Rhine, which have mean outflows larger than 1000 m3s−1), but also smaller rivers which can present a strong temporal variability (e.g., Guadalquivir, whose flow is usually below 200 m3s−1, but has peaks up to 6000 m3s−1). To compensate the input from missing rivers, whose individual debits can be negligible but in total have a significant impact on the coastal salinity, a monthly climatological runoff (hereafter FWF\_CRF) is applied all along the coast. It is derived from the same climatology used in the CMEMS GLOBAL configuration. The overall freshwater input, surface runoff and rivers, in the operational CMEMS IBI configuration during 2018, had a daily mean flow of 16 × 10<sup>3</sup> m3s−1, where river discharge accounted for 67%. The river discharge flows used in the reference CMEMS IBI configuration consist of a composite of observations and climatology (hereafter referred to as FWF\_REF) provided by the GRDC (Global Runoff Data Centre; http://www.bafg.de/GRDC (accessed on 8 April 2021)), the French "Banque Hydro" dataset (http://www.hydro.eaufrance.fr/ (accessed on 8 April 2021)), simulations from the Swedish Meteorological Service (SMHI) E-HYPE ((HYdrological Predictions for the Environment)) hydrological model (http://e-hypeweb.smhi.se (accessed on 8 April 2021)) and observations gathered by the Previmer project. For further technical details on

how freshwater forcing is included in the IBI operational forecast service, the reader is referred to Amo et al. [56].

Ten rivers discharge in the IBBIS area (Figure 1). Four of them are along the French coast (BISCA area), four in the Western Iberian shelf (WISHE area) and two in the Gulf of Cádiz (CADIZ area). Since April 2018, most of the river discharges were forced by a monthly climatology (FWF\_CLM) or by an annual mean flow (in the IBBIS region, only the Mondego river), calculated from the reference forcing. The new forcing set that will be assessed in this study comes from the LAMBDA project watershed product (FWF\_LAM), described above. Table 3 shows the mean daily and cumulated annual flow for each of the forcing datasets in the four areas. Mean daily flow of the LAMBDA forcing for year 2018 in the IBBIS area reaches 8910 m3s−1, that is 70% more than the river input used for the IBI reference run (5223 m3s−1) and 128% more than the climatological forcing (3893 m3s−1).

**Table 3.** Mean daily river discharge (m3s−1) and cumulated discharge (10<sup>3</sup> m3s−1) over 2018 for each forcing set, in the full study domain (IBBIS) and the three subregions (BISCA, WISHE, CADIZ).


The rivers prescribed as forcing in the subregions are: 1 Vilaine, Loire, Garonne, Adour. 2 Minho, Douro, Mondego, Tagus. 3 Guadiana, Guadalquivir.

> Figure 3 shows the daily river discharge of the reference forcing (FWF\_REF, blue line), LAMBDA river database (FWF\_LAM, red line), climatology (FWF\_CLM, cyan dashed line) and the extra coastal runoff (FWF\_CRF, grey line), in the IBBIS domain and the three subdomains (BISCA, WISHE, CADIZ). For the IBBIS region, in the first months of 2018, the variability and intensity are equivalent in both the reference and the new LAMBDA datasets, with peaks of freshwater in January and March, though peaks are not always simultaneous (Figure 3a). From April onwards, the FWF\_REF is mainly based on the climatology and shows a smooth decrease in total river debit from Spring up to August, when it reaches a minimum average daily flow of ≈2000 m3s−1. LAMBDA features a minimum in summer too, associated with a reduced variability, but shows two important peaks in June and November, that do not exist in the IBI reference forcing. The difference between the two datasets is particularly important in BISCA (Figure 3b), where the mean river discharge is 1.8 times more important in the new forcing (Table 3).

> The impact of using new river discharge forcing in CMEMS-IBI, and more specifically on the simulated salinity fields and frontal features linked to river plumes, will be evaluated through specific IBI-like model scenarios. To assess the performance of the numerical model and identify the main sources of uncertainty linked to the river runoff forcing in the model simulations, four scenarios forced by the previously described freshwater discharge parametrization were performed over the year 2018.

> Table 4 shows an overview of the four IBI model land boundary scenarios: IBI\_REF is the reference simulation forced by the reference river discharge forcing, and its parameterization is identical to the one used in the operational IBI configuration, IBI\_LAM differs from IBI\_REF by using the LAMBDA river database forcing, IBI\_NOR is like IBI\_LAM, but without adding the extra coastal runoff (i.e., the river discharge forcing is the only input of freshwater), and finally, IBI\_CLM is a run forced by climatology. The latter scenario was simulated only over the first 6 months of 2018, because the reference run from April is using as forcing the same climatology.

**Figure 3.** Average daily freshwater discharges imposed as forcing in different IBI sensitivity runs in (**a**) IBBIS, (**b**) BISCA, (**c**) WISHE and (**d**) CADIZ areas, in 2018. The freshwater forcing FWF\_REF, FWF\_LAM, FWF\_CLM and FWF\_CRF data sources are respectively depicted in blue, red, dashed blue and grey.



All the IBI scenarios were performed with the same model parametrization and forcing datasets, excluding the freshwater forcing, as the CMEMS-IBI operational configuration. A spin-up run of 110 days was performed to generate the same initial condition for all model scenarios. The model started from rest (the 6 September 2017) with initial three-dimensional temperature and salinity fields taken from the CMEMS-IBI operational outputs. After this spin-up period, the model solution was evaluated as stable, showing on the 1 January 2018 velocity fields similar in magnitude and pattern to the CMEMS-IBI operational solution.

All the outputs from these IBI model scenario runs are delivered at hourly frequency for the sea surface, and at a daily frequency for the rest of the water column. Thus, to validate model salinity outputs, surface hourly model data have been used for comparison with mooring observed timeseries to take advantage of the higher temporal frequency of these available in-situ salinity observations. On the other hand, to compare with other

observational data sources, such as the ones based on CTD or XBT profiles, daily model salinity outputs at the closest grid point and depth are collocated.
