Skill assessment of global, regional, and coastal circulation forecast models: evaluating the benefits of dynamical downscaling in IBI (Iberia–Biscay–Ireland) surface waters

In this work, a multi-parameter inter-comparison of diverse ocean forecast models was conducted at the sea surface ranging from global to local scales in a two-phase stepwise strategy. Firstly, a comparison of CMEMS GLOBAL and the nested CMEMS IBI regional system was performed against satellite-derived and in situ observations. Results highlighted the overall benefits of both the GLOBAL direct data assimilation in open water and the increased horizontal resolution of IBI in coastal areas. Besides, IBI (Iberia–Biscay–Ireland) proved to capture shelf dynamics by better representing the horizontal extent and strength of a river freshwater plume, according to the results derived from the validation against in situ observations from a buoy moored in NW Spain. Secondly, a multi-model inter-comparison exercise for 2017 was performed in the Strait of Gibraltar among GLOBAL, IBI, and SAMPA (Sánchez-Garrido et al., 2013) high-resolution coastal forecast systems (partially nested to IBI) in order to elucidate the accuracy of each system to characterize the Atlantic Jet (AJ) inflow dynamics. A quantitative validation against hourly currents from high-frequency radar (HFR) highlighted both the steady improvement in AJ representation in terms of speed and direction when zooming from global to coastal scales through a multi-nesting model approach and also the relevance of a variety of factors at local scale such as a refined horizontal resolution, a tailored bathymetry, and a higher spatio-temporal resolution of the atmospheric forcing. The ability of each model to reproduce a 2 d quasi-permanent full reversal of the AJ surface inflow was examined in terms of wind-induced circulation patterns. SAMPA appeared to better reproduce the reversal events detected with HFR estimations, demonstrating the added value of imposing accurate meteorologically driven barotropic velocities in the open boundaries (imported from the NIVMAR (Álvarez-Fanjul et al., 2001) storm surge model) to take into account the remote effect of the atmospheric forcing over the entire Mediterranean basin, which was only partially included in IBI and GLOBAL systems. Finally, SAMPA coastal model outputs were also qualitatively analysed in the western Alboran Sea to put in a broader perspective the context of the onset, development, and end of such flow reversal episodes.

1 Introduction

Over the last three decades, significant progress has been made in the discipline of operational oceanography thanks to the substantial increase in high-performance computational resources that have fostered seamless evolution in ocean modelling techniques and numerical efficiency (Cotelo et al., 2018) and given rise to an inventory of operational ocean forecasting systems (OOFSs) running in overlapping regions in order to reliably portray and predict the ocean state and its variability at diverse spatio-temporal scales.

Global circulation models have been steadily evolving in terms of complexity, horizontal resolution refinement, and process parameterization (Holt et al., 2017). Notwithstanding, such development involves compromises of scale and is subject to practical limits on the feasible spatial resolution (Greenberg et al., 2007). Although large-scale physical processes are properly resolved by the current state-of-the-art global model resolution (e.g. nominal 1∕12^∘), coastal and shelf phenomena are still poorly replicated or even misrepresented as the grid mesh is too coarse. This is especially true for complex-geometry regions such as sea straits, archipelagos, or semi-enclosed seas where the coastline, seamounts, and bottom topography are not well resolved. In this context, tides, vertical coordinates, mixing schemes, river inflows, and atmospheric forcings have been traditionally identified as five areas of further research in global ocean modelling (Holt et al., 2017).

Since the continental shelf is affected not only by natural agents (land–sea breezes, riverine discharges, bottom topography, coastline shape, etc.) but also by human-induced factors, an increased understanding of coastal circulation is essential for decision- and policy-making in the socioeconomically vital and often environmentally stressed coastal regions. Therefore, small-scale ocean features must be explicitly computed and accurately reproduced by means of regional models with finer horizontal grid spacing but for a particular delimited area. The success of this approach requires the seamless progress in several aspects, as previously identified by Wilkin et al. (2017) and Kourafalou et al. (2015): (i) a deep comprehension of the primary mechanisms driving coastal circulation, (ii) downscaling methods to adequately represent air–sea and land–sea interactions, and (iii) robust methods to embed high-resolution models in coarser-scale systems. Therefore, this approach implies the transfer of large-scale information from the global model to the interior of the nested regional domain by means of diverse methodologies. One of them is the so-called “spectral nudging” technique, adopted to ensure that the prevailing global conditions are not degraded in the open ocean, while allowing submesoscale processes to be resolved exclusively by the nested model in the continental shelf and coastal areas (Herbert et al., 2014). An alternative approach to computational time-demanding multiple nesting procedures consists of using unique unstructured grid models as they have proved to properly describe ocean processes at different spatial scales (Ferrarin et al., 2013, 2019; Federico et al., 2017; Cucco et al., 2012). They benefit from imposing higher resolution in shallow water areas to better resolve irregular coastlines, intricate bathymetries, and hence small-scale dynamics while applying a coarser grid resolution in open water to reproduce large-scale phenomena. Besides, unstructured grid models are particularly effective to achieve the seamless transition between adjoining basins interconnected by narrow straits (Ferrarin et al., 2018; Stanev et al., 2017; Zhang et al., 2016).

Additionally, the regional modelling strategy can include some fine-tuning of physical parameters, individually tailored to each chosen area, instead of the universally valid parameterizations associated with global OOFSs. The benefits of regional modelling over the driving global OOFS are generally assumed but to date only few studies have explored and quantified the potential added value of such approaches (Katavouta and Thomson, 2016; Rockel, 2015; Greenberg et al., 2007). The “parent–son” model inter-comparison is mandatory during both implementation and operational stages since it aids to (i) verify the most adequate nesting strategy, (ii) check the consistency of the nested model solution, and (iii) identify any potential problem that might be inherited from the coarser system.

In the framework of the Copernicus Marine Environment Monitoring Service (CMEMS), a global ocean model together with a wealth of nested regional OOFSs are currently running in different areas of the European seas and providing paramount oceanographic forecast products (Le Traon et al., 2018). Since the validation of OOFSs against independent measurements constitutes a core activity in oceanographic operational centres, the skill of Iberia–Biscay–Ireland (IBI) regional OOFS is routinely assessed by means of the NARVAL (North Atlantic Regional VALidation) system (Sotillo et al., 2015), a web-based toolbox that provides a series of skill metrics automatically computed and delivered in the QUality Information Document – QUID – (Sotillo et al., 2018). In this context, the first goal of this paper is to conduct a multi-parameter model inter-comparison between IBI regional OOFS and the coarser parent system, the CMEMS GLOBAL (Lellouche et al., 2018), with the aim of assessing their performance at the upper ocean layer. Their predictive skills to properly represent the sea surface temperature (SST) over IBI coverage domain and diverse subregions were evaluated by means of comparisons against remote-sensed and in situ observations. On the other hand, their prognostic capabilities to accurately reproduce the coastal surface circulation were assessed through the analysis of a single impulsive-type river outflow episode that took place in March 2018 in the Galician coast (NW Spain), a region of freshwater influence – ROFI – (Simpson, 1997).

Despite the recent advances in the development of CMEMS global and regional core products, many downstream services for user uptake require information on even smaller spatial scales, such as ocean forecasting for small island chains (Caldeira et al., 2016), intricate bights (Stanev et al., 2016), or port approach areas where sharp topo-bathymetric gradients pose special difficulties for accurate local predictions (Hlevca et al., 2018; Federico et al., 2017; Sánchez-Arcilla et al., 2016; Sammartino et al., 2014; Grifoll et al., 2012). A variety of operational products for harbours have been recently developed, although most of these coastal applications are wave and water-level forecasting systems (Lin et al., 2008; Pérez et al., 2013). By contrast, less attention has been devoted to hydrodynamic harbour conditions since its reduced dimensions and intricate layout confer harbour restrictions which are not present in the open sea. Besides, derivative products based on current forecasts, such as float trajectories, residence time maps, flushing patterns, and risk assessment of water quality degradation, can constitute additional assets for efficient harbour management (Álvarez-Fanjul et al., 2018; Sammartino et al., 2018). In order to overcome the existing gap between the scales effectively solved by the regional OOFSs and the coastal scales required to meet strong societal needs in support of blue and green growth, a number of downstream services are currently adopting different downscaling approaches. Dynamical downscaling takes regional boundary conditions to drive a high-resolution limited-area model in which coastal processes are calculated on a finer grid by resolving well-known hydrodynamic equations. However, uncertainties in the downscaling process must be evaluated since coastal model performance can be directly impacted by the propagation of any potential issue in the large-scale dynamics, inherited from the coarser system (Hernández et al., 2018).

As a representative example of downstream service developed by Puertos del Estado (PdE) in a hot spot area like the Strait of Gibraltar (GIBST), the operational PdE-SAMPA high-resolution coastal system (Sánchez-Garrido et al., 2013) is partially embedded in IBI and nowadays employed by the Port Authority of Algeciras Bay as predictive tool to support maritime policy and assist with high-risk decision-making related to marine safety, port operation optimization, and mitigation of both natural disasters and anthropogenic hazards. Previous research has unequivocally proved the ability of PdE-SAMPA to accurately capture basic circulation features of the GIBST area and Algeciras Bay (Sanchez-Garrido et al., 2014; Sammartino et al., 2014; Soto-Navarro et al., 2016). A preliminary model skill assessment was conducted within the framework of the MEDESS-4MS project (Sotillo et al., 2016a). However, the added value of this coastal OOFS with respect to the regional IBI system was only quantified from a Lagrangian perspective by using a wealth of drifters. The second goal of this contribution is thus to build upon previous model inter-comparison exercises, placing special emphasis on the characterization of the Atlantic Jet (AJ) inflow into the Mediterranean Sea in terms of speed and direction. This geostrophically adjusted jet fluctuates on a wide range of temporal scales and drives the main circulation in the Alboran Sea, feeding and surrounding the Western Alboran Gyre – WAG – (Macias et al., 2016). An inter-comparison exercise was conducted for 2017 among a global configuration (CMEMS GLOBAL), a regional application (CMEMS IBI) and a higher resolution coastal system (PdE-SAMPA) in order to characterize the AJ dynamics and their ability to adequately capture an extreme event: the quasi-permanent (up to ∼48 h long) full reversal of the AJ surface flow under intense and prolonged easterlies. To this end, a high-frequency radar (HFR) has been used as benchmark since it regularly provides quality-controlled hourly maps of the surface currents of the strait (Lorente et al., 2014). A detailed characterization of this unusual phenomenon is relevant from diverse aspects, encompassing search and rescue operations (to adequately expand westwards the search area), the management of accidental marine pollution episodes (to establish alternative contingency plans), or safe ship routing (to maximize fuel efficiency).

In summary, this paper serves one primary purpose: performing a multi-parameter model skill assessment in IBI surface waters, ranging from global to local scales in a two-phase stepwise strategy. This includes (i) a comparison between GLOBAL and IBI regional systems in the entire overlapping coverage domain, focusing special attention on regionalization and (ii) an event-oriented multi-model inter-comparison for 2017 with a focus on the complete inversion of the surface flow in the GIBST. This process-based validation approach, albeit commonly used in meteorology and weather forecasting, is rather novel in operational oceanography and mostly devoted to extreme sea level and wave height episodes.

This paper is organized as follows: Sect. 2 provides further details about the study areas, Sect. 3 describes the diverse models configuration, Sect. 4 outlines the observational data sources and methodology used in this study, and Sects. 5 and 6 present a detailed discussion of the results. Finally, main conclusions are summarized in Sect. 7.

2 Study areas

2.1 IBI area (and subregions)

From a pure physical oceanographic point of view, the IBI geographical domain is a very complex region (Fig. 1a), marked by a generally steep slope separating the deep ocean from the shelf. The western, and deeper, side of the IBI domain is affected by main large-scale currents, mainly the closure of the North Atlantic Drift. It splits into two major branches: the major one continuing northwards along the northwestern European shelves (North Atlantic Current, NAC, and North Atlantic Drift Current, NADC) and the other, the Azores Current (AC), which follows southeastwards and has continuity in the Canary Current (CaC). On the other hand, along the slope, a poleward slope current flows in the subsurface; it is observed as far north as at latitudes near Ireland. Instabilities in this slope current favour the occurrence of slope water oceanic eddies along the northern Iberian coast (Pingree and Le Cann, 1992). On the continental shelves, intense tidal motions provide the dominant source of energy (Álvarez-Fanjul et al., 1997): noticeable tidal mixing fronts arise on the most energetic tidal areas of the IBI region (i.e. English Channel, Celtic and Irish seas). Shelf and coastal areas of the region are also affected by strong storm surges (Pérez et al., 2012). Along the western Iberian and African coasts, strong summer upwelling of bottom cold and enriched waters takes place under predominant northerly wind conditions that trigger the Ekman-driven offshore deflection of the surface flux.

Figure 1 (a) Iberia–Biscay–Ireland service (IBISR) domain, which comprises nine subregions denoted by red squares. Filled red dots represent buoy locations. (b) Study area: coverage domain of SAMPA coastal ocean model, where the surface Atlantic Jet (AJ) flows through the Strait of Gibraltar into the Alboran Sea, feeding the Western Alboran Gyre (WAG); isobath depths are labelled every 200 m. Red dots indicate a topographic feature: Camarinal Sill (CS). (c) HFR hourly data availability for 2017: solid black squares represent radar sites, blue and red dots indicate Tarifa tide gauge and B7 buoy location, respectively. The black line denotes the selected transect and the white square represents its midpoint.

IBI is also a rather broad and heterogeneous area. In order to gain insight into the model skill assessment (as later explained in Sect. 5), the IBI service (IBISR) regional domain has been split into nine different subregions (Fig. 1a): the Irish Sea (IRISH), the English Channel (ECHAN), the Gulf of Biscay (GOBIS), the North Iberian Shelf (NIBSH), the West Iberian Shelf (WIBSH), the Western Mediterranean Sea (WSMED), the Gulf of Cádiz (CADIZ), the Strait of Gibraltar (GIBST), and the Canary Islands (ICANA).

2.2 Strait of Gibraltar

The Strait of Gibraltar (GIBST), the only connection between the semi-enclosed Mediterranean basin and the open Atlantic Ocean (Fig. 1b), is characterized by a two-layer baroclinic exchange which is hydraulically controlled at Camarinal Sill (Sánchez-Garrido et al., 2011). Whilst saltier Mediterranean water flows out at depth, an eastward surface jet of relatively fresh Atlantic water (AJ) flows into the Alboran Sea by surrounding the quasi-permanent Western Anticyclonic Gyre (WAG) and the more elusive Eastern Anticyclonic Gyre (EAG) in a wavelike path. As the WAG owes its existence to the input of new Atlantic waters provided by the AJ, both structures are widely considered to be coupled and usually referred to as the AJ–WAG system. A significant variety of analytical, field, and modelling studies have previously attempted to disentangle the AJ–WAG system and properly explain the underlying physical processes (Sánchez-Garrido et al., 2013; Macías et al., 2007a; Viúdez, 1997).

The position, intensity, and direction of the AJ fluctuate on a broad range of temporal scales, driving the upper-layer circulation of the Alboran Sea with subsequent physical and biological implications (Solé et al., 2016; Sánchez-Garrido et al., 2015; Ruiz et al., 2013). For instance, the presence of a strong AJ close to the northern shore of the Alboran Sea reinforces the coastal upwelling and therefore increases both the near-shore chlorophyll concentration and the spawning of fish in this region (Ruiz et al., 2013; Macías et al., 2008). By contrast, meteorologically induced inflow interruptions can trigger weakening and even the decoupling of the AJ–WAG system (Sánchez-Garrido et al., 2013), the subsequent eastward migration of the WAG, and the genesis of a new gyre that coexists with the other two, giving rise to a three-anticyclonic-gyre situation (Viúdez et al., 1998).

Within this context, the AJ pattern has been described to oscillate between two main circulation modes at seasonal scale (Vargas-Yáñez et al., 2002): (i) a stronger AJ flows northeastwards during the first half of the year and (ii) a weaker AJ flows more southwardly towards the end of the year. Sea level pressure (SLP) variations over the western Mediterranean basin and local zonal wind (U) fluctuations in the Alboran Sea are usually considered the main factors controlling and modulating the AJ variability (Macías et al., 2007b; García-Lafuente et al., 2002a, b). In particular, the second parameter has been largely invoked as the primary driving agent to explain both the intensification of the surface inflow during prevalent westerlies and also extreme AJ collapse events recorded when intense easterlies are predominant (Macías et al., 2016). The zonal wind intensity has been reported to follow an annual cycle with more westerly (easterly) winds during winter (summer) months (Dorman et al., 1995). The seasonal variability and occasional interruptions of the Atlantic inflow due to meteorological forcing have been earlier investigated with in situ data from fixed moorings (García-Lafuente et al., 2002a). More recently, a considerable number of satellite-tracked drifters were released on both sides of GIBST within the framework of the MEDESS-4MS project, hence providing a complete Lagrangian view of the Atlantic water inflow into the Alboran Sea (Sotillo et al., 2016-b).

3 Model descriptions

Whereas basic features of the three OOFSs employed in this work are gathered in Table 1, further details are provided in the following devoted subsections.

Table 1 Basic features of the ocean forecast systems employed in the present study. ^*The operational version of IBI used with spectral nudging. Assimilation scheme SAM2 was later introduced in version 4 (April 2018).

3.1 CMEMS GLOBAL system

The operational mercator global ocean analysis and forecast system provides 10 d of 3-D global ocean forecasts updated daily. This product includes daily mean files of temperature, salinity, currents, sea level, mixed layer depth, and ice parameters from the surface to seafloor over the global ocean. It also includes hourly mean surface fields for sea level height, temperature, and currents. The global ocean output files are displayed with a 1∕12^∘ horizontal resolution with regular latitude and longitude equirectangular projection. Some 50 vertical levels span from 0 to 5500 m.

The product is updated as follows: everyday, the daily configuration is run with updated atmospheric forcings, without assimilation, for days D−1 to D+9. The daily runs are initialized with the previous day's run, except on Thursdays, when they start from the weekly analysis run. Every week, on Wednesdays, the weekly configuration is run with assimilation for days D−14 to D−1. This run is separated into two parts: a best analysis for days D−14 to D−8 and an analysis for days D−7 to D−1. Therefore, every day, the time series is updated with new forecasts for days D−1 to D+9, erasing the previously available data for D−1 to D+8. In addition, on Thursdays, the analysis is also provided, replacing previously available files for days D−14 to D−1. For further details, the reader is referred to the GLOBAL Product User Manual – PUM – (Law Chune et al., 2019).

The system is based on the Nucleus for European Modelling of the Ocean (NEMO) v3.1 ocean model (Madec, 2008). The physical configuration is based on the tripolar ORCA (Madec and Imbard, 1996) grid type with a horizontal resolution of 9 km at the Equator, 7 km at Cape Hatteras (mid-latitudes), and 2 km towards the Ross and Weddell seas. The 50-level vertical discretization retained for this system has 1 m resolution at the surface decreasing to 450 m at the bottom, and 22 levels within the upper 100 m. The bathymetry used in the system is a combination of interpolated ETOPO1 (Earth TOPOgraphy 1 arcmin) and GEBCO8 (GEneral Bathymetric Chart of the Oceans, version 8) databases. The system was initialized on 11 October 2006 based on the temperature and salinity profiles from the EN4 monthly gridded climatology. The atmospheric fields forcing the ocean model are taken from the ECMWF (European Centre for Medium-Range Weather Forecasts) Integrated Forecast System. A 3 h sampling is used to reproduce the diurnal cycle. The system does not include tides nor pressure forcing. The monthly runoff climatology is built with data on coastal runoffs and 100 major rivers from the Dai et al. (2009) database (Lellouche et al., 2018). Altimeter data, in situ temperature and salinity vertical profiles, and satellite sea surface temperature data are jointly assimilated to estimate the initial conditions for numerical ocean forecasting. Moreover, satellite sea ice concentration data are now assimilated in the system in a mono-variate or mono-data mode. More information can be found in Lellouche et al. (2018).

3.2 CMEMS IBI regional system

The IBI OOFS provides a real-time short-term 5 d hydrodynamic 3-D forecast (and 1 d of hindcast as best estimate) of a range of physical parameters (currents, temperature, salinity, and sea level) since 2011 (Sotillo et al., 2015). IBI is based on an eddy-resolving NEMO model application (v3.6) that includes high-frequency processes required to characterize regional-scale marine processes. The model application is run at 1∕36^∘ horizontal resolution and final products are routinely delivered in a service domain extending between 26–56^∘ N and 19^∘ W–5^∘ E. The NEMO model (Madec, 2008) solves the three-dimensional finite-difference primitive equations in spherical coordinates discretized on an Arakawa-C grid and 50 geopotential vertical levels (z coordinate), assuming hydrostatic equilibrium and Boussinesq approximation. Partial bottom cell representation of the bathymetry (a composite of ETOPO2 and GEBCO8) allows an accurate representation of the steep slopes characteristic of the area. The model grid is a subset of the Global 1∕12^∘ ORCA tripolar grid used by the parent system (the CMEMS GLOBAL system) that provides initial and lateral boundary conditions (west – the Atlantic Ocean, east – the Western Mediterranean Sea), but refined at 1∕36^∘ horizontal resolution.

The IBI run is forced every 3 h with up-to-date high-frequency (1∕8^∘ horizontal grid resolution) meteorological forecasts (10 m wind, surface pressure, 2 m temperature, relative humidity, precipitation, shortwave and longwave radiative fluxes) provided by ECMWF. CORE empirical bulk formulae (Large and Yeager, 2004) are used to compute latent sensible heat fluxes, evaporation, and surface stress. Lateral open boundary data are interpolated from the daily outputs of the GLOBAL system. These are complemented by 11 tidal harmonics built from FES2004 (Lyard et al., 2006) and TPXO7.1 (Egbert and Erofeeva, 2002) tidal model solutions. Atmospheric pressure component, missing in the large-scale parent system sea level outputs, is included assuming pure isostatic response at open boundaries (inverse barometer approximation). Fresh water river discharge inputs are implemented as lateral open boundary condition for 33 rivers. Flow rate data imposed are based on a combination of daily observations from PREVIMER, simulated data from E-HYPE hydrological model, and monthly climatological data from Global Runoff Data Centre (GDRC) and French “Banque Hydro” datasets. Further details can be found in Sotillo et al. (2015).

Originally, the operational IBI system was based on a periodic reinitialization from the GLOBAL parent solution. Afterwards, IBI has steadily evolved: by April 2016, an upgrade of the downscaling methodology was implemented, substituting the periodic reinitialization by a spectral nudging technique in order to avoid temporal discontinuity inherent to the periodic reinitialization and to minimize dependency from the GLOBAL parent solution on the shelf.

The spectral nudging aims at forcing the regional child model solution (IBI) to be close to the parent system (GLOBAL) in those areas where the latter is supposed to be accurate thanks to data assimilation, mainly in deep waters outside the continental shelf (Herbert et al., 2014). After each forecast cycle of IBI, the increment between the parent model analysis X_p and the child model forecast X_c is calculated for a chosen state variable X (typically currents, temperature, and salinity). Such an increment consists of a space and time low-pass filter of the differences X_p–X_c in order to keep the characteristic scales so the parent GLOBAL system can properly resolve (from large scale to mesoscale). After each forecast cycle of the child IBI system, a new cycle (called analysis) is relaunched where X_c is nudged to X_p at each time step by the weekly mean of the daily increments previously computed in the parent system grid during the forecast cycle. The 1-week time smoothing window was imposed as this is the typical timescale associated with mesoscale structures.

Furthermore, the nudging is spatially limited in those areas where the parent system can not improve the regional model (e.g. where there is no data assimilation of altimetry or where the physics is missing, for instance on the shelf) or where the spatial filtering processes are potentially detrimental (close to the bottom or the open boundaries). This spatial weight function is a 3-D mask (for further details, see Herbert et al., 2014) that differentiates the zones where IBI system is nudged and thus reconciled with the parent system (typically in the open sea) from those where IBI remains free (continental shelf, coastal areas, and regions close to the open boundaries) to fully compute higher-frequency processes such as thermal tidal fronts, river discharges, etc.

Finally, a SAM2-based data assimilation scheme (Lellouche et al., 2013; Brasseur et al., 2005) was recently introduced (April 2018) in order to enhance IBI predictive skills but will not be further described here as only outputs from 2017 have been used in the present work.

3.3 PdE SAMPA coastal system

The PdE-SAMPA operational forecast service started in April 2012 (Sammartino et al., 2014; Sánchez-Garrido et al., 2014). It routinely provides a daily short-term forecast (72 h horizon) of currents and other oceanographic variables in the Strait of Gibraltar and its surroundings (Gulf of Cádiz and Alboran Sea). The PdE-SAMPA model application was developed by the University of Málaga in collaboration with PdE in order to provide a tailored forecasting service to one of their main stakeholders, the harbour of Algeciras Bay (Fig. 1b). It is based on the Massachusetts Institute of Technology global circulation model – MITgcm – (Marshall et al., 1997). The domain, which extends from the Gulf of Cádiz to the Alboran Sea (Fig. 1b), is discretized with an orthonormal curvilinear grid of variable horizontal resolution, sparser close to the boundaries (∼8–10 km), and higher in the strait (∼300–500 m). In the vertical dimension, SAMPA has 46 unevenly spaced z levels with maximum resolution of 5 m near the surface, exponentially decaying towards the seafloor. The shallower level is at 2.5 m depth. The bathymetry is derived from a combination of the GEBCO bathymetry dataset and fine-resolution bathymetric charts of the Strait of Gibraltar and the continental shelf of the Gulf of Cádiz and northern coast of the Alboran Sea. The bottom topography is represented as partial vertical cells. In the two lateral open boundaries (west – the Atlantic Ocean, east – the Western Mediterranean Sea) the model is partially forced by daily mean temperature, salinity, and velocity fields from CMEMS IBI regional model (Sotillo et al., 2015). Since such frequency is not suitable to resolve barotropic flows through the strait either (García-Lafuente et al., 2002b), tidal and meteorologically driven barotropic velocities are prescribed across the open boundaries: the former extracted from the Mog2D model described by Carrere and Lyard (2003) and the latter from the storm surge operational system developed by Álvarez-Fanjul et al. (2001), which accounts for the remote effect of the atmospheric forcing in the barotropic flow through GIBST. This nesting strategy ensures that the SAMPA model captures a realistic variability in inflow and outflow currents through the strait. At the sea surface, the model is forced by hourly values of wind stress, air humidity and temperature, fresh water, and heat surface fluxes provided by the Spanish Meteorological Agency through the operational Forecast System based on the HIRLAM model (Cats and Wolters, 1996). Further details on the SAMPA model configuration are provided in Sanchez-Garrido et al. (2013).

4 Validation of OOFSs

4.1 Framework

The validation of OOFSs against independent measurements constitutes a core activity in oceanographic operational centres since it aids (i) to infer the relative strengths and weaknesses in the modelling of several key physical processes, (ii) to compare different versions of the same OOFS and evaluate potential improvements and degradations before a new version is transitioned into operational status, (iii) to compare coarse-resolution “father” and nested high-resolution “son” systems to quantify the added value of downscaling.

With regards to the third aspect, IBI forecast products are regularly compared not only against other CMEMS regional model solutions (e.g. NWS and MED) in the overlapping areas but also against its parent system (GLOBAL) by means of NARVAL (North Atlantic Regional VALidation) login-protected web-based application (Sotillo et al., 2015). This tool has been implemented to routinely monitor IBI performance and to objectively compare model reliability and prognostic capabilities. Both real-time validation (“online mode”) and regular-scheduled “delayed-mode” validation (for longer time periods) are performed using a wealth of observational sources as benchmark: in situ observations from buoys and tide gauges, SST satellite-derived products, temperature and salinity profiles from ARGO floats, and HFR. Product quality indicators and skill metrics are automatically computed in order to infer IBI accuracy and the spatio-temporal uncertainty levels. The evaluation metrics regularly generated by NARVAL are delivered online in the QUID, which is periodically updated and freely available on the CMEMS website (http://marine.copernicus.eu/, last access: 1 July 2019).

Complementarily, opportunistic inter-comparisons are conducted in the frame of diverse EU-funded projects such as MEDESS-4MS (Sotillo et al., 2016a): 35 satellite-tracked drifters were released on both sides of the Strait of Gibraltar and the quality-controlled in situ data of sea surface temperature and currents were collected to build the MEDESS-GIB database (Sotillo et al., 2016b), hence providing a complete Lagrangian view of the surface inflow of Atlantic waters through the GIBST and the Alboran Sea. Such valuable oceanographic information was subsequently used to compare IBI and SAMPA forecast products to identify strengths (realistic simulation of the Atlantic Jet and the Algerian Current) and shortcomings (position and intensity of the Alboran gyres, especially the western one) in both model performances. This exercise reflected the effectiveness of the dynamical downscaling performed through the SAMPA system with respect to the regional solution (in which SAMPA is partially nested, as it is also embedded in Mog2D and NIVMAR models), providing an objective measure of the potential added value introduced by SAMPA.

Eventually, ancillary validation approaches have been recently adopted focused on the evaluation of ocean model performances in specific situations and on their ability to accurately reproduce singular oceanographic processes (Hernández et al., 2018). Since the NARVAL tool is devoted to compare model solutions on a monthly, seasonal, or annual basis, part of the picture is missing due to traditional time averaging. Hence the quality indicators computed, albeit valid, mask somehow model capabilities to replicate ocean phenomena of particular interest at shorter timescales. This event-oriented multi-model inter-comparison methodology allows us to better infer the ability of each system to capture small-scale coastal processes. In this context, the recurrent question “Which model is the best one?” should be reformulated by firstly admitting that one system can outperform the rest of the OOFSs for a particular event but by contrast can be also beaten when attempting to reproduce and characterize some other distinct ocean phenomenon.

Those oceanographic events subject of further insight might encompass: (i) coastal upwelling, dowelling, and relaxation episodes; (ii) submesoscale eddies (Mourre et al., 2018); (iii) extreme events; and (iv) complete flow reversals. Particularly, in the present work, attention has been devoted to the full and permanent reversal of the surface AJ in the GIBST during, at least, 48 h. This unusual episode has been detected by means of HFR current estimations and further examined with OOFS outcomes. The agreement between both in situ and remote-sensing instruments as well as the ocean forecasting system has been evaluated by means of computation of a set of statistical metrics traditionally employed in this framework: histograms, bias, root-mean-square differences (RMSDs), scalar and complex correlation coefficients, current roses, histograms, quantile–quantile (QQ) plots, and the best linear fit of scatter plots. In the following subsection, all the in situ and remote-sensed observations employed in the present work are described.

4.2 Observational data sources

4.2.1 In situ observations

The study domain includes an array of buoys operated by Puertos del Estado and the Irish Marine Institute (Fig. 1a), providing quality-controlled hourly-averaged observations of SST, sea surface salinity (SSS) and currents. To ensure the continuity of the data record, occasional gaps detected in time series (not larger than 6 h) were linearly interpolated. Basic features of each in situ instrument are described in Table 2.

Table 2 Description of the network of directional buoys used in this work. Year label stands for year of deployment. Subregions are defined in Fig. 1a.

4.2.2 Satellite-derived observations

The European Ocean Sea Surface Temperature L3 Observations is a CMEMS operational product which provides a daily fusion of SST measurements from multiple satellite sensors over a 0.02^∘ resolution grid. The L3 multi-sensor (super-collated) product is built from bias-corrected L3 mono-sensor (collated) products. If the native collated resolution is N and N<0.02^∘ the change (degradation) of resolution is done by averaging the best quality data. If N>0.02^∘ the collated data are associated to the nearest neighbour without interpolation nor artificial increase in the resolution. A synthesis of the bias-corrected L3 mono-sensor (collated) files remapped at resolution R is done through a selection of data based on the following hierarchy: AVHRR_METOP_B, SEVIRI, VIIRS_NPP, AVHRRL-19, AVHRRL-18, MODIS_A, MODIS_T, AMSR2. This hierarchy can be changed in time depending on the health of each sensor. Further details can be found in the Product User Manual (PUM), freely available on the CMEMS website (http://cmems-resources.cls.fr/documents/PUM/CMEMS-SST-PUM-010-009.pdf, last access: 1 July 2019)

4.2.3 HFR-derived observations

The HFR system employed in the present study consists of a three-site shore-based CODAR SeaSonde network, installed in GIBST (Fig. 1b and c). Hereafter the sites will be referred to by their four-letter site codes: CEUT, CARN, and TARI (Fig. 1c). Each site is operating at a central frequency of 26.8 MHz, providing hourly radial current measurements which are representative of the upper 0.5 m of the water column. The maximum horizontal range and angular resolution are 40 km and 5^∘, respectively. Radial current measurements from the three stations are geometrically combined with an averaging radius set to 3 km, in order to estimate hourly total current vectors on a Cartesian regular grid of 1 km×1 km horizontal resolution.

A source of error to be considered in the computation of the total vectors is the so-called geometrical dilution of precision (GDOP). The GDOP is defined as a dimensionless coefficient of uncertainty that characterizes how radar system geometry may impact on the measurement accuracy and position determination errors, owing to the angle at which radial vectors intersect. Maps of east and north GDOP for this HFR system (not shown) follow a pattern where their values increase with the distance from the radar sites and along the baselines (lines connecting two HFR sites), as the combining radial vectors are increasingly parallel and the orthogonal component tends to zero. Further details can be obtained from Lorente et al. (2018).

The accuracy of HFR measurements, which are affected by intrinsic uncertainties (radio frequency interferences, environmental noise, etc.) have been previously assessed by comparing against in situ observations provided by a pointwise current meter (Lorente et al., 2014), yielding correlations above 0.7 and RMSD below 13 cm s⁻¹. Such results revealed that this HFR network has been operating within tolerance ranges, properly monitoring the surface circulation in near-real time of this geostrategic region.

Recent work relying on this HFR system has successfully investigated the water exchange between Algeciras Bay and the Strait of Gibraltar (Chioua et al., 2017), the impact of the atmospheric pressure fluctuations on the mesoscale water dynamics of the Strait of Gibraltar and the Alboran Sea (Dastis et al., 2018), the dominant modes of spatio-temporal variability in the surface circulation (Soto-Navarro et al., 2016), or the characterization of the Atlantic surface inflow into the Mediterranean Sea (Lorente et al., 2018).

In the present work, quality-controlled hourly HFR current measurements collected during the entire 2017 were used as benchmark to elucidate the skill of a number of OOFSs. The data availability was significantly high: almost 100 % in the selected transect (solid black longitudinal line shown in Fig. 1c), decreasing in the easternmost sectors. The transect used to examine the AJ surface inflow was readily chosen as the associated total GDOP, reported in Lorente et al. (2018); was reduced (below 1.3); and the spatial and temporal data availability were optimal during 2017. From an oceanographic perspective, the election of such a transect was also convenient to better characterize both the intensity and direction of the AJ, since its midpoint covers the area where the highest peak of current speed is usually detected and also where the inflow orientation is not influenced yet by the water exchange between Algeciras Bay and the Strait of Gibraltar.

5 Comparison between CMEMS model solutions in IBI waters

5.1 Temperature

The CMEMS L3 satellite-derived daily data were used to validate the SST fields predicted by both GLOBAL and IBI systems. The map of annual availability of this remote-sensing product for 2017 (Fig. 2a) reveals that the highest percentages of available observations (above 80 %) were found in the south of Canary Islands, the Gulf of Cádiz, and in the Western Mediterranean Sea. Equally, there was also a significant data provision in the west coast of the Iberian Peninsula and Morocco, although it decreased to 70 % in nearshore areas over the shelf such as the Iberian and African upwelling or the Strait of Gibraltar. By contrast, the lowest data availability (below 40 %) was detected in the northernmost latitudes, including the Irish and North seas.

Figure 2 (a) Availability of satellite-derived L3 SST daily data for 2017; (b, c) annual mean absolute difference (MAD) of SST for 2017 in open water where the spectral nudging is applied: GLOBAL versus satellite and IBI versus satellite; (d, e) Analogous to (b, c) but in coastal waters where no spectral nudging is applied. Zones of worse model performance (i.e. higher MAD) are delimited with blue rectangles. Spatially averaged MAD values are provided in the lower-right corner.

Maps of annual mean absolute differences (MADs) were firstly computed in open water where the spectral nudging technique was applied (Fig. 2b and c). Apparently, both models behaved similarly during 2017, although GLOBAL performance was slightly better due to the direct data assimilation scheme implemented, as reflected by a lower spatially averaged MAD of 0.13^∘. In the open sea, IBI benefited indirectly from the data assimilation conducted in its parent system thanks to the spectral nudging technique (MAD = 0.15^∘), thereby allowing the regional and global model states to be reconciled with each other. IBI even outperformed GLOBAL locally in specific zones (delimited with blue rectangles in Fig. 2b) such as the western Canary Islands, the periphery of Madeira Island, or a portion of the African coastal upwelling system. Furthermore, a narrow belt of SST anomalies could be observed along the continental shelf break in the case of GLOBAL (Fig. 2b). Since IBI presents a higher grid resolution, it could partially resolve the internal waves breaking, which leads to turbulence and energy for increased vertical mixing with cooler waters beneath the pycnocline, ultimately contributing to the reduced SST differences observed in IBI estimations over the continental shelf break (Fig. 2c). By contrast, GLOBAL appeared to outperform IBI in the Gulf of Cádiz, the NW Iberian open water and the region between Canary and Madeira islands (Fig. 2c).

Maps of annual MAD in coastal areas (Fig. 2d and e), where no spectral nudging was imposed and the IBI system ran freely, exhibited some similarities such as the pronounced SST differences encountered in the Iberian and African coastal upwelling systems (ICUS and ACUS, respectively). According to the skill metrics, spatially averaged over the entire coastal waters domain, IBI performance was, on average, more accurate than GLOBAL (0.17^∘ versus 0.20^∘) thanks to several factors (i.e. the more refined mesh and the tidal solution included in IBI). GLOBAL outputs exhibited relevant SST differences over the continental shelf, especially in energetic tidally dominated areas such as the English Channel, the Irish Sea, and southern part of the North Sea (Fig. 2d). In these regions, monthly maps of SST bias (not shown) revealed the alternation between winter cold anomalies and summer warm anomalies, in accordance with earlier results of Graham et al. (2018). Tidally driven vertical mixing could account for a portion of the discrepancies encountered between the coarser detided GLOBAL and IBI model solutions, where the former seems to predict an over-stratification in shelf seas. IBI also appeared to slightly outperform GLOBAL in some delimited areas of the Portuguese coast, the Gulf of Cádiz, and the ACUS (Fig. 2d). Besides, the higher SST anomalies observed for GLOBAL in the Strait of Gibraltar and the western Alboran Sea were likely linked to the inadequate representation of the speed and direction of the Atlantic Jet inflow, something that will be addressed in the following sections. On the contrary, it is also true that GLOBAL was to a little extent more precise in some parts of Western Mediterranean Sea, the ICUS, and ACUS (Fig. 2e).

The MAD metrics for each subregion within the IBI regional domain (defined in Fig. 1a) are shown in Fig. 3. As above mentioned, IBI generally outperformed GLOBAL system in those subregions where no spectral nudging was applied (English Channel, Irish and North seas, North Iberian Shelf), with the exception of the Western Iberian Shelf (WIBSH). In addition, better metrics were obtained for IBI in the Gulf of Biscay (GOBIS), where the spectral nudging was only applied in the westernmost off-shelf area. On the contrary, GLOBAL seemed to better replicate the SST field in the open sea thanks to the direct data assimilation. IBI metrics, albeit consistent, were slightly worse than those obtained for GLOBAL in this region.

Figure 3 Annual (2017) mean absolute difference (MAD) of SST, spatially averaged for each subregion (defined in Fig. 1a), between IBI (GLOBAL) and satellite L3 observations, represented by blue (black) columns. IBI-SN-ON (IBI-SN-OFF) represents the area where the spectral nudging technique was (was not) applied.

Although both models appeared to better fit to observations in the open sea, a portion of the SST differences detected in this region might be attributed to the fact that the satellite data assimilated into the GLOBAL system (CMEMS OSTIA gap-free product) are different from the independent (not-assimilated) satellite-derived observations (CMEMS L3 gappy product) used as benchmark to comprehensively validate the models. Furthermore, satellite products are affected by intrinsic uncertainties. The quality information documents (QUIDs) focused on the accuracy assessment of OSTIA (McLaren et al., 2016) and L3 (Saux-Pickart et al., 2019) satellite products have reported a RMSD of 0.4^∘ (0.2^∘) between OSTIA (L3) and drifting buoys observations for 2012 (2018). More specifically, in the case of OSTIA estimations, the RMSD is 0.40^∘ globally with regional values ranging from 0.28^∘ in the South Pacific to higher values in the North Atlantic (0.47^∘) or in the Mediterranean Sea (0.89^∘). By contrast, the L3 SST product shows reasonably good difference statistics against drifting buoy measurements. This bias in satellite estimations should be taken into account when interpreting the results of the model validations.

On the other hand, the highest SST discrepancies for both models were generally located in very coastal areas (African and Iberian upwelling systems, Strait of Gibraltar, etc.) and over the continental shelf (the Irish Sea and the English Channel) where satellite remote sensing could be complicated by weather patterns and dissolved organic compounds of terrestrial origin that may attenuate signals and yield unreliable results (Thakur et al., 2018). As a consequence, many processed remote sensing products apply a land mask that excludes mixed pixels in nearshore areas, use temporal averaging to account for missing observations, or even impose an optimal interpolation process to transform the original satellite tracks into a regular grid (i.e. OSTIA product). In this context, the availability of the SST L3 data (Fig. 2a) was lower on the shelf near coastal areas featuring a complex bathymetry (likely due to poor satellite coverage, application of a land mask, or cloud cover), justifying to some extent the predominance of SST anomalies near the shorelines.

For the sake of completeness, supplementary validation work in the entire 3-D water column with Argo floats is regularly conducted to assess model vertical structure. For the period 2012–2016 and the entire IBI domain, the averaged RMSD for full profiles of temperature and salinity are 0.51^∘ and 0.13 PSU, respectively (for further details, we refer the reader to the IBI QUID). Both models perform fairly well in open water, given the fact that GLOBAL assimilates these types of in situ observations and subsequently transfers the information to the nested IBI system thanks to the aforementioned spectral nudging technique. Nevertheless, validation focused on smaller scales and high-frequency processes is still crucial to analyse in detail the performance of both modelled products in intricate coastal regions.

Hourly in situ observations from eight buoys, moored within specific subregions (Fig. 1a), were used as benchmark to validate both GLOBAL and IBI outputs. The annual time series of SST exhibited a significantly high resemblance, properly reproducing the expected annual-cycle shape (Fig. 4). According to the consistent skill metrics derived from the comparison against three deep-water buoys (B3, B4, and B5, in Table 2), both models had a rather similar performance during 2017 with RMSD and correlation coefficients in the ranges 0.44–0.96 ^∘C and 0.86–0.99, respectively (Fig. 4c–e). While the similar behaviour observed off the shelf is partially attributable to the aforementioned spectral nudging technique, the model–observation comparison in near-shore areas revealed noticeable discrepancies.

Figure 4 Annual (2017) time series of hourly sea surface temperature (SST) at eight different locations within IBISR area. In situ observations from moored buoys (blue dots), GLOBAL model predictions (green line), and IBI model outputs (red line) are depicted. Skill metrics derived from model–observation comparison are gathered in black boxes.

On the one hand, IBI appeared to outperform GLOBAL system in the Irish Sea (Fig. 4b), Gulf of Cádiz (Fig. 4f), and GIBST subregions (Fig. 4g), as reflected by lower (higher) RMSD (correlation) values obtained. Particularly, the results for the Strait of Gibraltar are not in complete accordance with the statistics previously derived from the comparison against L3 satellite-derived data (Fig. 3j), likely due to the fact that remote-sensed SST estimation areas might be affected by higher intrinsic uncertainties (i.e. land contamination and cloud cover). Although both comparisons against remote and in situ observations confirmed the model SST overestimation in GIBST, especially during summertime, the former (latter) indicated that IBI precision was significantly lower (higher). Another relevant aspect is the notable ability of IBI to capture sharp summer SST rise (steeper than 3 ^∘C) during prevalent easterlies (Fig. 4g), as a result of the surface inflow reversal and subsequent intrusion of warmer Mediterranean waters into GIBST (this phenomenon will be subject of further analysis in Sect. 7). However, GLOBAL appeared to overestimate SST in this area during the entire year, as reflected by a RMSD of 1.64 ^∘C.

On the other hand, GLOBAL seemed to behave slightly better at B1 location – IBISR area – (Fig. 4a) and substantially more accurately at B8 buoy location – in the Canary Islands, ICANA, – where a permanent SST overestimation from June to December was evidenced in IBI predictions (Fig. 4h), thereby yielding a RMSD twice higher than that obtained for GLOBAL estimations, in agreement with Figs. 2i–j and 3i. The lower performance of IBI in ICANA subregion was previously reported by Aznar et al. (2016) when inter-comparing IBI forecast and 1∕12^∘ reanalysed solutions. At this point it is worth recalling that GLOBAL includes a data assimilation scheme, whereas IBI takes realistic ocean conditions from weekly global analyses. This fact shows up the possible benefits of the observational data assimilation in these areas, at least in terms of surface variables. Furthermore, a fraction of observed model–buoy discrepancies in SST can be explained in terms of disparate depth scales. Whereas IBI and GLOBAL daily outputs are representative of the temperature in the upper 1 m of the water column, moored buoys provide temperature estimations at a deeper nominal depth (between 1 and 3.5 m, depending on the brand). Future validation exercises should include the interpolation of model outputs to both the exact buoy location and also to the specific depth level in order to more accurately assess the model skilfulness. It is worthwhile mentioning that a variety of previous work has focused on the comparison of in situ and remote-sensed SST data, reporting significant differences across different geographical regions (Thakur et al., 2018; Stobart et al., 2015; Smit et al., 2013). Such discrepancies could also be observed in this work as lower biases were generally derived from the pointwise comparisons, thereby highlighting the intrinsic uncertainties associated with satellite data in nearshore areas.

Complementarily, a quarterly analysis was performed to infer any potential degradation in model performances during a specific season of the year (Fig. 5). Overall, both GLOBAL and IBI predictions seemed to be more reliable in winter (except at B1 location, Fig. 5a) in terms of lower RMSD. They also emerged to be less realistic during summer, as denoted by abrupt decreases in quarterly correlation indexes (from 0.9 down to 0.5) at B2, B4, and B6 locations and the relevant rise of RMSD (up to 2.5 ^∘C) at B7 location (GIBST subregion). This SST overestimation could be partially explained in terms of imprecise latent sensible heat fluxes and excess of evaporation, although additional efforts should be devoted to shed light on it. Once again, IBI performance appeared to be more accurate in coastal zones featuring a more complex bathymetry (at B2, B4, B6, and B7 locations), whereas GLOBAL fitted better to in situ observations in off-shelf regions such as at B1 and B8 locations. In the rest of the cases, both model solutions were rather similar. It is noteworthy that each pointwise buoy is not representative of the entire subregion in which it is deployed, thus explaining to some extent the discrepancies arisen between Sect. 6.1 and 6.2.

Figure 5 Annual evolution (2017) of seasonal skill metrics derived from the comparison of GLOBAL and IBI models against in situ SST hourly observations provided by eight buoys. RMSD and correlation coefficient represented by coloured bars and lines, respectively.

5.2 Salinity

As pointed out in the introduction, the enhancement of riverine forcing is still a priority in ocean modelling as the estuarine circulation is mainly driven by horizontal density gradients which are ultimately modulated by freshwater inputs. In this context, previous work has investigated the potential benefits of replacing old climatologies by data from hydrological model predictions (O'Dea et al., 2017). Here we provide a specific example to illustrate the discrepancies between GLOBAL and IBI performances in the Galician coast (NW Spain), as a consequence of the different horizontal resolution and distinct runoff forcing implemented in the operational chain. While both model performances are rather similar in open water (according to the results derived from the validation against 3-D Argo-float profiles and exposed in the QUID), higher discrepancies are expected to arise in coastal and shelf areas as they are governed by small-scale processes such as land–sea breezes, runoff (and the resulting stratification and buoyancy-driven circulation), and transport materials (nutrients, sediments, pollutants, etc.).

As shown in Fig. 6a, hourly in situ SSS data collected by B4 buoy during March 2018 experienced an abrupt decrease from a standard value around 36 PSU down to almost 33 PSU in just a few hours during 20 March, likely due to a noticeable filament of freshwater discharged by Miño River. IBI outputs at the closest grid point appeared to properly capture both the sharp drop in SSS values and the persistent low salinity values for the next 4 d period. By the end of the month, the modelled salinity field seemed to steadily recover to usual levels in the range of 35.5–35.8 PSU, whereas in situ observations revealed a steeper rise to 34.8 PSU by 23 March. Nevertheless, the skill metrics confirmed the rather accurate IBI performance, with a correlation coefficient of 0.92 and a RMSD of 0.33 PSU. By contrast, although GLOBAL outputs could replicate the mean SSS, it did not reproduce satisfactorily the freshwater episode and barely showed any temporal variability, as reflected by a negligible correlation coefficient (0.09) and a higher RMSD (0.84 PSU).

Figure 6 (a–d) Monthly inter-comparison (March 2018) between GLOBAL (green line), IBI (red line / red triangles) and B4 buoy (blue dots): sea surface salinity (SSS), temperature (SST), current direction (SCD), and current speed (SCS); (e, f) daily maps of SSS and SST derived from GLOBAL outputs for the 21 March. Red filled dot represents B4 buoy location; (g, h) daily maps of SSS and SST derived from IBI outputs for the 21 March; (i) monthly surface current roses. Monthly skill metrics derived from model–observation comparisons are provided.

Consequently, the impact of colder freshwater river inputs on the SST was also evaluated (Fig. 6b). Once again, while the sudden cooling of 1.5 ^∘C denoted by in situ observations was fairly well replicated by IBI, the GLOBAL system could only correctly predict the overall decreasing trend along with the SST values immediately before (13.5 ^∘C) and after (13 ^∘C) the analysed event. As a consequence, the monthly correlation coefficient (RMSD) obtained for IBI is higher (lower): 0.79 versus 0.20 (0.25 ^∘C versus 0.35 ^∘C). It is worth mentioning that we could not find any relevant evidence that could justify that the observed SST drop was also partially attributable to a large-scale process (i.e. coastal upwelling event). The wind rose (derived from B4 buoy observations) and map of HFR-derived surface currents for May 2018 (not shown) did not exhibit the well-documented conditions associated with NW Iberian upwelling episodes, such as the predominant northerly winds and the subsequent offshore deflection of the surface flux, respectively.

The buoyancy input introduced by large freshwaters fluxes (particularly during the spring freshet), together with topographic effects, contributed to the development of the well-documented Western Iberian Buoyant Plume (Peliz et al., 2002; Otero et al., 2008), which strongly influenced the shelf circulation, forming an averaged veering to ∼270^∘ (measured clockwise from the north) during 20–21 March, as reflected by in situ observations and IBI outputs (Fig. 6c). However, GLOBAL could only partially reproduce the prevailing surface flow as modelled currents were mainly advected to the south-southwest (180–270^∘). Equally, IBI appeared to correctly replicate the acceleration of the upper-layer stream from 10 to 45 cm s⁻¹ due to impulsive-type freshwater river outflow already observed in situ estimations of sea surface currents (Fig. 6d). Notwithstanding, GLOBAL current intensity remained moderate (below 20 cm s⁻¹) during most of March, including the selected episode, as reflected by the poorer skill metrics obtained. The current speed underestimation observed in this tidal environment is mainly attributable to the fact that GLOBAL system provides a detided solution, so barotropic tidal velocities do not contribute to the final prescribed total velocity.

Daily-averaged maps of modelled SSS and SST were computed for 21 March (Fig. 6e–h) to infer the differences between GLOBAL and IBI. As it can be seen, the former showed a relatively smoothed and spatially homogeneous decrease in the salinity and temperature fields along the entire coastline (Fig. 6e and f), while the latter exhibited more intricate patterns with many filaments together with a significant drop in SSS and SST (Fig. 6g and h) in the periphery of the three main local river mouths (from north to south: Miño, Douro, and Tagus) as a result of freshwater plumes outflowing over saltier Atlantic waters. In the three cases, the SST field could effectively act as a tracer for the salinity stratification.

There is a significant resemblance between the monthly current roses derived from in situ observations and IBI predictions in terms of speed and mean direction (Fig. 6i), showing the predominance of the so-called Iberian Poleward Current, flowing northwards and circuiting the western and northern Iberian margins under prevailing southerly winds (Torres and Barton, 2006). GLOBAL current outputs differed from observations, exhibiting an overall tendency for eastward directions. The skill metrics derived from time series comparison at B4 buoy location confirmed that the regional OOFS outperformed the global one during March 2018, hence postulating the benefits of improved horizontal resolution to better resolve the plume dynamics and its extension off shelf. In addition, the increased horizontal resolution of IBI allows us to better resolve individual frontal fluctuations and horizontal salinity gradients by preserving the signal of river plume closer to the coast and with a more complex structure. The impact of model resolution in both the horizontal extent of the plume and the strength and position of the freshwater front has been subject of previous studies (Bricheno et al., 2014). Since both models present 50 depth levels and similar vertical discretization, the horizontal resolution and the riverine forcing are assumed to play a primary role when attempting to explain the differences encountered in model performance for this specific test case.

6 Circulation in the Strait of Gibraltar: multi-model comparison from global to coastal scales

Proved the relevance of the intensity and orientation of the AJ in determining the surface circulation of the Alboran Sea, the ability of each OOFS to portray the upper layer circulation in the GIBST area has been evaluated. The annually averaged surface pattern provided by the HFR network revealed northeastward speeds around 100 cm s⁻¹ in the narrowest section of the strait (Fig. 7a). SAMPA coastal model seemed to capture well the time-averaged intensity and orientation of the Atlantic inflow (Fig. 7b), whereas IBI regional model clearly overestimated the mean surface circulation speed (Fig. 7c). Finally, the coarser OOFS (GLOBAL) barely captured the most basic features on the incoming flow and its subsequent propagation towards the northeast (Fig. 7d).

Figure 7 Annual mean circulation pattern in GIBST for 2017, derived from hourly estimations provided by (a) HFR, (b) SAMPA coastal model, (c) IBI regional model, and (d) GLOBAL model. For the sake of clarity, only one vector for every two was plotted in the HFR map.

As this qualitative model intercomparison on a yearly basis was insufficient to infer the skilfulness of each system, a quantitative validation at the midpoint of the selected transect (white square in Fig. 1c) was assessed. The scatter plot of HFR-derived hourly current speed versus direction (taking as reference the north and positive angles clockwise) revealed interesting details (Fig. 8a): firstly, the AJ flowed predominantly eastwards, forming an angle of 78^∘ with respect to the north. The current velocity, on average, was 100 cm s⁻¹ and reached peaks of 250 cm s⁻¹. Speeds below 50 cm s⁻¹ were registered along the entire range of directions. Westwards currents, albeit in the minority, were also observed and tended to predominantly form an angle of 270^∘.

Figure 8 (a–d) Annual (2017) scatter plot of hourly AJ current speed versus direction (angle measured clockwise from the north); estimations provided by (a) HFR, (b) SAMPA, (c) IBI, and (d) GLOBAL. Mean and standard deviation values of both AJ speed and direction are gathered in black boxes; (e–g) annual scatter plot of differences (observation minus model) in AJ speed and direction between (e) HFR and SAMPA, (f) HFR and IBI, and (g) HFR and GLOBAL.

The scatter plot of SAMPA estimations presented a significant resemblance in terms of prevailing current velocity and direction (Fig. 8b). Although the time-averaged speed and angle were slightly smaller (90 cm s⁻¹) and greater (88^∘), respectively, the main features of the AJ were qualitatively reproduced: maximum velocities (up to 250 cm s⁻¹) were associated with an eastward flow and an AJ orientation in the range of 50–80^∘. Besides, surface flow reversals to the west were properly captured.

By contrast, noticeable differences emerged in the scatter plot of regional IBI estimations (Fig. 8c): surface current velocities below 30 cm s⁻¹ were barely replicated and the AJ inversion was only observed very occasionally. Despite the fact that IBI appeared to properly portray the mean characteristics of the eastwards flow, the model tended to favour flow directions comprised between 60 and 180^∘ and to overestimate the current velocity, with averaged and maximum speeds around 117 and 280 cm s⁻¹, respectively.

In the case of the scatter plot derived from GLOBAL estimations, even more substantial discrepancies were detected as the variability in both the AJ direction and speed were clearly limited to the range 65–80^∘ and 50–200 cm s⁻¹, respectively (Fig. 8d). No flow reversals were detected and peak velocities of the eastward flow were underestimated.

The scatter plots of observation–model differences provided relevant information (Fig. 8e–g). In the case of SAMPA, discrepancies were clustered around zero for both parameters, with an asymptotic distribution along the main axes (Fig. 8e). On the contrary, a negative bias to negative differences was observed for both IBI (Fig. 8f) and GLOBAL (Fig. 8g), especially for the latter. In other words, the regional and global OOFSs overestimated both the current speed and the angle of the AJ, reflecting a tendency for more southeasterly directions (clockwise rotated with respect to the north). Overall, a steady improvement in the AJ characterization is evidenced in model performance when zooming from global to coastal configurations, highlighting the benefits of the dynamical downscaling approach along with other relevant factors such as a more detailed bathymetry, a higher spatio-temporal resolution of the atmospheric forcing, or the inclusion of accurate tidal and meteorologically driven barotropic velocities prescribed across the open boundaries. These results highlight the relevance of the remote forcing in the flow reversals and how subtidal barotropic lateral forcing imposed in SAMPA, obtained from NIVMAR storm surge model, ensures that SAMPA model captures a realistic variability in inflow and outflow currents through the strait. Therefore, the remote barotropic effect of the meteorological forcing over the entire Mediterranean basin, which is only partially included in IBI and GLOBAL systems, plays a major role not only in the regulation of the seasonal cycle of the AJ but also in the occasional inflow inversions, in accordance with previous work (Macias et al., 2016). Since the atmospheric pressure forcing is missing in the GLOBAL sea level outputs, this system only considers the barotropic response to wind forcing. In the case of IBI system, both atmospheric forcings are taken into account but the inverted barometer approximation is solely imposed over the Western Mediterranean Sea (not over the entire basin). Therefore, only a portion of the subtidal variability in the flow through the Strait of Gibraltar can be adequately explained (García-Lafuente et al., 2002b).

Additional statistical indicators were computed: two histograms illustrated the number of hourly zonal (U) and meridional (V) velocity data per class interval (Fig. 9a and b). HFR-derived zonal velocity estimations exhibited a Gaussian-like shape clustered around 84 cm s⁻¹ and slightly shifted to lower values in the case of SAMPA coastal model (79 cm s⁻¹). Both datasets show similar positive skew and variability, with the standard deviation around 56–57 cm s⁻¹ for 2017 (Fig. 9a). IBI and GLOBAL presented narrowed histograms, with distributions positively shifted and constrained to zonal velocities above 0 and 40 cm s⁻¹, respectively. In the case of meridional currents, each distribution exhibits a nearly symmetrical Gaussian-like shape but shifted towards different values (Fig. 9b). Whilst SAMPA and IBI exhibited a similar distribution (and moderately similar to that revealed for HFR estimations), GLOBAL histogram emerged again dramatically shortened and restricted only to positive values, revealing a recurrent predominance of the AJ to flow northeastwards.

Figure 9 Annual (2017) histogram of hourly (a) zonal current velocities and (b) meridional current velocities, as provided by HFR, SAMPA, IBI, and GLOBAL. Mean and standard deviation values are gathered in black boxes. Quantile–quantile plots of hourly (c) zonal current velocities and (d) meridional current velocities, as derived from the observation–model comparison. The 10 %–99 % quantiles were established (red filled dots).

Based on the QQ plot for the zonal velocity component (Fig. 9c), it can be concluded that SAMPA estimations were consistent despite the slight overestimation observed for the highest velocities (95th–100th percentiles). The general IBI overestimation along the entire range of percentiles was also clearly evidenced. In accordance with its histogram, GLOBAL system overestimated (underestimated) zonal currents below (above) the 90th percentile. A similar behaviour was also observed for GLOBAL meridional velocities, this time around for the 20th percentile (Fig. 9d). On the contrary, both SAMPA and IBI appeared to generally underestimate the meridional surface current speed, even more for higher percentiles.

Class-2 skill metrics, gathered in Table 3, were also computed in order to provide a quantitative perspective of model performance at the midpoint of the selected transect (white square in Fig. 1c). SAMPA clearly outperformed both systems, as reflected by lower RMSD values for both velocity components together with a complex correlation coefficient (CCC) and phase (CCP) of 0.79 and −8^∘, respectively, which means that SAMPA predictions were highly correlated with HFR current observations although slightly clockwise rotated (i.e. more southeastwards). The agreement between HFR hourly data and IBI and GLOBAL estimations, albeit significant (CCC above 0.6), was lower as the related phase values decreased substantially (especially for GLOBAL: CCP below −20^∘), indicating a more zonal surface flow.

Table 3 Skill metrics derived from the 1-year (2017) validation of sea surface currents estimated by three operational forecasting systems against HFR-derived observations at the midpoint of the selected transect in the Strait of Gibraltar (Fig. 1c).

The three systems predicted more precisely the zonal velocity component than the meridional one, with scalar correlations emerging in the ranges 0.68–0.83 and 0.15–0.56, respectively. Notwithstanding, RMSD values were more moderate for the latter (below 37 cm s⁻¹) than for the former (below 53 cm s⁻¹). This could be attributed to the extremely intense and predominant west–east zonal exchange of Atlantic-Mediterranean waters through GIBST, with the meridional flow playing a residual role.

The statistical results derived from SAMPA–HFR comparison, gathered in Table 3, are in line with those earlier obtained in a 20-month validation performed by Soto-Navarro et al. (2016), which reported correlations of 0.70 and 0.27 for the zonal and meridional velocities, respectively. The observed model–radar discrepancies might be attributed to the fact that the uppermost z level of SAMPA model is 2.5 m, while HFR observations are representative of the first 0.5 m of the water column and thus more sensitive to wind forcing. This might explain some model drawbacks detected in relation to the reduced energy content in surface current speeds, as reflected by the positive bias between HFR estimations and SAMPA outputs (Table 3).

Complementarily, the multi-model inter-comparison exercise in the GIBST region focused on the ability to adequately reproduce an extreme event: the quasi-permanent full reversal of the AJ surface flow during, at least, 48 h when intense easterly episodes were prevalent. Under this premise, only four episodes were detected and categorized during the entire 2017 (Fig. 10). The prevailing synoptic conditions were inferred from ECMWF predictions of sea level pressure (SLP: Fig. 10a–d) and zonal wind at 10 m height (U−10: Fig. 10e–h). A significant latitudinal gradient of SLP was observed in three episodes (February, March, and December), with high pressures over the Gulf of Biscay and isobars closely spaced in GIBST, giving rise to very strong easterlies (above 10 m s⁻¹), channelled through the strait (Fig. 10e, f, and h). In August, the typical summer weather type was observed with Azores High pressures governing the Atlantic area and moderate but persistent easterly winds blowing through the entire Western Mediterranean Sea (Fig. 10c and g).

Figure 10 Two-day averaged synoptic maps of (a–d) sea level pressure (SLP) and (e–h) zonal wind at 10 m height (U-10), provided by ECMWF, corresponding to each of the four Atlantic inflow reversal events analysed during 2017. (i–l) Monthly time series of SLP, spatially averaged over the Western Mediterranean Sea (WSMED) subregion, marked with a big red box in the maps of the first row; (m–p) monthly time series of U-10, spatially averaged over the Strait of Gibraltar (GIBST) subregion, marked with a small red box in the maps of the second row; (q–t) monthly Hovmöller diagrams of HFR-derived zonal current velocity at the selected transect. Red (blue) colour represent eastward (westward) flow; (u–x) monthly time series of hourly sea surface height (SSH) provided by Tarifa tide gauge, represented by a blue dot in Fig. 1. Two-day episodes of permanent flow reversal are marked with black boxes in (i–x).

Both atmospheric variables were spatially averaged over specific subregions (WSMED and GIBST, respectively, indicated by a red square in Fig. 10a–h) and 3-hourly monitored along the selected months (Fig. 10i–p). Very high SLP values and extremely high (negative) U−10 (i.e. intense easterlies) led to a complete inversion of the surface flow from the prevailing eastward direction to a westward outflow into the Atlantic Ocean, as reflected in the Hovmöller diagrams computed for the HFR-derived zonal currents (Fig. 10q–t). In February, a brief 24 h inversion (related to less intense easterlies) preceded the full reversal of the surface flow (Fig. 10q). Likewise, the event detected in March consisted of an abrupt interruption and complete reversal of the eastwards AJ (Fig. 10r). By contrast, in August and December, the classical AJ intense inflow (above 100 cm s⁻¹) into the Mediterranean was only observed in the southern part of the transect, whereas in the northern sector some fluctuations between weaker eastward and westward currents were evidenced, mainly associated with changes in the prevalent wind regime (Fig. 10s and t). Under persistent easterlies, a weaker coastal counter current was detected flowing westwards and bordering the Spanish shoreline (Fig. 10s). Such coastal flow inversion has been previously reported and subject to further analysis by Reyes (2015). Particularly, the flow reversal detected in August was not triggered by high SLP (Fig. 10k) but induced by moderate and persistent easterlies (5 m s⁻¹, Fig. 10o).

Short-lived reversals of the surface inflow have been previously reported to occur almost every tidal cycle in Camarinal Sill (western end of GIBST: Fig. 1b) mainly due to the contribution of the semidiurnal tidal component M₂ (Reyes, 2015; Sannino et al., 2004; García-Lafuente et al., 1990; La Violette and Lacombe, 1988). Since the mean inflow of Atlantic water is modulated by barotropic tidal currents, hourly-averaged sea surface height (SSH) observations provided by Tarifa tide gauge (Fig. 1c) were used to elucidate if the four 2 d inflow reversal events in the eastern end of the strait could have been mostly influenced by spring–neap tidal cycle fluctuations (Fig. 10u–x). Although the fortnightly variability was clearly observable in a monthly time series of SSH, no cause–effect relationship could be visually inferred from the inspection of zonal velocities at the selected transect (Fig. 10q–t). Apparently, evidence of preference for a specific tidal cycle was not observed as the four flow reversal episodes took place under strong easterlies but during different tidal conditions, ranging from neap tides (Fig. 10u) to spring tides (Fig. 10v and x). As shown in Lorente et al. (2018), tides seemed to play a secondary role by partially speeding up or slowing down the westward currents, depending on the phase of the tide. These results are in accordance with previous modelling studies (Sannino et al., 2004) where the contribution of the semidiurnal tidal component to the transport was proved to be relevant over the Camarinal Sill, (incrementing the mean transport by about 30 % for both the inflow and the outflow), whereas it was almost negligible at the eastern end of the strait.

The observed 2 d averaged HFR-derived circulation patterns associated with the four events here studied are depicted in Fig. 11a, e, i, and m. Some common peculiarities were exposed, such as the overall westward outflow through the narrowest section of GIBST or the subtle anticyclonic inflow into the Algeciras Bay. Three study cases revealed a predominant circulation towards the west together with a marked acceleration of the flow in the periphery of Algeciras Bay, reaching speeds above 70 cm s⁻¹ (Fig. 11a, e, and i). The fourth case (December 2017) was substantially less energetic and exhibited a rather counter-clockwise recirculation in the entrance to GIBST (Fig. 11m). On the other hand, two episodes illustrated how the circulation in the easternmost region of the study domain followed a clockwise rotation (Fig. 11e and m).

Figure 11 Two-day averaged maps of the surface circulation in GIBST, corresponding to each of the four Atlantic inflow reversal events detected in 2017 (from top to bottom). Maps derived from hourly estimations were provided by (from left to right): HFR, SAMPA coastal model, IBI regional model, and GLOBAL model. For the sake of clarity, only one vector for every two was plotted in the HFR map.

From a qualitative perspective, SAMPA was able to reproduce fairly well at least two of the four inversion episodes in terms of overall circulation pattern in GIBST and adjacent waters (Fig. 11f, j, and n). In the event of March, SAMPA replicated the intense eastern anticyclonic gyre, with velocities up to 80 cm s⁻¹, along with the inflow into the Algeciras Bay. However, the model could only partially resolve the AJ inversion, exhibiting a counter-clockwise recirculation with the outflow restricted to the northwestern Spanish shoreline (Fig. 11f). In the episode corresponding to 4–5 December (Fig. 11n), the upper-layer dynamic was rather similar to the one previously described for March, albeit less vigorous. The visual resemblance with the HFR map (Fig. 11m) was generally high, according to common features observed: the eastern anticyclonic gyre, the central belt of currents circulating towards the northwest and eventually the cyclonic recirculation structure in the entrance to GIBST. On the contrary, in the event that occurred between 14–15 August (Fig. 11j), a moderate observation–model resemblance was deduced in the northeastern sector of the domain: SAMPA was able to resolve the observed southwestward stream, the inflow into the Algeciras Bay, and the weak intrusion of Mediterranean waters into GIBST bordering the northern shoreline but, by contrast, it was ultimately impelled to join the general AJ inflow governing the strait and propagating towards the east. Finally, although SAMPA predicted the occurrence of AJ reversal by 20—21 February (Fig. 11b), the simulated circulation structure partially differed from that observed with HFR estimations (Fig. 11a). Whereas the former prognosticated a meander-like circulation, a predominant cross-shore stream within the channel and a flow inversion uniquely circumscribed to the entrance of GIBST; the latter provided an overall westward outflow from the Mediterranean Sea into the Atlantic Ocean.

In the case of IBI, the Atlantic inflow was always present. In two episodes, the intense AJ was directed towards the northeast (Fig. 11g and o), converging with the overall clockwise gyre that dominated the easternmost region, which was already observed in the HFR estimations (Fig. 11e and m). By contrast, in the two remaining episodes the surface inflow was predominantly zonal (Fig. 11c) and directed southeastwards (Fig. 11k). Whereas in the former event no common features could be observed between the HFR and IBI, in the latter a moderate observation–model resemblance was deduced in the northeastern sector of the domain, as similarly occurred with SAMPA estimations (Fig. 11j). Leaving aside the counter-clockwise eddy observed in IBI pattern (Fig. 11k), absent from HFR map (Fig. 11i), IBI partially resolved the observed southwestward flow, the circulation into the Algeciras Bay, and the westward penetration of surface waters along the northern shoreline of the strait. Finally, GLOBAL system barely replicated the HFR-derived circulation patterns as the northeastward stream was permanently locked, showing further reduced speed variations from one episode to another (Fig. 11d, h, l, and p).

The reversal of the surface inflow is caused by meteorology-driven flows through the strait associated with the passage of high-pressure areas over the Mediterranean (García-Lafuente et al., 2002a). Because these flows originate in the far field and not in the strait itself, the different grid resolutions of IBI and SAMPA do not appear to be the likely explanation for these events to not show up in the IBI model. Instead, their different skill in capturing such extreme events seems to be associated with their different forcing.

Among the physical implications of the surface inflow reversal, abrupt increases in the SST field were revealed, especially during summertime when warmer surface waters flowed out into the Atlantic from the Mediterranean (Fig. 12). During August 2017, the aforementioned CCC raised on day 11 and lasted until the end of the month, confined at higher latitudes except for the already analysed 2 d event of 14–15 August, coinciding with the full reversal mentioned (Fig. 12a). The monthly inter-comparison of the zonal currents at the midpoint of the selected transect (represented by a black square in Fig. 12a) confirmed the progressive improvement in the skill metrics obtained (Fig. 12b, right box) thanks to both the multi-nesting strategy and the inclusion of accurate tidal forcing. SAMPA and IBI were able to accurately reproduce the wide tidal oscillations, although only the former could properly capture the flow inversions represented by negative zonal velocities that took place between the 14 and 15 August and between 21 and 24 August (Fig. 12b), as SAMPA properly resolves the meteorologically driven (barotropic) currents through the strait, imported from NIVMAR storm surge model. GLOBAL detided outputs only reproduced basic features of the surface flow, showing always smoothed eastward velocities. As a consequence, skill metrics for the coastal OOFS were better than for the regional system, and recursively skill metrics for IBI were in turn better than global ones in terms of higher (lower) correlation (RMSD) values. Analysis for the meridional velocity component (not shown) revealed similar results with the SAMPA outperforming the coarser models. Notwithstanding, the three OOFS proved to be more skilled to forecast zonal rather than meridional currents. The complex correlation coefficient and the related phase were 0.85 and −7.37^∘, respectively (Fig. 12b, black box in the right side), indicating both the relevant SAMPA–HFR agreement and the slight veering of model outputs with respect to HFR estimations: a negative value denoted a clockwise rotation of modelled current vectors (i.e. a more southwardly direction). In the case of IBI, although the phase was similar (−7.92^∘) the complex correlation was lower (0.72). GLOBAL current vectors were, on average, significantly veered clockwise (−25.71^∘), despite the high complex correlation coefficient (0.70).

Figure 12 (a) Monthly Hovmöller diagram of HFR-derived zonal current velocity at the selected transect in the Strait of Gibraltar for August 2017. Red (blue) colour represent eastward (westward) flow. A complete Atlantic inflow reversal episode marked with black box for the 14–15 August; (b) monthly times series of zonal current velocity at the midpoint of the transect (represented by a black square in the Hovmöller diagram) provided by HFR (blue dots), SAMPA (black dots), IBI (red dots), and GLOBAL (green dots); (c) monthly time series of SST provided by B7 buoy (blue dots), SAMPA (black line), IBI (red line), and GLOBAL (green line). Monthly skill metrics derived from observation–model comparison are gathered in black boxes on the right.

From 11 to 17 August, a progressive warming of 7.5 ^∘C at the upper ocean layer of the northern shoreline was observed (Fig. 12c) according to the in situ estimations provided by B7 buoy (whose latitude is located with a solid black dot in Fig. 12a). As easterly winds progressively dominated the study area and persisted enough, the CCC broadened and the complete inflow reversal transported warmer Mediterranean waters to the west through the entire transect, as reflected by the pronounced SST maximum (∼25 ^∘C) detected soon afterwards by 18 August. A secondary peak of SST was monitored by 25 August, before the CCC started weakening. In accordance with previous statements about model behaviour for the zonal currents, once again SAMPA outperformed the coarser systems as reflected by a significantly high correlation of 0.89 and a lower but statistically relevant RMSD of 1.22 ^∘C. IBI presented a general bias (positive the first week of August and negative the rest of the month) but adequately reproduced the temporal variability in the SST field (correlation of 0.67). In the case of GLOBAL, the system could not benefit from data assimilation in this intricate coastal area with a low level of available observations: worse skill metrics were subsequently obtained with a correlation of 0.65 and a RMSD above 2 ^∘C.

Figure 13 Sequence of SAMPA daily surface circulation maps covering the period from 29 November to 5 December 2017. General map on the right and zoom over the Strait of Gibraltar on the left. An inflow reversal through the narrowest section of the Strait of Gibraltar is evidenced on 5 December, as a result of a change in the wind regime, from westerlies to easterlies.

Finally, outputs from SAMPA high-resolution coastal model were used to provide further insight into the entire AJ–WAG system and how diversity from the classical picture of the Alboran Sea surface circulation emerged from changes in the intensity and direction of the AJ. Although only one episode (corresponding to December 2017) is shown here (Fig. 13), the four events followed a similar scenario.

• i.

  Prelude: the AJ was observed flowing vigorously (with velocities clearly above 80 cm s⁻¹) into the Alboran Sea with a rather zonal direction (Fig. 13a), heading northeast later on, and surrounding and feeding the WAG, which appeared to be slightly detached from its traditional position in the Western Alboran Sea (Fig. 13b).

• ii.

  Onset: as westerly wind lost strength, the AJ speed became progressively weaker and tended to flow more southwardly, giving rise to a weakening and subsequent decoupling of the AJ–WAG system along with the reinforcement of an already existing small-scale coastal eddy that coexisted with the WAG (Fig. 13c and d). Circulation snapshots with three gyres (including the EAG, out of the pictures) have been previously reported in the literature (Flexas et al., 2006; Viúdez et al., 1998). This coastal eddy could be either cyclonic and confined northeast of Algeciras Bay (February 2017, not shown) or be anticyclonic, starting to grow, detach from the coast, and migrate eastwards as a result of both the change in AJ orientation and the WAG displacement (Fig. 13e and f). Meanwhile, the WAG presented different configurations: from an almost-symmetric aspect (August 2017, not shown) to a more elongated shape in the cross-shore direction (December 2017: Fig. 13f) or in the along-shore direction (March 2017, not shown).

• iii.

  Development: the AJ velocity reached a minimum (below 50 cm s⁻¹) associated with a sharp change in the predominant wind regime from westerlies to easterlies (Fig. 13g and h). A branch of the eddy, neighbouring the strait, was wind weakened and deflected from the main rotating pathway and started to flow westwards to the GIBST.

• iv.

  Full establishment of the inflow reversal: complete westward outflow from the Mediterranean Sea into the Atlantic Ocean through the narrowest section of GIBST, reaching a peak of velocity over Camarinal Sill (Fig. 13i). The migratory eddy and the WAG started merging into one single anticyclonic gyre (Fig. 13j).

• v.

  Epilogue: afterwards, in three of the cases the resettlement of predominant westerlies (Fig. 10m, n, p) favoured the return of the northeast-oriented Atlantic inflow and the consequent reactivation of the usual AJ–WAG system (not shown). By contrast, in the fourth episode (August 2017), summer easterly winds kept blowing moderately for two extra weeks (Fig. 10o) but were too weak to preserve the induced reversal, thus the Atlantic inflow reappeared again.

7 Conclusions

The current generation of ocean models have undergone meticulous tuning based on several decades of experience. The ever-increasing inventory of operational ocean forecasting systems provides society with a significant wealth of valuable information for high-risk decision-making and coastal management. Some of them are routinely operated in overlapping regions, offering the opportunity to compare and judge the strengths and weaknesses of each system, and eventually evaluate the added value of high-resolution coastal models with respect to coarser model solutions.

In this work, a multi-parameter model inter-comparison was conducted at the sea surface, ranging from global to local scales in a two-phase stepwise strategy. Firstly, a comparison of CMEMS products (GLOBAL and the nested IBI regional system) was performed against remote-sensed and in situ observations. In terms of temperature, results highlighted the overall benefits of both the GLOBAL direct data assimilation in open water and the increased horizontal resolution of IBI in coastal areas. GLOBAL seemed to replicate slightly better the SST field likely thanks to recent progresses in data assimilation schemes and to the growing wealth of available observational data. In this context, the assimilation of new types of observations (drifting buoys SST) should improve constrains in modelled variables and overcome the deficiencies of the background errors, in particular for extrapolated and/or poorly observed variables (Gasparin et al., 2018; Lellouche et al., 2018). In the open sea, IBI benefited indirectly from the data assimilation conducted in its parent system thanks to the spectral nudging technique, thereby allowing the regional and global model states to be reconciled with each other. IBI even outperformed GLOBAL locally in specific open-water zones such as the continental shelf break. IBI performance was also more accurate in those coastal regions characterized by a jagged coastline and a substantial slope bathymetry. As GLOBAL has a smoothed bathymetry and does not resolve many narrow features of the real sea floor, the depths where mixing takes place could be biased. Besides, those mixing processes acting at scales smaller than the grid cell size might substantially affect the resolved large-scale flow in the coarser GLOBAL system.

On the other hand, since GLOBAL is a detided model solution, tidally driven mixing could account for a portion of the discrepancies found between GLOBAL and satellite-derived SST estimations in energetic tidal areas such as the English Channel, the North Sea, and the Irish Sea. Whereas GLOBAL seemed to predict an over-stratification in shelf seas, IBI could better reproduce the vertical stratification and hence the SST field in the aforementioned subregions.

Complementarily, an isolated but rather illustrative example of the impact of impulsive-type river freshwater discharge on local surface circulation in NW Spain was provided. The increased horizontal resolution of IBI allowed a more accurate representation of horizontal salinity gradients, the horizontal extent of the plume, and the strength and position of the freshwater front, according to the results derived from the validation against in situ observations of SSS, SST, and currents provided by a moored buoy. Since both GLOBAL and IBI present 50 depth levels, similar vertical discretization and comparable climatological runoff forcing, the horizontal resolution is assumed to play a primary role when attempting to explain the differences encountered in model performance for this specific test case. Notwithstanding, the authors are fully aware that this single isolated example does not suffice and additional events over the entire IBI coastal domain should be examined in future work.

Finally, a 1-year (2017) multi-model inter-comparison exercise was performed in the Strait of Gibraltar between GLOBAL, IBI, and SAMPA coastal systems in order to elucidate the accuracy of each OOFS to characterize the AJ dynamic. A quantitative comparison against hourly HFR estimations highlighted both the steady improvement in AJ representation when moving from global to coastal scales through a multi-nesting model approach and also the relevance of a variety of factors at local scales.

• i.

  A sufficiently detailed representation of bathymetric features: the very high horizontal resolution of SAMPA (∼400 m) and, consequently, the tailored bathymetry employed in order to capture small-scale ocean processes and resolve sharp topographic details.

• ii.

  A better representation of air–sea interactions: the adequate refinement of the spatio-temporal resolution of the atmospheric forcing used in SAMPA, especially in a complex coastal region where topographical steering further impacts on flows.

• iii.

  The nesting strategy implemented in SAMPA, where the coastal model is forced not only by daily fields provided by the regional IBI system but also by outputs from Mog2D and NIVMAR models. The additional inclusion of these accurate tidal and meteorologically driven barotropic velocities, respectively, allowed a detailed characterization of the variability in the surface Atlantic inflow, including persistent full reversal episodes. Although the matching between HFR observations and SAMPA outputs is mainly found in two of the four reversal events detected, this result demonstrates its added value as modelling tool towards the comprehension of such singular oceanographic events. A detailed characterization of this phenomenon is relevant from diverse aspects, encompassing search and rescue operations, the management of accidental marine pollution episodes, or safe ship routing.

Finally, SAMPA coastal model outputs were analysed in order to put in a broader perspective the context of the onset, development, and end of such flow reversal and its impact on the AJ–WAG coupled system. The synergistic approach based on the integration of HFR observing network and SAMPA predictive model has proved to be valid to comprehensively characterize the highly dynamic coastal circulation in the GIBST and the aforementioned episodic full reversals of the surface inflow. In this context, data assimilation (DA) would provide the integrative framework for maximizing the joint utility of HFR-derived observations and coastal circulation models. A DA scheme could be incorporated in future operational versions of SAMPA in order to improve its predictive skills, since similar initiatives are currently ongoing with positive results (Vandenbulcke et al., 2017; Stanev et al., 2015). Although DA is a powerful technique, advances in coastal ocean modelling should also encompass an improved understanding of high-frequency small-scale physical processes; the accurate model parameterization of the effects triggered by such sub-grid phenomena; and the integration of air–sea, wave–current, and biophysical interactions by means of coupled forecasting systems (Wilkin et al., 2017). Complementarily, as part of possible next improvements of SAMPA system, we will explore the possibility of a more direct nesting strategy into IBI since the current operational version of IBI includes the delivery of 3-D hourly outputs.

Future efforts are planned to improve CMEMS global and regional OOFSs in several aspects already addressed in the present work. While the GLOBAL system will be evolved towards a 1∕36^∘ model application, a substantial refinement will be accomplished for the regional IBI system in both vertical and horizontal resolutions: from 50 to 75 depth layers and from 1∕36 to 1∕108^∘ (∼1 km), respectively. Whereas the first feature will be incorporated during CMEMS Phase-2 (2018–2021), the second milestone will be achieved in the frame of the IMMERSE (Improving Models for Marine EnviRonment SErvices) H2020 project (https://cordis.europa.eu/project/rcn/218810/factsheet/en, last access: 1 July 2019) and is expected to positively impact on a more accurate representation of coastal processes: submesoscale shelf break exchanges and connectivity, fronts, river plumes, or topographic controls on circulation. In addition, a more detailed bathymetry is expected to be introduced in future operational versions of IBI in order to better resolve those regions with complex coastline and intricate bottom topography. Other factors that could be potentially improved but still deserve further analysis are the air–sea and the land–sea interactions, i.e. the meteorological and riverine forcings. With regards to the former, a more skilful atmospheric forecast model with a higher spatio-temporal resolution (i.e. hourly prediction over a more refined grid) could aid to better represent the coastal circulation by a more accurate discrimination of the topographic structures and the replication of the inertial oscillations and mesoscale processes. On the other hand, each main river basin hydrology should be more accurately represented with daily-updated outputs from tailored hydrological models. Finally, refined mixing schemes might also produce notable improvement in the representation of water masses, resulting in a substantial reduction of temperature and salinity bias relative to model solution.