Multi-platform model assessment in the Western Mediterranean Sea: impact of downscaling on the surface circulation and mesoscale activity

In numerical ocean modeling, dynamical downscaling is the approach consisting in generating high-resolution regional simulations exploiting the information from coarser resolution models for initial and boundary conditions. Here we evaluate the impacts of downscaling the 1/16o (~ 6–7 km) CMEMS Mediterranean reanalysis model solution into a high-resolution 2-km free-run simulation over the Western Mediterranean basin, focusing on the surface circulation and mesoscale activity. Multi-platform observations from satellite-borne altimeters, high-frequency radar, fixed moorings, and gliders are used for this evaluation, providing insights into the variability from basin to coastal scales. Results show that the downscaling leads to an improvement of the time-averaged surface circulation, especially in the topographically complex area of the Balearic Sea. In particular, the path of the Balearic current is improved in the high-resolution model, also positively affecting transports through the Ibiza Channel. While the high-resolution model produces a similar number of large eddies as CMEMS Med Rea and altimetry, it generates a much larger number of small-scale eddies. Looking into the variability, in the absence of data assimilation, the high-resolution model is not able to properly reproduce the observed phases of mesoscale structures, especially in the southern part of the domain. This negatively affects the representation of the variability of the surface currents interacting with these eddies, highlighting the importance of data assimilation in the high-resolution ocean model in this region to constrain the evolution of these mesoscale structures.


Introduction
Monitoring and modeling the ocean and coastal areas are crucial for both society and science (Visbeck 2018), with uses that include safety-at-sea (e.g., ocean-rescue, storm surges), economics (e.g., sustainable fisheries, offshore energy designs, support to navigation), and the environment (e.g., oil spill tracking, plastic debris control). Numerical models are essential tools since they cover a wide range of spatial and temporal scales and are complementary to ocean observations, which are generally sparse. The spatial scales of the monitoring and modeling system must match the end-user requirements as much as possible. In particular, high-resolution models are required for regional and local applications and for the understanding of processes over a wide range of scales. Short spatial scales are especially important in coastal areas, as well as in straits, channels, and transition zones connecting coastal and open seas, where they have the potential to provide refined representations of transports and exchanges of heat, fresh water, and biogeochemical tracers.
In order to reach the required small scales, the "downscaling" approach is employed to generate high-resolution regional simulations. The downscaling techniques, also widely used in atmospheric sciences (Navon, 2009), use initial and boundary conditions from a lower resolution model previously implemented in the same domain ("parent model") to drive higher resolution regional simulations (e.g., Mason et al. 2010), usually forced by higher-resolution atmospheric forcings. The downscaling approach permits the resolution of small-scale processes and representation of coastal processes of interest (e.g., topographically constrained coastal circulation, sea breezes), while exploiting the large scale information provided by the parent model (Capó et al. 2018). The question arises whether this increase in resolution really leads to a better realism of the simulations, and to what extent. This question is addressed in this paper by considering a high-resolution multi-year freerun simulation in the Western Mediterranean Sea, by means of a multi-platform evaluation approach.
In the Western Mediterranean Sea (WMS) the internal Rossby radius is around 10-15 km (Font et al. 1988;Robinson et al. 2001;Escudier 2015), smaller than its values in the Atlantic and Pacific oceans at the same latitudes. Processes associated with small spatial scales, such as the mesoscale (around 10-100 km) and the submesoscale (around 1-10 km), account here for an important percentage of the total energy and dominate most of the movements of the sea (Robinson et al. 2001). The current and growing interest in understanding these small-scale features and their relevance in ocean energetics converts the Mediterranean Sea into a reduced-scale ocean laboratory for the examination of processes of global importance (e.g., the CALYPSO Research Initiative of the Office of Naval Research (CALYPSO-ONR, 2017)). The implementation of high-resolution models that resolve the mesoscale and permit some submesoscale variability is therefore particularly interesting in this area.
In this context, the Balearic Islands Coastal Observing and Forecasting System (SOCIB, Tintoré et al. 2013) has developed a 2-km-resolution free-run hindcast simulation based on the high-resolution Western Mediterranean Operational System (WMOP, Juza et al. 2016;Mourre et al. 2018). This system is nested in the larger scale CMEMS Med Rea reanalysis simulation (Simoncelli et al. 2014) through a downscaling approach. The WMOP high-resolution hindcast simulation under consideration does not include data assimilation. This kind of simulations is suitable for process and variability studies (e.g., formation and propagation of eddies, eddy statistics, energy budgets, water mass formation) since they provide continuous dynamics, contrary to data-assimilating simulations, whose fields are perturbed at regular time intervals according to the information provided by observations. Moreover, they may also serve as a base for efficient data assimilation if they are proved unbiased, sufficiently realistic, and able to statistically represent the main processes detected in observations. However, due to the lack of constraint by observations, free-running simulations are also expected to be barely able to represent mesoscale structures at their proper times and positions.
The general objective of this paper is to perform a multiplatform assessment of both the parent reanalysis model and the downscaled free-run simulation spanning the period 2009-2015, to report the effects of the downscaling procedure on the time-averaged surface circulation and its temporal variability, as well as the mesoscale activity over the WMS. Multi-platform observations from high-frequency radar (HFR), satellite-borne altimeters, fixed moorings, and gliders are used for this purpose.
Specifically, this paper aims (1) to evaluate to what extent the high-resolution model modifies the time-averaged circulation, both over the entire basin and in the Balearic Channels and coastal areas monitored by HF radars, (2) to report the effects of downscaling on the variability of the surface slope currents along the Iberian margin under the influence of evolving mesoscale structures, and (3) to estimate the impact of the higher resolution on the mesoscale activity, evaluating the abundance and size of eddies over the whole domain, and their propagation velocities around the Algerian basin.
The paper is organized as follows. Section 2 presents a brief description of the study area, the models, and the observation datasets; section 3 reports the effects of downscaling on the time-averaged surface circulation, considering both basin and coastal scales, the variability of the slope currents, and the mesoscale eddy activity. Finally, section 4 discusses the results and presents the main conclusions and future work.

Study area
The Western Mediterranean Sea (WMS, here considered between 35-44.5°N and 6°W-9°E, so excluding the Tyrrhenian Sea) is a semi-enclosed sea connected to the Atlantic Ocean through the Gibraltar Strait ( Figure 1). The surface circulation in the WMS is mainly cyclonic around the basin. The northern part of the domain is dominated by the Northern Current (NC), which flows westward along the Gulf of Lions and the Iberian coast until the Balearic Sea. Approaching the Ibiza Channel (IC), the NC splits into two branches (Font et al. 1988), one deviating eastwards along the Balearic northern shore and forming the Balearic Current (BC, Violette et al. 1990;Ruiz et al. 2009;Mason and Pascual, 2013) and the other flowing southward through the IC (Pinot et al. 1994;Pinot et al. 1995). In the southern WMS, the Alboran Sea and the Algerian basin are characterized by the presence of permanent and semi-permanent eddies (Tintore et al. 1988;Testor et al. 2005;Renault et al. 2012;Escudier et al. 2016a) and by the Algerian Current, flowing along the African coast. The circulation in this sub-basin is mainly driven by the Atlantic Jet (AJ, Parrilla and Kinder, 1987;Tintoré et al. 1991;Viúdez and Haney, 1997), which feeds the permanent Western Anticyclone Gyre (WAG, Viúdez et al. 1996;Baldacci et al. 2001;Flexas et al. 2006) located between the Gibraltar strait and the 3.5°W longitude (Juza et al. 2015), and the less intense and semi-permanent Eastern Anticyclonic Gyre (EAG) located between 2.5°and 1.5°W . The IC was shown to be an important circulation choke point with highly variable meridional water mass transports (Heslop et al. 2012;Lafuente et al. 1995;Pinot et al. 2002). Important exchanges between the fresher Atlantic waters and the saltier Mediterranean waters take place in this channel, impacting the transports of heat, salt, and nutrients and with relevant impacts on regional ecosystems (Balbín et al. 2014).

Numerical models
In this work, we use a free-run hindcast simulation spanning the period 2009-2015 and based on the WMOP forecasting system configuration (Juza et al. 2016;Mourre et al. 2018). WMOP is a regional configuration of the Regional Ocean Modelling System (ROMS, Shchepetkin and McWilliams, 2005) with a horizontal resolution around 2 km (more specifically, the resolution ranges from 2.16 to 1.88 km from South to North). It uses 32 stretched sigma-levels resulting in a vertical resolution varying from 1 to 2 m at the surface, 30-40 m at 200 m depth, and around 250 m at 1000 m depth. The vertical terrain-following sigma levels are controlled by the following parameters: transformation equation Vtransform = 2 and stretching function Vstretching = 4 with θs = 7 (stretching coordinate for the surface boundary layer), θb = 0.2 (stretching in the bottom boundary layer), and hc = 60 m (critical depth). The bottom topography is derived from a 1min resolution database (Smith and Sandwell, 1997). The WMOP free-run hindcast simulation used in this study is nested in the larger-scale reanalysis CMEMS Med Rea (MEDSEA_REANALYSIS_PHYS_006_004, Simoncelli et al. 2014). The WMOP simulation uses surface forcing from the high-resolution (3 h-5 km) HIRLAM model from the Spanish Meteorological Agency AEMET (Undén et al. 2002). The simulation includes observed daily river discharges from the six major rivers of the domain (Var, Rhône, Aude, Hérault, Ebro, and Júcar) over the period 2009-2015 provided by the French HYDRO database and the Spanish hydrographic confederations of Ebro and Júcar rivers. WMOP uses the generic model of the two-equation Generic Length Scale turbulence closure scheme (Umlauf and Burchard, 2003, with parameters p = 2.0; m = 1.0; n = − 0.67 as in line 1 of their Table 7) for turbulent vertical mixing of momentum and tracers. The initial and open boundary conditions are imported from daily outputs of the CMEMS Med Rea reanalysis simulation and implemented in ROMS using Flather conditions (Flather, 1976) on the 2-D momentum equations and Chapman conditions (Chapman, 1985) for the surface elevation. Mixed radiation-nudging boundary conditions (Marchesiello et al. 2001) are used for 3D momentum equations. A particular treatment including the alignment of the CMEMS model and ROMS bathymetries and a correction of interpolated velocities is applied at the Strait of Gibraltar to ensure that the boundary forcing field properly represents the original inflow and outflow transports of the large scale model across the strait. More concretely, after interpolation of the daily CMEMS Med Rea velocities onto the WMOP grid along the open boundary section in the Strait of Gibraltar, positive (resp. negative) zonal velocities are corrected by addition/ subtraction of a depth-averaged velocity over the inflow (resp. outflow) section so that the corresponding transports after interpolation on the WMOP grid are equal to the transports provided by the CMEMS Med Rea model.
The larger-scale reanalysis CMEMS Med Rea model has a 1/16°horizontal resolution (around 6-7 km) and 72 unevenly spaced vertical z-levels (Oddo et al. 2009) resulting in a surface vertical resolution around 3 m. The model uses vertical partial cells to fit the shape of the ocean bottom. The bottom topography is derived from a standard 1-min data set of the U.S Navy Digital Bathymetric Data Base (DBDB1). CMEMS Med Rea is based on the NEMO model (Madec 2008) with surface forcing from ERA-Interim (ECMWF; 6 h, 70 km, Dee et al. 2011). It uses a daily variational data assimilation scheme (OceanVAR) developed by Dobricic and Pinardi (2008). The assimilated data include the following: sea level anomaly, sea surface temperature, in situ temperature profiles by VOS XBTs (Voluntary Observing Ship-Expendable Bathythermograph), and in situ temperature and salinity profiles from Argo floats and CTDs (conductivitytemperature-depth). This model uses a mixed up-stream/ MUSCL parameterization for the advection scheme for active tracers. The vertical diffusion and viscosity terms are function of the Richardson number as parameterized by Pacanowski and Philander (1981).

Satellite altimetry
Daily delayed time gridded datasets of sea level anomaly (SLA) with 1/8°(around 12 km) spatial resolution over the Mediterranean region are distributed by CMEMS and referenced as, SEALEVEL_MED_PHY_L4_REP_OBSERVATION-S_008_051. In this work, we use data extracted from 2009 to 2015 in an area extending from 6°W-9°E to 35°N-44°N. This product is created by the SL-TAC multimission altimeter data processing system. It processes data from all altimeter missions: HY2, Saral/AltiKa, Cryosat-2, Jason-2, Jason-1, T/P, ENVISAT, GFO, and ERS1/2. All altimeter fields are interpolated at crossover locations and dates. Data are then cross-validated, filtered for residual noise and small-scale signals, and sub-sampled. Finally, an optimal interpolation is applied in order to compute gridded SLA by merging data from all of the satellites. The mean dynamic topography (MDT; Rio et al. 2014) is added to the SLA to compute the absolute dynamic topography (ADT).

Fixed moorings
In situ observations of surface currents are extracted from six oceanographic buoys installed and maintained by Puertos del Estado (PdE; http://puertos.es) and located along the Iberian coast ( Fig. 1) from January 2008 to December 2016 with a 1-h sampling frequency. The exception is the IC mooring, installed and maintained by SOCIB, available since 2013 (see Table 1) and with a 10-min sampling frequency. The mooring depths and the availability period of each mooring are displayed in Table 1.

Gliders
Several glider missions (31) were undertaken in the IC between 2011 and 2015 by SOCIB using Slocum deep gliders (Heslop et al. 2012). The gliders performed saw-tooth patterns through the water column between depths of 20 and 950 m, while navigating a standard transect line at 39°N approximately perpendicular to the main currents (subset in Fig. 1). Each latitudinal transect (80 km) took an average of 3 days to complete. The thermal wind equation applied to the zonal hydrographic section allows us to infer the meridional geostrophic velocities along the section performed at 39°N following the method of Pinot et al. (1995b). It assumes a zero velocity reference depth of 800 m or the seafloor where the profile depth is shallower than 800 m (Heslop 2015). Current time-series were then extracted at 47 points over the glider transects.

High-frequency radar
The high-frequency radar (HFR) provides high-resolution synoptic observations of sea surface currents, with a limited geographical coverage but high spatiotemporal resolution, similar to that of the high-resolution model. Three different operational HFR systems monitor strategic coastal areas of the WMS: the Ebro delta region, IC, and Strait of Gibraltar. The Ebro delta HFR system was deployed in December 2013 and is currently maintained and operated by PdE. The HFR network consists of three Coastal Ocean Dynamics Applications Radar (CODAR) SeaSonde systems (Paduan and Rosenfeld, 1996), transmitting at a central frequency of 13.5 MHz that provide hourly radial measurements at 3 km spatial resolution.
In the eastern part of the IC, the monitoring system has been working since June 2012 and is based on two CODAR monostatic medium-range operated by SOCIB , also transmitting at a central frequency of 13.5 MHz. The area extends from the coast to 60 km offshore. The spatial resolution is 3 km for zonal and meridional velocity components. In the Strait of Gibraltar, the HFR system comprises three shore-based CODAR SeaSondes also owned and operated by PdE, transmitting at a central frequency of 26.28 MHz. The data collection started in May 2011 with two antennas; it was completed with a third one in 2012 to gain accuracy and insight over the entire radar coverage. The HFR hourly derived surface current maps are provided at 1 km spatial resolution (Mader et al. 2016).

Basin-scale
The 2009-2015 time-averaged geostrophic surface currents from altimetry and both models are represented in Fig. 2. The spatial patterns show good consistency between them, with velocities of the same order of magnitude. In general terms, both models represent reasonably the mean flows (NC, Algerian current, Balearic current, exit flow south of Sardinia) and the permanent WAG structures in the Alboran Sea. Specifically, both models describe quite well the NC along the French shelf break until reaching the latitude of Barcelona around 41.5 o N, with two acceleration zones at both ends of the Gulf of Lions, although the amplification produced at 7°E is a bit underestimated in WMOP (maximum speeds of 33, 25, and 36 cm s −1 , for altimetry, WMOP and CMEMS Med Rea, resp.). Further south, the two models represent differently the NC retroflection flow pattern. In the CMEMS Med Rea model, this pattern is marked by a strong time-averaged current with a detachment from the Iberian shelf break which occurs around 2 o E, before reaching the Tarragona mooring location (northernmost black point in Fig. 2). This retroflection is too intense and too far to the east compared with altimetry and to the WMOP model. The WMOP model correctly represents this retroflection, with a moderate southward flow centered around 1.5°E before flowing northeastward along the Balearic shelf break, in good agreement with altimeter estimates. This circulation bias in CMEMS Med Rea was already pointed out in Juza et al. (2016) and in Pinardi et al. (2015). The early NC detachment in CMEMS Med Rea in turn affects the circulation in the Ibiza Channel, where the parent model describes a homogeneous northward flow, especially intense on both sides of the channel, in contrast to altimetry and WMOP where this northward flow is confined to the eastern part of the channel. Compared with CMEMS Med Rea, WMOP also provides an improved representation of the mean currents along the northeastern coast of the Menorca Island, representing southeastward flows rounding the island, which are in agreement with altimeter estimates. In the Balearic Sea (39.5-41°N; 1-5°E), the improvement provided by WMOP in comparison with the parent model is around 2 cm s −1 for both velocity components (u, v). The root mean square difference (RMSD) with respect to altimeter estimates for WMOP and CMEMS-Med Rea is respectively: for u 3.6 vs 5.1 cm s −1 ; for v 3.8 vs 5.6 cm s −1 .
Regarding the southern part of the domain, both models represent the WAG, with a slightly better representation in WMOP than in CMEMS Med Rea. The latter is not able to reproduce the high velocities reported by WMOP and altimeter observations in its southern part. Indeed, CMEMS Med Rea overestimates the maximum velocities in the northern part of the gyre (55, 69, and 90 cm s −1 , for altimetry, WMOP, and CMEMS Med Rea, resp.). Regarding the EAG, the two models significantly differ in their time-averaged representation. While the EAG dimensions are underestimated in CMEMS Med Rea, they are significantly overestimated in WMOP. Overall, both models present significant errors in the Alboran sea (35-37°N; 3-1°W), these being slightly larger in WMOP than in the parent model (RMSD: for u 12 vs 10 cm s −1 ; for v 9 vs 8 cm s −1 ) due to the overestimation of the EAG in WMOP.
The Algerian current (36-37.5°N; 0-6.5°E) has similar mean characteristics in both models and in altimetry, with maximum values of 49, 43, and 36 cm s −1 , for altimetry, WMOP and CMEMS Med Rea, resp. However, it penetrates slightly further east in WMOP than in CMEMS Med Rea, in better agreement with the altimeter estimates. This is reflected in the RMSD computed in this area, especially considering the zonal component (RMSD: for u 4.5 vs 6.3 cm s −1 ; for v 3.4 vs 3.7 cm s −1 ).
The models also depict a south-eastward time-averaged current flowing from the Balearic Sea (5°E-41°N) to the southern part of Sardinia (6°E-40°N). This current, aligned along the mistral wind corridor, is not as pronounced in altimetry. In the southeastern part of the WMOP domain, both models have slightly different circulation patterns. In particular, the mean WMOP surface velocities show a marked anticyclonic eddy centered around 6.5°E-40°N and a strong coastal detachment of the Algerian current around 6°E, which are not present in the altimeter estimates. At the same time, WMOP properly represents the northward flow off Sardinia Island north of 40°N, which is not the case for CMEMS Med Rea. CMEMS Med Rea properly represents the southward circulation pattern along the southern part of the Sardinian coast but with an overestimation of the magnitude of the mean flow. WMOP also represents this southward coastal flow but with a reduced width. The exit flows south of Sardinia and north of Corsica are represented in both models while the resulting RMSD found in this area (37-40°N; 7-9°E) are very slightly lower in CMEMS Med Rea (RMSD: for u 7 vs 6.1 cm s −1 ; for v 6 vs 5.7 cm s −1 ), the maximum values are better described by CMEMS Med Rea (25, 26, and 32 cm s −1 , for altimetry, WMOP and CMEMS Med Rea, resp.) While the temporal average surface currents represented in Fig. 2 allows the characterization of the permanent structures and the mean flows, the mean kinetic energy (KE) and the eddy kinetic energy (EKE) provide complementary diagnostics accounting for the temporal variability along the period of the simulation. Mourre et al. 2018 compared the mean EKE deduced from altimetry and from the WMOP simulation when considering surface total velocities and geostrophic estimates with and without spatiotemporal filtering. They showed that the EKE of the surface geostrophic velocities of the WMOP model without filtering was around three times larger than the EKE derived from altimeter estimates. We extend here these diagnostics to CMEMS Med Rea. The comparison of the averaged KE and EKE values shows that the higher-resolution model represents a much more energetic dynamics than its parent model and altimetry (EKE WMOP 202.1 cm 2 s −2 ; EKE CMEMS-Med-Rea 99.6 cm 2 s −2 ; EKE ALT 83.7 cm 2 s −2 ; KE WMOP 258 cm 2 s −2 ; KE CMEMS-Med-Rea 151.7 cm 2 s −2 , and KE ALT 152.8 cm 2 s −2 ).

Coastal-scale
Surface circulation in HF radar coastal areas In this section, the performance of both models is assessed in terms of surface currents in three coastal areas monitored with HFR: the Strait of Gibraltar, the Ibiza Channel, and the Ebro delta region. Figure 3 shows the time-averaged surface currents obtained from HFR, WMOP, and CMEMS Med Rea data in the three areas.
In the Gibraltar Strait area (Fig. 3a), which is also the western boundary of the WMOP domain, all datasets describe the Atlantic Jet (AJ), although with some differences. WMOP provides mean velocities with higher magnitude than CMEMS Med Rea. These velocities are in better agreement with the HF radar estimates in the strait itself but are overestimated once entering the Mediterranean Sea. Moreover, the average AJ direction in WMOP is oriented too far to the north, in contrast to CMEMS Med Rea which describes an average direction to the east in agreement with the HF radar velocities. The RMSD with respect to HFR data shows a large error for both models, this error being significantly higher in WMOP than in its parent model, especially in the meridional direction (RMSD: for u 21.3 vs 16.1 cm s −1 ; for v: 20.1 vs 13.8 cm s −1 for WMOP and CMEMS Med Rea, resp.) In the Ibiza Channel (Fig. 3b), WMOP represents the average meridional surface flow reversal (northwards in the eastern part and southwards in the western part) described by HFR, although slightly displaced northwestward and with higher intensities. It improves the circulation compared with CMEMS Med Rea, which shows a homogeneous northward flow over the whole HF radar area, confirming the indications provided by altimetry in the previous section. This remains reflected in the RMSD computed the meridional component, which results to be 1.4 cm lower for WMOP than for CMEMS Med Rea (RMSD: for v 4.9 vs 6.3 cm s −1 ). Regarding the zonal component, the error obtained by the parent model is lower for CMEMS Med Rea due to an overestimation of the flow in WMOP in the southern part (RMSD for u: 4.7 vs 2.1 cm s −1 ).
In the Ebro delta region (Fig. 3c), WMOP improves the average position of the NC compared to its parent model, which locates it outside of the HFR domain due to the retroflection issue mentioned in the previous section. Moreover, WMOP qualitatively reproduces the small-scale coastal flow intensification at the mouth of the Ebro River, which is not achieved in the coarser resolution CMEMS Med Rea model. This improvement is depicted in the zonal and meridional RMSD of the two models (RMSD: for u 2.5 vs 3.3 cm s −1 ; for v: 6.6 vs 7.6 cm s −1 ). Despite this improvement, WMOP clearly underestimates the intensity of the NC as measured by the Ebro delta HFR. This behavior was also observed by Lorente et al. (2016) using the Iberia-Biscay-Ireland (IBI) regional circulation model (Sotillo et al. 2015), a NEMO-based model from PdE, which has a similar spatial resolution as WMOP. These authors relate these shortcomings to the low resolution and realism of the atmospheric forcing used in their model, which is the same as in WMOP. However, comparisons between the winds extracted from the moored buoys in Tarragona and Ibiza Channel and the wind used as forcing in WMOP do not support this hypothesis (not shown). The Ebro delta region is characterized by strong instabilities produced by the river discharges and by significant temporal and spatial variability induced by the shelf/slope fronts (Wang et al. 1988;Pinot et al. 1994), which makes the realistic modelling of this area challenging.
Surface meridional flows inferred from glider data The routine glider measurements carried out by SOCIB since 2011 provide complementary information on the meridional transports in the IC. The time-averaged meridional surface geostrophic velocities inferred from the glider are compared here to HFR measurements, altimetry data and both model estimates along the glider section (Fig. 4). Notice that the glider velocities are not directly observed but are deduced from the application of the thermal wind equation to the hydrographic section, which is dependent on the choice of a level of no motion (800 m in this study). This process introduces some uncertainty in the computed velocities (Heslop 2015). Surface velocities are only averaged over the glider sampling periods. Glider observations provide a pattern characterized by a central southward flow in the deep part of the channel, a secondary southward flow along the 500-m isobaths on the western side, and two northward flows on the shallower parts of the section on each side of the channel. Both WMOP and altimetry show the presence of both northward and southward flows, but with a spatial pattern differing from the glider estimates, northward velocities being found on the eastern part of the channel and southward velocities on the western side. This comparison illustrates the limitation of conventional altimetry in these kinds of coastal areas. WMOP improves the timeaveraged velocities compared with the CMEMS Med Rea model, which shows a northward flow over the entire channel section. This bias is very probably linked to the Balearic Current retroflection issue which was highlighted previously. Yet, the zonal modulation of the flow with more intense northward flows on the sides of the channel is captured by CMEMS Med Rea. Observations from gliders and HFR describe very similar spatial patterns in the eastern part of the channel, with a flow reversal from southwards in the deep part of the channel, to northwards when approaching the coast. Glider velocity estimates are found to have a slightly larger magnitude than HFR velocities despite the fact that we are comparing geostrophic velocities from the glider versus total velocities from HFR.

Variability of surface currents at fixed mooring locations
3.2.1 Main axes of variability Figure 5 illustrates variance ellipses computed at mooring locations. The quantities represented in this figure are close, but not exactly the same as the classical variance ellipses presented in e.g., Morrow et al. (1994). In both cases, the represented ellipses are oriented along the axis of largest variance (major axis). While in the classical variance ellipses, the length of the two axes scales with the standard deviation of the velocities along these axes, in our case it is proportional to the percentage of explained variability. The ellipses in Fig. 5 represent (1) the direction of highest variability of the flow denoted by the major axis of the ellipse, (2) the percentage of explained variability along this direction, and (3) the percentage of explained variability along the minor axis providing insights into the directional variability of the flow. The length of the ellipse axes are proportional to the percentage of variability explained along these axis. If σ x is the velocity variance along the major axis and σ y the variance along the minor axis, these percentages are calculated as: (σ 2 x = σ 2 x þ σ 2 y Þ⋅100 and σ 2 y = σ 2 x þ σ 2 y ⋅100 . Ellipses elongated in the direction of the major (resp. minor) axis denote relatively low (resp. high) directional variability. Results show good agreement between both models and moorings in TAR and DRA, both in terms of axis and directional variability, in contrast to altimetry, which does not get the right direction of the main axis and also underestimates the directional variability. In general terms, the obtained directional variability using altimetry is lower than that reported by moorings, which suggests that altimetry is not capturing all the variability existing in the region. This fact is probably related to coastal restrictions and the limitations of altimetry when approaching the coast (Cipollini et al. 2010;Pascual et al. 2015;Vignudelli et al. 2011).
The improvement produced by WMOP is noteworthy at the VAL mooring, where neither the CMEMS Med Rea nor altimetry are able to capture the variance shown by the moorings. At that location, surface velocities from WMOP are in very good agreement with that of the mooring in terms of the main axis of variability and its directions.
In the Ibiza Channel (IC), the variance showed by the mooring reveals the high directional variability of this area, already described in some other works (Heslop et al. 2012;Lafuente et al. 1995;Mason and Pascual 2013;Pinot et al. 2002). This variability is also represented by the other datasets, although underestimated, especially by altimetry. While the direction of the main axis of variability provided by both models is in good agreement with the mooring and oriented along the meridional direction approximately, it has an artificial zonal deviation in altimetry.
At the PAL mooring location, both WMOP, CMEMS Med Rea, and altimetry provide quite accurate estimates of the main axis of variability and similar directional variability.
In the southernmost position (GAT), all the datasets show the same zonal direction for the main axis of variability, linked to the variability of the EAG. However, none of the models is able to capture the high variability reported by mooring and altimetry observations.

Inter-annual changes
The inter-annual changes of the surface circulation in the four available datasets are evaluated at the locations where the percentage of gaps for the mooring data was lower than 10%, namely at TAR, PAL, and GAT moorings (Fig. 6). Surface current time series have been projected into their main axis (see Fig. 5) and then high-frequency filtered using a moving average of 182 days to smooth out short-term fluctuations and focus on the seasonal and inter-annual variability. Positive/negative values denote northeastward/southwestward flows, respectively. Figure 6 confirms the systematic southward flow underestimation from the CMEMS Med Rea at TAR location (previously reported in section 3.1) during all the study period. This effect is corrected in WMOP, which shows better agreement with the mooring observations and with altimetry. On the contrary, at the southern positions (PAL and GAT), CMEMS Med Rea properly represents the observed inter-annual variability, which is not the case for WMOP, in particular during the years 2012 and 2013 (northward flow not captured at PAL and eastward flow overestimated at GAT). The variability reported by the altimetry is in very good agreement with mooring locations, especially at these two southern positions.
To better understand the WMOP errors highlighted at PAL and GAT, Fig. 7 illustrates the surface currents averaged from July to September of 2012 from WMOP and altimetry. The anticyclonic eddy depicted by altimetry in the area of PAL mooring (Fig. 7b), which produces intense positive (northeastward) velocities at the mooring location (as seen in Fig. 4 Time-averaged meridional (N-S) surface geostrophic flows along the Ibiza Channel glider section, from gliders (GLI; yellow), altimetry (ALT; green), WMOP (red), and CMEMS Med Rea (black) data. A surface current map from HFR averaged over the same period is also represented (HF; blue). The isobaths of 500 and 1000 m depth are also shown. Surface glider velocities are relative to 800 m depth (or the bottom if it is shallower) Fig. 6, PAL), is not present in WMOP (Fig. 7a). WMOP represents an eddy with similar characteristics, but it is located further northeast. As a consequence, the cyclonic eddy between the EAG (in GAT) and this second eddy is much more elongated in WMOP compared with altimetry. Low velocities are thus found in WMOP when focusing on the PAL mooring location. In the southernmost position (GAT), WMOP (Fig.  7a) is able to generate and reproduce the EAG as showed by altimetry (Fig. 7b). However, WMOP generates the gyre slightly bigger than altimetry and locates it too far north, directly affecting the comparison at this mooring location by producing locally overestimated high zonal velocities (as seen in Fig. 6, GAT). This analysis highlights how a small spatial mismatch in the model representation of eddies with respect to observations can lead to very large errors in point-to point velocity comparisons at specific locations.
This fact highlights the importance of the mesoscale structures in this area for the proper representation of the variability of surface currents. Data assimilation, which is applied in CMEMS Med Rea but not in the WMOP simulation, is therefore also necessary in the high-resolution model to properly locate mesoscale structures in space and time and so improve the comparisons with the mooring velocity observations.

Statistical abundance and size
The recent development of oceanic eddy tracking algorithms allows the automatic identification and tracking of eddies from SLA altimeter maps (Isern-Fontanet et al. 2003;Morrow et al. 2004;Chelton et al. 2007;Chaigneau et al. 2008). These are useful not only to improve our knowledge of eddy properties and their variability but also as a validation tool when comparing data from models and observations. In this work, we have used an automated eddy tracker, the so-called py-eddy-tracker (Mason et al. 2014), which detects closed contours in SLA maps to determine the presence, position, and size of eddies. This code has been applied to both models and altimetry databases from 2009 to 2015 to assess the effects of model downscaling on the mesoscale eddy statistics. As a first validation step, the eddy tracker algorithm was applied to spatially filtered SLA fields Fig. 5 Main axis of variability of the surface ocean currents at the fixed mooring locations (as shown in Fig. 1), obtained from WMOP (red), CMEMS Med Rea (black), altimetry (green), and moorings (blue). The axis represented within the ellipse denotes the direction of highest variability of the flow and its length is proportional to the percentage of explained variability. The minor axis represents the directional variability of the main axis, and it is also proportional to its length from both models and altimetry. The model fields were previously interpolated on the altimeter spatial grid of~12 km resolution prior to apply the algorithm. A 48 km moving average was also applied to match the effective resolution of altimetry (~40-50 km; Chelton et al. 2011;Escudier et al. 2016b) and to remove the effects of small-scale structures present in the models but not represented in altimetry. In the second part of the section, the same comparison was made using as inputs of the eddy tracker the original SLA fields from the three datasets (without filtering and considering the original spatial resolutions of each source).
When using the same spatial resolution in all datasets (Fig. 8, continuous lines), both models are able to generate a similar quantity of eddies over 7 years (1225 and 1253 eddies for WMOP and CMEMS Med Rea) of analogous sizes (between 20 and 30 km) as detected in the altimeter observations (1254 eddies). This fact demonstrates that the high-resolution model does not statistically generate more large eddies (those that can be observed by altimetry) than its parent model. In terms of the eddy daily-density, results are also similar for the three datasets (8 eddies per day for altimetry and 7 eddies per day for both models) with an averaged radius of 38 km for models and of 40 km using the altimeter observations. The comparison of the eddy tracker results when applied to the original model fields (Fig. 8, dashed lines) illustrates that the higher spatial resolution of WMOP leads to a larger number of detected eddies (WMOP 37 eddies/day, CMEMS Med Rea 20 eddies/day and altimetry 15 eddies/day) due to its capacity to represent eddies of smaller size (mean radius for WMOP 16 km, CMEMS Med Rea 26 km, and altimetry 33 km). These results are in agreement with Escudier et al. (2016b), who reported mean radii from 30 km (altimetry) to 25 km (~3 km resolution model) with densities from 15 eddies/day (altimetry) to 34 eddies/day (model), considering three different eddy detection methods applied to 1/32°R OMS model outputs in the WMS from 1993 to 2012.

Propagation of Algerian eddies
The propagation of Algerian eddies is evaluated by means of Hovmöller diagrams showing the SLA along the propagation paths previously identified in Escudier et al. (2016a). Figure 9 presents these diagrams for altimetry and model outputs over the period 2009-2015 along these paths. The comparison reveals that the number of eddies generated in this section by WMOP is quite consistent with the observations in terms of Fig. 6 Inter-annual and seasonal variability of surface currents at TAR, PAL, and GAT locations from fixed moorings (blue), altimetry (green), WMOP (red), and CMEMS Med Rea model (black) data numbers of eddies per year (approximately one large eddy per year) with significant associated inter-annual variability. However, WMOP eddies are not time-synchronous with observed eddies, in contrast to those generated by the CMEMS Med Rea model, whose timing is properly constrained by data assimilation. The most frequent propagation speed of eddies is around 5 km/day (Figure 9, yellow dashed lines) in both sections and both models. However, in the first 200 km of Section 1 and for specific years (e.g., 2012 and 2013), altimetry and CMEMS Med Rea data suggest lower propagation velocities (around 1 km/day, green dashed lines in Figure 9) that is not replicated by WMOP. This difficulty for the highresolution model to "retain" Algerian eddies in the first half of Section 1 generates significant discrepancies with observations in the eastern Alboran Sea in the presence of recentlyformed Algerian eddies, as analyzed previously for the period July-September 2012. Here again, data assimilation might help to improve the description of Algerian eddy propagation in the high-resolution model.

Discussion and concluding remarks
We have evaluated a high-resolution free-run simulation of WMOP, nested in the larger spatial resolution model of CMEMS Med Rea using various sets of multi-platform observations (HFR, satellite-borne altimeters, fixed moorings, and gliders). The effects of downscaling were assessed considering (1) the surface mean circulation both at basin and coastal Results show that WMOP improves the time-average circulation in the Balearic Sea and coastal areas of the Ebro delta and Ibiza Channel. Concretely, in the Ebro delta region, WMOP has been able to improve the location of the NC compared to its parent model, which represents a detachment too far to the east compared to observations. This improvement is confirmed using both large and coastal scale assessments. WMOP is also able to qualitatively reproduce the small-scale coastal flow intensification at the mouth of the Ebro River, contrarily to CMEMS Med Rea. However, although WMOP also slightly improves the intensity of the NC compared with the parent model, this current remains underestimated according to the HFR observations. The downscaled model was also found to improve the circulation rounding the northeast coast of Menorca Island.
The refinement achieved by WMOP in the characterization of the circulation in the Ibiza Channel is probably related to the better positioning of the NC located further north in WMOP in comparison with CMEMS Med Rea model. WMOP is able to describe the reversal of the flow produced in the channel along the zonal direction, in contrast to its parent model which generates a homogeneous northward flow over the whole channel. However, despite this improvement, the downscaled model is not able to reproduce all the finescale structures revealed by routine glider measurements. The possible errors associated with the computation of the glider geostrophic velocities from the thermal wind equation assuming the level of no motion at 800 m must be considered. In this regard, first analysis using the depth averaged velocity (DAV) approach based on the average drift of the glider have been performed at SOCIB and should help to mitigate this issue in the future. The WMOP also overestimates the intensity of the currents relative to the Ibiza Channel HFR observations. Fig. 9 Propagation of mesoscale eddies represented using Hovmöller diagrams along two sections in the Alborán Sea, shown in the inset map at the top-right corner, from a altimetry, b WMOP, and c CMEMS Med Rea. Yellow/green-dashed lines denotes "typical" propagation velocities of eddies of 5 km/day and 1 km/day, respectively. Note that for WMOP, propagation of eddies tend to be faster than in the other datasets In the southern part of the domain, WMOP provides a better time-averaged representation of the Algerian current in comparison with CMEMS Med Rea: WMOP allows the current to reach further east (< 6°E), in agreement with observations. However, it overestimates the coastal detachment of this current when reaching 6°E, in comparison with both altimetry and CMEMS Med Rea. Moreover, WMOP improves the representation of the surface flow divergence when approaching Sardinia Island around 40°N.
Regarding the validation procedures using automated oceanic eddy tracking algorithms, it should be emphasized that SLA fields from models must be re-gridded to altimeter resolution and spatially filtered (48 km) prior to comparison with altimetry in order to discard small-scale structures not reproduced by altimetry. Our results confirm that, on average, 8 eddies with a lifetime averaged radius of 40 km are produced in the WMS per day. Higher resolution observations would be needed to confirm the rates reported by the high-spatial resolution model WMOP (37 eddies with a mean size of 16 km per day).
The improvements achieved by WMOP in topographically complex areas might be linked to the use of sigma-levels compared with z-levels. The smoother representation of the topography in terrain-following coordinate models, in comparison with the step-like representation of the z-levels models, even with the use of vertical partial cells like in CMEMS Med Rea, allows better simulation of the interactions between flows and topography (Ezer and Mellor 2004) which may in turn be reflected in a better characterization of the surface flows, especially in coastal and slope regions. The use of finer temporal and spatial resolution atmospheric forcings in WMOP compared with CMEMS Med Rea may also contribute to a better representation of the general dynamics. Specifically, Schaeffer et al. (2011) demonstrated for instance how improved spatial resolution in atmospheric forcings provided more realistic upwelling events in the Gulf of Lions. Lebeaupin Brossier et al. (2011) also conducted some sensitivity tests on ocean models using different spatial resolution of atmospheric forcings in the Western Mediterranean Sea, showing that an increased resolution was able to locally strongly modify the thermohaline circulation in the Gulf of Lions.
WMOP, in contrast to its parent model which is constrained by data assimilation, is not able to correctly locate the mesoscale structures, which negatively affects the representation of the variability of surface currents. This is especially critical in the southern part of the domain, where these structures interact with the general circulation producing, in turn, significant errors in the local temporal variability of the surface currents at the locations of Cabo de Gata and Cabo de Palos moorings. The downscaling approach, through the application of initial and boundary conditions, high-resolution atmospheric forcings, and observed daily freshwater discharges, is not sufficient here to constrain the mesoscale field inside the domain and the variability of the circulation interacting with these mesoscale eddies. This limitation probably reveals the role of the chaotic intrinsic variability in this region, which produces mesoscale fluctuations and associated inter-annual and longer-scale variability (e.g., Fernández et al. 2005;Arbic et al. 2014;Penduff et al. 2018), which are not entirely determined by the external forcings of the WMOP free run simulation. In this context, data assimilation should help to improve the positioning of these mesoscale features and the variability of the currents.
In spite of this limitation, the WMOP free run simulation was shown to reproduce a realistic mean surface circulation with correct eddy statistics, therefore providing a useful highresolution inter-annual dataset to study regional ocean processes, and an appropriate basis for the future generation of reanalysis simulations that include data assimilation. Along this line, the first experiments assimilating observations in WMOP (from satellites, ARGO, CTDs, and gliders) were successfully conducted in the framework of specific oceanographic campaigns (Pascual et al. 2017;Hernandez-Lasheras and Mourre, 2018).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.