Retrieval and assimilation of velocities at the ocean surface

Ocean currents play a key role in Earth’s climate, they are of major importance for navigation and human activities at sea, and impact almost all processes that take place in the ocean. Nevertheless, their observation and forecasting are still difficult. First, direct measurements of ocean currents are difficult to obtain synoptically at global scale. Consequently, it has been necessary to use Sea Surface Height and Sea Surface Temperature measurements and refer to dynamical frameworks to derive the velocity field. Second, the assimilation of the velocity field into numerical models of ocean circulation is difficult 5 mainly due to lack of data. Recent experiments assimilating coastal-based radar data have shown that ocean currents will contribute to increase the forecast skill of surface currents, but require to be applied in multi-data assimilation approaches to allow better identification of the thermohaline structure of the ocean. In this paper we review the current knowledge on these fields and provide global and systematic view on the technologies to retrieve ocean velocities in the upper ocean and the available approaches to assimilate this information into ocean model. 10


Introduction
Surface ocean currents are recognized to critically contribute to characterize the Earth's climate (WMO, 2015).Knowledge of ocean surface velocities is a key and cross-cutting issue with impact on many societal challenges far beyond the ocean and climate research context.As such, ocean surface currents are been included in the list of essential climate variables (Bojinski et al., 2014).Indeed, ocean currents transport and redistribute heat, dissolved salts, sediments, plankton, nutrients and ocean pollutants.Strong ocean currents define corridors used by marine mammals, birds and fishes, and sustain their migrations in search for food, breeding sites and spawning areas.As a result, knowledge of the detailed structure and variability of ocean currents is required for fisheries and environmental management.On the other hand, surface currents directly affect many important socio-economic activities as global marine trade and shipping or marine pollution and safety, to mention a few.
Ocean surface currents appear as the result of a non-trivial combination of different types of periodic and aperiodic phenomena whose ranges span a continuous spectra of space and time scales, from basin-wide motions ( 1000 km) to fast narrow currents and mesoscale eddies (30-100 km wide), submesoscale features (1-10 km), and turbulence scales (1-100 m).Different components of the current ocean observing system capture different parts of this range.Land-based coastal HF radars are able to resolve rapid changes.However, although the number of HF radars has rapidly increased in the last decades, they coverage remains limited.Currents are also observed at a few moorings.Until now, Lagrangian drifters and satellite altimeter-derived surface geostrophic currents have been the only sources of information able to provide global coverage.Drifters are able to 1 Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017Discuss., doi:10.5194/npg- -14, 2017 Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.provide hourly observations but with irregular coverage with approximately one point within a 5 degree box (Dohan and Maximenko, 2010).At global scale observations using moored instruments in both the ocean surface and the water column are distributed mostly near and along the coasts, particularly in the northern hemisphere (see figure 1 Holloway, 2008;Scott et al., 2010).These are often arrays of point-based currentmeters or current profilers from ocean floor that provide limited temporal extent and concentrated inn the upper 100 m (Holloway et al., 2011).Observations in coastal regions benefit from additional efforts to keep networks for risk assessment, environmental monitoring of marine protected areas and marine security.However, spatial coverage with acoustic currentmeters configurations is still quite limited.Drifters can be deployed to reveal circulation patterns including the coastal zone however deployments there tend to be sparse in particular due to the risk of beaching and/or losing the equipment.Altimetry observations can be used to infer the geostrophic portion of surface currents on scales of several hundred kilometers and 5-10 days.Together, drifting buoys and altimetry remain the backbone of operative/operational synthesis products such as AVISO (CLS, 2016).However, as the Rossby deformation radius (providing the preferred horizontal scale of ocean structures) rapidly decreases with increasing latitude, the variability and interaction of currents with winds at mesoscale and submesoscale are not well captured as current observing systems fail to capture horizontal gradients at these scales.In short, serious gaps still exist to properly determine some key climate processes and operational applications influenced by the spatial variability of ocean currents.
As pointed out by Hogg (1996), following the original work of Levitus (1982) there have been a number of efforts to compile all past and present oceanic observations to construct a reference state of the world ocean.Most of these efforts have concentrated on the historical data of water properties as temperature, salinity, dissolved oxygen and nutrients.Large amounts of past observations have been collected, validated and calibrated from diverse sources and international archives, enhancing the geophysical consistency of the resulting climatological fields.These climatologies have shown to be useful to initialize numerical models, constrain their evolution, and validate their results.In the case of the ocean currents, the first measurements were already done by the HMS Challenger expedition (1872-76).However, during decades, the main source of ocean currents had been the ship-drift reports.Using about four million observations of ship-drift data, Richardson (1989) calculated annual and monthly mean surface currents in a 2 • × 5 • grid.Their charts served to identify large gaps in the databases and requesting the recuperation of data being withheld after the second world war.Maximenko et al. (2009) generated new global maps using the position of a global array satellite-tracked, near-surface drifters deployed since the 1980s (Niiler, 2001), altimetry and ECMWF wind reanalysis.However, due to the complexity of the current power spectra, the meaning and representativeness of the averaged and residual current measurements depend on the averaging period, its time and location (Neumann, 1968).As such, the mesoscale variability of ocean currents is still not well captured by today's global observing system.
In this context, increasing efforts are being devoted to infer complementary information of ocean currents from alternative satellite sensors.For example, ocean currents have been estimated from Sea Surface Temperature (SST) imagery using different approaches as retrieving surface temperature from the heat conservation equation (Kelly, 1989;Vigan et al., 2000b;Chen et al., 2008), neural networks (Côté and Tatnall, 1997) and the Maximum Cross Correlation technique (Bowen et al., 2002;Afanasyev et al., 2002;Dransfeld et al., 2006) based on monitoring the motion of temperature patterns in consecutive images of SST.Alternatively, a better understanding of the dynamics in the upper layers of the ocean has allowed to propose a new  Holloway (2008) framework based on the Surface Quasi-geostrophic equations (SQG Held et al., 1995;Lapeyre and Klein, 2006) able to derive surface currents from a single SST image (LaCasce and Mahadevan, 2006;Isern-Fontanet et al., 2006a;González-Haro and Isern-Fontanet, 2014).This manuscript aims to provide a review of the different approaches able to produce estimates of the surface currents from remote sensing and the gaining terrain by the synergistic combination of data from multiple sensors.Finally, the gaining experience with the assimilation of coastal current data will provide information about the gains to be expected by global data assimilation systems if real-time, quasi-synoptic maps of ocean currents were available.As a consequence, the paper is organized into two main blocks: one devoted to the approaches used to retrieve ocean velocities with remote sensing techniques (section 2) and one to the approaches used to assimilate the velocity field (section 3), which are summarized in section 4.

Retrieval from satellite observations
Surface ocean currents can be directly measured using the Doppler effect, i.e. the frequency shift of an emitted electromagnetic wave due to the relative motion between the emitter and the sea surface.This phenomenon has been exploited to retrieve them from the satellite measurements provided by Synthetic Aperture Radar (SAR, see Chapron et al., 2005).The basic idea is to separate this shift in two components: one due to the Earth rotation, and the other due to the sea surface motion.This methodology has two important advantages: it is not limited by the presence of clouds and its high spatial resolution allows measurements close to the coast.Some studies have already shown great potential in areas with very intense currents (Chapron et al., 2005;Rouault et al., 2010).There are, however, some limitations: only one component of the velocity is derived, the narrow swath limits the coverage and the retrieved current speed may contain contributions other than the ocean current.Indeed, in weak currents the dominant contribution is the wind-induced wave motion (Mouche et al., 2012).Moreover, at the moment 3 Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017Discuss., doi:10.5194/npg- -14, 2017 Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License. of writing this review, several missions able to measure the Dppler shift are under consideration by space agencies such as NASA and ESA.Some of these missions propose to use altimeters (e.g SKIM) and scatterometers (e.g.DopSCAT) to this end, other are new instruments (e.g.SeaStar).Consequently, the direct measurement of surface measurements by satellite is still quite limited at present, which makes necessary to exploit existing measurements of variables that contain information about ocean currents.
Surface currents can also be retrieved from measurements of Sea Surface Height (SSH), Sea Surface Temperature (SST) or tracers such as chlorophyll concentration invoking two conservation laws: the conservation of momentum (sections 2.1, 2.2 and 2.5) and the conservation of heat or other tracers (section 2.4).Here, we focus on the two-dimensional velocity field at the ocean surface, though this may correspond to layers at different depths depending on the approach used to derive velocities (see sections 2 and 2.2).Moreover, surface velocities have three main contributions: the velocity fields generated by density gradients (sections 2.1 and 2.5), wind (also known as Ekman current) and wave-currents interaction (section 2.2).Dynamically, currents are dominated by Earth rotation at scales large enough, implying small Rossby numbers, which is defined as with U and L being characteristic velocity and length scales and f 0 the Coriolis parameter.This motivates the expansion of dynamical variables in terms of Ro, i.e.
where v(x, z) = (u, v) is the horizontal velocity, w(x, z) is the vertical velocity, p(x, z) is the pressure and x = (x, y).At leading order, the geostrophic balance, i.e. the equilibrium between the Coriolis and pressure forces, dominates pointing to classify the different contributions into geostrophic and ageostrophic contributions.The later include both, the effect of winds and waves and O(Ro) corrections.It is worth mentioning that, although at first approach the different contributions can be computed separately, direct measurements is the result of all of them.This can make difficult to asses their relative importance and the accuracy of the approaches used.

Geostrophic currents: Sea Surface Height
At zeroth order the horizontal velocity field is non-divergent and it is possible to define a stream function ψ(x, z), that only depends parametrically on z, such that the velocity field is given by (e.g.Vallis, 2006) where e z is the unit vector in the z directions and ∇ z ≡ (∂ x , ∂ y , 0).This stream function is proportional to pressure at zeroth order, p 0 (x, z), then  with ρ 0 and f 0 been a reference density and value of the Coriolis parameter respectively.Close to the surface, the pressure field along an equipotential surface is related to the Sea Surface Height (SSH) η(x), through the hydrostatic equation.Then, surface velocity at zeroth order becomes This provides the fundamental framework that allows to retrieve surface ocean currents from the satellite measurements of 5 SSH given by altimeters (see Robinson, 2004, for more details).
Current altimeters provide measurements of SSH along the satellite track with a sampling frequency of 20 Hz, implying a spatial resolution of the order of ∼300 m.The Power Spectral Density (PSD) of these measurements show the presence of white noise (e.g.Le Traon et al., 2008;Xu andFu, 2011, 2012;Dibarboure et al., 2014;Zhou et al., 2015), which is a major limiting factor for the estimation of ocean currents.Since noise has a stronger effect on small scales it is common to low-pass 10 filter altimetric measurements before computing velocities.Nevertheless, this approach does not remove noise at large scales, which can be important in low energetic areas (e.g.Xu and Fu, 2012).Moreover, noise level strongly depends on the sea state, which makes it highly variable in space and time implying that the effective resolution of altimetric measurements is also variable.In a recent study, Dufau et al. (2016) Fontanet et al., 2016a).During the last years there have been major improvements in radar altimetry technology that not only have reduced noise levels (Dufau et al., 2016) but also have improved the capability to limit the impact of inhomogeneities in measurements (Dibarboure et al., 2014).Nevertheless, current altimeters still present strong limitations in observing small scale features O(10 km) not only due to noise but also due to temporal sampling (Chavanne and Klein, 2010).Finally, it is worth mentioning that, current altimeters still have difficulties in providing measurements at distances between 10-50 km from the coast in spite of the significative advances done during the recent years (Cipollini et al., 2017).
Altimeter measurements only allow to retrieve the velocity perpendicular to the satellite track (equation 3), as it is evident in the example shown in figure 2. Consequently, two-dimensional fields are typically obtained through the interpolation of measurements in space and time using the classical Optimal Interpolation schemes (e.g.Le Traon et al., 1998).This procedure has a two-side effect.On one side, the separation between tracks and the time sampling reduce the spatial resolution that can be achieved in comparison with along-track measurements.Indeed, Chelton et al. (2011) estimated that the shortest wavelength that can be achieved trough the interpolated two-dimensional fields is λ ∼150-200 km implaying that vortices with diameters smaller than 75 -100 km cannot be observed.This gives rise to the so-called altimetric gap, i.e. the range of scales that cannot be currently observed by altimeters.Figure 2 evidence that altimetric maps are unable to capture the signature of the smallest structures.On the other side, the limited amount of altimeters as well as the rapid evolution of some structures may induce errors in the location and geometry of ocean vortices.Pascual et al. (2006) showed that the difference between using 2 or 4 altimeters induces RMS difference in Sea Level Anomalies up to 10 cm and differences in the Eddy Kinetic Energy as big as 400 cm 2 s −2 and the comparison with drifting buoys unveils important errors in the location of some vortices (see figure 3 in Pascual et al., 2006).Moreover, Isern-Fontanet (2016) and Isern-Fontanet et al. (2017b) have shown that SSH maps derived from altimetry does not capture the fast evolving structure seen in SST of the Alboran Sea.Several efforts have been done during the last years to improve the capability to obtain two-dimensional velocities from along-track data.On one side, Ubelmann et al. (2015) have recently proposed a new approach to interpolate the sparse altimetric measurements into a regular grid based on the advection of Potential Vorticity (see section 2.5 below) during short periods of time (< 20 days) of scales smaller than the Rhines scale (∼ 300 km).This method has been recently adapted to the interpolation of along-track altimetric measurements improving the performance of the classical Optimal Interpolation schemes (Ubelmann et al., 2016).On the other side, other proposed approaches attempt to improve altimetric maps based on a two-step approach.After the standard maps are computed, the residuals to respect along track data are reinterpolated using different correlation functions that may include bathymetric constrains (see Escudier et al., 2013, and references therein).It is expected that, in the following years, the two-dimensional SSH field will be directly measured by the Surface Water and Ocean Topography (SWOT) mission using swath altimetry (Durand et al., 2010).
The geostrophic flow is O(1) and, therefore it misses the advective term in the momentum equation, i.e. v • ∇v.If the flow is considered to be axisymmetric, v(r) = v θ e θ , the advection term becomes −r −1 v 2 θ e r , where r is the radius of curvature and e r and e θ are the radial and tangential unit vectors.Momentum equations can be then easily solved giving rise to the Gradient 6 Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017Discuss., doi:10.5194/npg- -14, 2017 Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.
Wind solution (e.g.Holton, 1992).This provides a correction to the geostrophic velocities derived from altimetry that can be up to 50% of the geostrophic velocity in intense vortices (Penven et al., 2014).This corrrection, which depends on the curvature of the streamlines, can be implemented using the iterative method proposed by Endlich (1961) and Arnason et al. (1962) in the atmosphere consisting on the iteration point by point using until the improved velocity falls below a threshold or it starts to increase Penven et al. (2014).

Ageostrophic currents: wind and waves
Altimeter-derived geostrophic currents only account of a part of the surface circulation.In particular, they does not provide wind-driven currents which, in principle, can be derived from satellite measurements of wind provided by scatterometers.Since Ekman (1905), different models have been proposed to derive wind-driven currents.Many of such models are based on the momentum equations for a steady, Boussinesq flow in hydrostatic balance, i.e where v(x, z) = (u, v) is the total horizontal velocity field, v S (x, z) = (u S , v S ) the stokes drift, τ (x, z) = (τ x , τ y ) the turbulent stress, b(x, z) = −gρ/ρ 0 is buoyancy and p(x, z) and ρ(x, z) a perturbation pressure and a perturbation density to respect a reference density ρ 0 , such that |ρ| ρ 0 and which has associated a reference pressure given by ∂ z p 0 = −gρ 0 and g gravity.
As in section 2.1, the Rossby number is assumed to be small allowing to remove non-linear terms from the equation.Vertical boundary conditions are with τ w been the surface wind stress and z = −h the no-stress depth.Turbulent stress is commonly parametrized as a simple gradient transfer eddy-viscosity model with A v (z) been the eddy viscosity (e.g.Polton et al., 2005;Cronin and Kessler, 2009;Wenegrat and McPhaden, 2016).It is common to write equation 5 in its complex form as with Ṽ (x, z) = u + iv, ṼS (x, z) = u S + iv S , τ (x, z) = τ x + iτ y and ∇ = ∂ x + i∂ y . 7 Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017-14,2017 Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.
At the ocean surface, the variables that appear in equations 9 and 10 can be retrieved from satellite observations.Indeed, the perturbation pressure at the ocean surface can be derived from SSH measurements provided by altimeters p s (x) = ρ 0 gη; buoyancy can be written in terms of SST T s (x) provided by infrared and microwave radiometers and SSS S s (x) also by microwave radiometer as where ρ 0 = ρ(T 0 , S 0 ), α T is the thermal expansion coefficient and β S the saline contraction coefficient; and the wind stress τ w can be derived from scatterometer measurements.This approximation, which does not includes Stokes drift ( ṼS term in equation 9), is used to generate ocean current products such as OSCAR by NOAA (Lagerloef et al., 1999;Bonjean and Lagerloef, 2002;Johnson et al., 2007) and GEKCO by LEGOS (Sudre and Morrow, 2008;Sudre et al., 2013).The Stokes drift contribution, on the contrary, is difficult to retrieve from satellite observations.It is defined as the difference between the Eulerian velocity and Lagrangian velocity due to wave motion averaged over a wave period.And, for a monochromatic wave, it is given by (Phillips, 1977) where a w is the wave amplitude, ẽk the direction of propagation in complex notation, k w is the wavenumber and σ w the wave frequency (improved approximations to the Stokes drift can be found in Breivik et al., 2016, and references therein).Some necessary information about surface waves can be retrieved from altimeters which provide, in addition to SSH, the Mean Square Slope (MSS, related to returned power) and Significant Wave Heigh (SWH, related to returned signal shape) and from these values it is possible to empirically retrieve the mean period (Gommenginger et al., 2003).These measurements, however, miss the direction of propagation.Synthetic Aperture Radars provide information on the wave spectra but only for wavelengths longer than 150 m, typically.In practice, wave parameters are estimated from wave models (e.g.Hui and Xu, 2016).
Since the momentum balance of equations 9 and 10 is linear and, assuming that pressure gradients are not related to local wind nor waves, they are often separated into a geostrophic velocity field Ṽg , which depends on the pressure gradients and can be derived from SSH measurements (equations 1 and 3), and an ageostrophic field Ṽa driven by wind and waves.Ekman This solution only depends on the wind stress and the constant value given to A 0 , where d E = 2A v f −1 is the Ekman depth (see figure 3).On the contrary, if turbulent stress (equation 8) is assumed to be a linear function of depth, i.e. the resulting ageostrophic velocities are which is the so-called slab model characterized by a vertically homogeneous ageostrophic velocity field.Both solutions depend on τ w (x), which be retrieved from satellite measurements, and some parameters, i.e.A 0 and h, that have to be determined.
Notice, however, the key differences between these two solutions.Ekman solution has the ageostrophic velocity field that decrease with depth and surface velocities are at π 4 rad to the right (left) of wind in the Northern (Southern) Hemisphere while in the slab model solution velocities are vertically homogeneous in the upper layer and surface velocity is at π 2 rad to the right (left) of wind in the Northern (Southern) Hemisphere.The main approaches to retrieve wind-induces currents usually does 9 Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017-14,2017 Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.not attempt to reconstruct the vertical profile of velocities but focus on determining the average motion of a layer and may take into account the singularity at the equator due to the Coriolis parameter (e.g.Lagerloef et al., 1999).Notice that other parameterizations of turbulent shear, e.g. through the dependence of the eddy viscosity on wind stress or shear of turbulence fluxes, are possible (see Wenegrat et al., 2014, and references therin).
In practice, the wind-induced ageostrophic contribution of the velocity field is estimated from satellite observations using a physically-based statistical model calibrated with independent observations of the velocity field, typically surface drifters (e.g.Lagerloef et al., 1999;Rio and Hernandez, 2003;Poulain et al., 2012) .The most widely used model is of the type where B and θ are the constants that have to be fitted to the observed velocities (Ralph and Niiler, 1999;Rio and Hernandez, 2003;Poulain et al., 2009;Chiswell, 2016).As a consequence, the resulting velocities derived from satellite wind measurements will be representative of the motion at the depth of measurements.The angle θ observed using SVP drifters, which represents the motion of a 10 m layer centered at 15 m deep (see Lumpkin et al., 2017, and references therein) has been found to be within the range between 20 • and 60 • for the global ocean and the Eastern Mediterranean sea (Rio and Hernandez, 2003;Poulain et al., 2009).Moreover, Poulain et al. (2009) found very small differences in the direction between the SVP and CODE (drogued at ∼ 1m Lumpkin et al., 2017) buoys in the Mediterranean sea when fitting the model given by equation 16.On the contrary, Rio et al. (2014) found large differences between angles using SVP and Argo drifters with a geographical and seasonal dependence.
These approaches, in general, does not take into account the contribution of waves.The interaction of the Stokes drifts with planetary vorticity introduce and additional force on the momentum equations known as the Coriolis-Stokes force.As a consequence, the solution of the ageostrophic component of the velocity has additional terms given by (Polton et al., 2005) see, in general, differences between different types of drifters (Rio et al., 2014).Interestingly, this differences are very small in the Mediterranean (Poulain et al., 2009(Poulain et al., , 2012)).The parametrization of turbulent stress in terms of the velocity field points to combine equations equation 5 and 8 and solve the resulting second-order linear equation for the velocity.However, an alternative approach is obtained differentiating equation 5 and manipulate it to obtain an equation for the turbulent stress τ (x, z) known as the Generalized Ekman Model (Bonjean and Lagerloef, 2002;Cronin and Kessler, 2009;Wenegrat and McPhaden, 2016) or the Turbulent Thermal Wind Balance (Gula et al., 2014;McWilliams et al., 2015): Once shear has been retrieved, velocity can be computed using equation 9.This is the approach used by the OSCAR product without the Coriolis-Stokes term.This approach improves the solution of Lagerloef et al. (1999) and has been extensively validated (e.g.Bonjean and Lagerloef, 2002;Johnson et al., 2007).Recently, Wenegrat and McPhaden (2016) have provided and approximate general solution to this equation based on Green's function given by where and G(z, s) is the Green's function given by equation 9 in Wenegrat and McPhaden (2016).This solution is quite general and admits different parameterizations of turbulent viscosity coefficient A v (z).In addition, this solution can also include the forcing from buoyancy and the effect of wave.
Finally, it is worth mentioning that the use of forcing data (SST and SSH) with different effective resolutions in equation 5 may induce unphysical imbalance associated to the different spatial resolution of products such as SST (of the order of 1 km for IR radiometers) and SSH (of the order of 50-100 km for altimetric maps).Consequently, the spatial resolution of this approach is limited by the field with lower effective resolution.A possible approach to increase the spatial resolution of altimetric maps (see the discussion in section 2.1) consists in merging altimetric maps with Lagrangian measurements.Indeed, Taillandier  and altimetric maps, who found that not only it is possible to restore some of the variability missed in altimetric maps but also ageostrophic contributions beyond the simple Ekman model.Obviously, this approach is limited by the availability of enough drifter data.

Tracer phase: singularity analysis
One of the most daunting problems in fluid mechanics is that of turbulence: to have a complete description of turbulent flows the knowledge of an infinite number of degrees of freedom is required, what makes this problem untractable following a classic deterministic approach.Nevertheless, by introducing a statistical formulation the properties of turbulent flows can be described by means of appropriate scaling relations derived from the so-called structure functions (Frisch, 1995).The realization of this fact can be traced back up to the seminal works by Kolmogorov and all the extensive literature on turbulence developed afterwards (Novikov, 1994;Frisch, 1995).But it is not until the last decades of the 20th century that the abstract entity (multifractal hierarchy) used for describing the observed anomalous scaling in the structure functions (Parisi and Frisch, 1985) started to be interpreted as something more geometrical and directly linked to the properties of each realization of the flow (Bacry et al., 1993).The application of the so-called multifractal formalism to geophysical flows was introduced by about the same time (Davis et al., 1994), including ocean variables (Seuront et al., 1999).But it was not until the beginning of the 21th Century that multifractal formalism was not formulated in more geometrical terms, with the introduction of specific techniques 12 to compute singularity exponents from 2D maps of scalars, typically remote sensing images (Turiel et al., 2005;Isern-Fontanet et al., 2007;Turiel et al., 2008).Singularity exponents are dimensionless (no physical unit) variables with measure the local degree of regularity (if positive) or irregularity (if negative) of the scalar at each point.The set of singularity exponents do not only provide information about the statistics of changes of scale in the scalar, but also the specific geometrical arrangement of the structure explaining those changes in scale.
A striking feature of singularity exponents is that singularity isolines, especially those associated to more singular (i.e., more negative) values, seem to delineate with remarkable accuracy the streamlines of the flow, much more closely than the isolines of the scalar from which they are derived (see, for instance, figure 8 in Turiel et al. ( 2009)) -notice however that no theoretical proof of this observed property has been given so far.In Figure 4  This correspondence between singularity lines and streamlines motivated the introduction of a simple method (called Maximum Singular Stream function Method or MSSM) that provides an estimate of a unitarized stream function from the singularity exponents obtained from a map of a given ocean scalar.However, the MSSM is not very useful for dynamic studies, as it just gives information on the geometry of the flow, but neither the modulus of the velocity vector nor the sense of the circulation (upstream or downstream the depicted streamlines) are known; in particular, the MSSM cannot be used for Lagrangian studies.
Besides, by construction the MSSM relies in the capability of the so-called Most Singular Manifold (MSM) to describe the full geometry of the flow, something that introduces a certain degree of quality loss in the method due to numerical degradation.
It has not been until quite recently that the importance of knowing all the singularity exponents, both more and less singular, has been put in evidence.It is by the precise knowledge of all singularity exponents that the structural correspondence in the phase of different scalars can be used to put the underlying into correspondence and, for instance, improve the quality of a noise or damaged maps by fusing it with another map of a different scalar acquired with better quality (Umbert et al., 2014).Additionally, it has been shown that the statistical properties derived from the singularity analysis can be translated to Lagrangian studies (Hernandez-Carrasco et al., 2011), so singularity analysis can provide information about horizontal mixing and particle dispersion.The power of singularity analysis lies on its capability of extracting the phase of different ocean

Tracer conservation: sequence of images
The apparent motion of surface tracers such as SST and chlorophyll concentration suggest the use of sequences of satellite images to retrieve the velocity field that originated this motion.This is being done using two main approaches: feature tracking and inverting the conservation equation for the tracer, which, in general is given by where c(x, t) can be SST or chlorophyll concentration or even the MSS and Ċ are the sources and sinks of this tracer, including the vertical advection contribution, i.e. −w∂ z c.It is important to realize the advection term v•∇ z C is the inner product between velocity and tracer gradients, which implies that only the velocity component parallel to tracer gradient can be retrieved by inverting equation 25.This is what is known as the aperture problem.On the other side, the wealth of satellite measurements of SST points to the their use for estimating ocean currents although this is not necessarily the best choice in certain situations.
The skin depth of SST is of the order of a few µm implying that air-sea interactions can mask the presence of oceanic structures.
Moreover, the algorithm used to retrieve SST introduce additional noise.Therefore, in some situations Brightness Temperature (BT) is better suited than SST for the estimation of currents (e.g.Bowen et al., 2002;Isern-Fontanet and Hascoët, 2014).
Notice, howevere, that BT does not contains the atmospheric corrention implying that temperatures are lower and atmospheric patterns may contaminate the image.Chlorophyll concentration, on the other hand integrates information of the upper tens of meters and is able to outline ocean patterns better than SST.Nevertheless, less images are available sionce Ocean color data can only be used during daytime.In any case, the use of ocean color and SST data are limited by the need of have cloud-free sequences of images.
The standard approach used in feature tracking is the so called Maximum Cross-Correlation method (Emery et al., 1986;Bowen et al., 2002;Barton, 2002).The underlying idea is quite simple: given a template of N x × N y grid points in an image at time t 0 , it consists in searching which subwindow of size N x × N y has the maximum cross-correlation within a larger search window in an image at time t 0 + ∆t and take the displacement vector between images as the velocity field.This approach has been mainly applied to SST (e.g.Dransfeld et al., 2006;Castellanos et al., 2013;Doronzo et al., 2015) although recently it has been also applied successfully to ocean color data (e.g.Yang et al., 2015;Hu et al., 2017).An alternative approach consists on tracking the biogenic surface slicks.These slicks form monomolecular slicks that modify surface tension and therefore affect capillary waves reducing the backscatter or microwave radar emissions.This allows to observe such slicks in MSS images provided by SAR and use the MCC technique to retrieve currents.This approach was successfully tested by Qazi et al. (2014), who used SAR data from Envisat and ERS-2 separated by only 30'.Although the use of SAR data allows to overcome the limitation imposed by cloud-coverage, the interpretation of MSS is strongly dependent on weather conditions (Robinson, 2004;Kudryavtsev et al., 2005) implying that this approach can only be applied for winds within the range 2-7 m/s (Qazi et al., 2014).Marcello et al. (2008) proposes to improve the MCC approach using a two-step procedure: in the first step image segmentation is used to unveil the patterns present in the image, which are tracked in the second step.This tracking combines MCC vectors and Optical Flow methods, i.e. inversion of equation 25 with Ċ = 0.In general, the resulting velocity field is sparse and is post-processed to retrieve a smoother field (e.g.Afanasyev et al., 2002) or it is combined with altimetric measurements (e.g.Abraham, 1998;Wilkin et al., 2002).Notice that, the MCC approach requires high resolution data such as the observations provided by infrared and visible radiometers (resolutions ∼ 1 km) but the resulting velocity field has spatial resolutions of the order of the window used to track features (∼ 20 km, e.g.Bowen et al., 2002).
An alternative to feature tracking is to solve the heat equations, which provides an equation for the evolution of SST.
Integrating over the Mixed Layer (ML), the heat equation can be written as where Q(x, t) are the heat fluxes, κ is the thermal diffusion, w e is the entrainment velocity at the base of the ML which is non-zero only if there is a deepening of the ML (e.g.see Klein and Hua, 1990) and T d is the temperature below the ML.In the ocean the Péclet number is smaller than one implying that the diffusion term can be removed from equation 26.As outlined above, only the cross-isotherm component of the velocity can be retrieved unless additional constrains are taken into account.
To solve this problem, Kelly (1989) and Kelly and Strub (1992) used horizontal divergence ∇ z • v and the vertical component of vorticity (∇ × v) z as regularizing constrains for the cost function given by with Ṫ (x, t) been the source terms in equations 26 and a and b penalty parameters to tune the the influence of vorticity and divergence, which has been solved using a wide variety of numerical schemes (Kelly, 1989;Vigan et al., 2000a;Chen et al., 2008).An alternative approach to solve the aperture problem consists on using background velocity information (Piterbarg, 2009) such as altimetry (Rio et al., 2016).In that case, the velocity field is given by where v alt (x) is the velocity field given by altimeters.This methodology has the same problem than MCC, it requires a sequence of cloud-free images.Nevertheless, since it does not attempt to track features it could be applied to low resolution SST data such as the measurements provided by microwave radiometers, which are not affected by clouds.This could contribute to correct the topology of SSH fields if not enough altimeters are available (see discussion in section 2.1).

Potential vorticity inversion: synergy of sensors
The above methods rely on altimetric measurements to obtain the geostrophic component of the velocity field.As it was discussed in section 2.1, altimeters are limited by current technology (noise level, distance to coast) and sampling geometry (difficulty to retrieve two-dimensional currents), which motivates the development of approaches that exploit the characteristics of SST measurements.The necessary framework can be found at O(Ro) in the so-called Quasi-Geostrophic approximation.
Within this framework, Potential Vorticity (PV) anomaly q(x, z) is related to the geostrophic stream function (equation 2) through Alternativelly, where we assume that the bottom is far enough.Then, using the principle of invertibility of PV (Hoskins et al., 1985), the geostrophic stream function can be computed from the knowledge of surface buoyancy, that can be retrieved from SST and SSS measurements (see equation 11); N (z) that can be obtained from climatologies or density profiles from Argo buoys and the knowledge of PV.Unfortunately, PV is not known and cannot be derived from satellite measurements.Nevertheless, Lapeyre and Klein (2006) showed that the large-scale forcing in density and PV can lead to the property that the interior PV mesoscale anomalies are correlated to the surface buoyancy anomalies in the upper ocean.In that case, the PV anomaly can be separated as with ξ(z) being a function that specifies the amplitude of PV anomaly.As consequence, equation 29 can be used to retrieve the stream-function from surface buoyancy, i.e. from SST and SSS measurements.Bretherton (1966) and Lapeyre and Klein (2006) proposed to solve this problem by splitting it into two solutions: an interior solution ψ int (x, z) obtained assuming zero surface buoyancy and non-zero interior PV anomaly (b s = 0 and q = 0) and a surface solution ψ srf (x, z) obtained assuming non-zero surface buoyancy and zero interior PV (b s = 0 and q = 0).
Assuming a constant stratification N (z) = N 0 and an ocean of depth H, the surface solution is (Tulloch and Smith, 2006) ψsrf (k, z) = bs where ˆstands for the Fourier transform, k = (k x , k y ) is the wavevector, k = k its modulus and n 0 ≡ f −1 0 N 0 .If H → ∞, the surface solution becomes the classical Surface Quasi-Geostrophic (SQG) solution given by ψsrf (k, z) = bs (Held et al., 1995;Lapeyre, 2017).On the other side, the interior solution is (e.g.Klein et al., 2010), which corresponds to the baroclinic mode.The relative dominance of each solution can be separated by a critical wavelength that depends on the large scale properties of the flow (Lapeyre, 2009;Klein et al., 2010).Additional expressions can be obtained taking, for example, an exponential stratifications (e.g.LaCasce, 2012).At the ocean surface, ψ srf dominates and it is not orthogonal to ψ int (Lapeyre and Klein, 2006;LaCasce, 2012), which was used by Lapeyre and Klein (2006) to propose to approximate the total solution by the surface solution introducing an effective Brunt-Väisälä frequency n e that has to be adjusted using independent observations.Then, the three-dimensional geostrophic stream function and buoyancy can be retrieved from satellite measurements of SST as (Isern-Fontanet et al., 2008) b These equations are known as the effective SQG (eSQG) model.It is worth mentioning that, the parameter n e contains the contribution of interior PV as well as the effect of SSS, if salinity measurements are not used to derive the geostrophic velocities (see Isern-Fontanet et al., 2008).Moreover, using the relationship between SSH and the stream function (section 2.1), the above equations can be written for SSH (Isern-Fontanet et al., 2008) bs (k, z) = n e gk η exp(n 0 kz) (40) Notice that, within this framework, SST and SSH contain the same information and, once buoyancy and the stream function are known at all depths, vertical velocities can be estimated (Lapeyre and Klein, 2006;LaCasce and Mahadevan, 2006;Klein et al., 2009;Isern-Fontanet and Hascoët, 2014).LaCasce and Mahadevan (2006) and Isern-Fontanet et al. (2006b) demonstrated, for the first time, that this approach can be used to derive ocean currents from real SST measurements and, later, Isern-Fontanet et al. (2016b) that it is also possible to use SSS from SMOS.Moreover, the eSQG approach has shown to provide good results in highly variable areas such as the Alboran Sea (Isern-Fontanet, 2016;Isern-Fontanet et al., 2017b) and for small (∼ 10 km) coastal eddies (Isern-Fontanet et al., 2017a).The validity of the SQG approach has been extensively investigated using both, numerical models and real data (Lapeyre and Klein, 2006;Isern-Fontanet et al., 2006b, 2008, 2014;González-Haro and Isern-Fontanet, 2014;Qiu et al., 2016).Results show that the Mixed Layer (ML) depth is a good indicator of the periods in which the phase shift between SSH and SST is minimal, but different from zero, and, consequently, the eSQG approach can be applied (Isern-Fontanet et al., 2014).The best situations correspond to deep ML, that are typically found in winter when smaller stratification favors the deepening of the the ML (see Klein and Hua, 1990, for a discussion on the effect of ML deepening on SST).Notice that this approximation has a limited capability to reconstruct the vertical structure of the ocean (e.g.Isern-Fontanet et al., 2008;LaCasce, 2012) which has lead to propose improved models of the upper ocean dynamics (Wang et al., 2013;Ponte and Klein, 2013;Chavanne and Klein, 2016).These models, however, require to know the geostrophic stream function at the ocean surface, which is the sought field here.
The comparison between altimetric measurements of SSH and SST images (e.g.figure 2) unveils the synergy between these two measurements.In general, current SST images provide information about the location and geometry of ocean structures 18 where F (k) can be empirically estimated combining SST and SSH measurements as This idea has been analyzed in Isern-Fontanet et al. (2014) and González-Haro and Isern-Fontanet (2014) showing that the transfer function can be approximated by a Butterworth filter with γ = 1, k c a cut-off frequency and A an amplitude that has to be determined from other measurements such as altimetric data, drifters, etc (equivalently to the n e parameter in the eSQG approach).This approach is well suited to simultaneous measurements of SST and SSH such as the ones provided by Sentinel-3 satellite from ESA.
During the recent years there have been some efforts to include the ageostrophic effects into the SQG framework.On one side, Ponte et al. (2013) included wind-driven ageostrophic contributions into the SQG dynamics.As in Lagerloef et al. (1999), they integrated equation 5 (without the buoyancy and Stokes terms) over a ML of depth h, using the parameterization of the turbulent stress given by equation 8 and using the SSH to derive the pressure but, they used SQG solution (equation 41) to retrieve the vertical variation of pressure and obtained that the total velocity is given by where v 0 (x) is the geostrophic velocity at the surface.Interestingly, the effect of wind does not appears explicitly in the above equation and is contained in the ML depth.Moreover, this solution implies that at small scales smaller than wind stress, i.e a few hundreds of km, the total averaged velocities is in phase with the geostrophic velocity.On the other side, Badin (2013) also included ageostrophic effects by re-writing the SQG using the two-dimensional semi-geostrophic equations allowing to extend the this approach to scales smaller than the Rossby radius of defotmation.

Data assimilation of ocean currents
At any time, the value of any surface ocean current measurement results from the contribution of the geostrophic component  time variability), tidal motions and the wind-driven and wave-induced turbulent dynamics.At global scales, satellite altimeters provide information about the geostrophic component every five days on a 1/3 degree grid (WMO, 2015).From the approximate 1250 surface drifting buoys of the Global Drifter Program (GDP), an hourly estimate is available on a 5 • grid (Dohan and Maximenko, 2010).However, due to the complexity of the current power spectra, the meaning and representativeness of the averaged and residual current measurements depend on the averaging period, its time and location (Neumann, 1968).As such, the mesoscale variability of ocean currents is still not well captured by today's global observing system and, with the exception of the work of Santoki et al. (2013), ocean currents have not been assimilated in global or basin-wide simulations.

High-Frequency Radar observations
Alternatively, direct measurements using HF radar technology is adequate to satisfy the necessary requirements of synoptic coverage and time resolution although being restricted to surface fields.The technology is based on the analysis of the Doppler shift associated to the scatering produced by the surface wave field (Bragg scattering).Two methodologies are presently in use: the CODA Seasonde Barrick (2008) and the Wellen radar (WERA Gürgel et al., 1999), being the differences between them the configuration for retrieving both the speed and direction.HF radar systems in coastal areas have rapidly evolved during the first decade of this century and presently the global network is composed of roughly 170 sites mostly in the west and east coast of the US and with lesser extent in Europe (Rubio et al., 2017) and Australia (figure 6).Data from radar have been extensively used for oceanographic studies in coastal regions (see an exhaustive review by Paduan and Washburn (2013) and references therein).HF radar fields combined with drifter positions may result in better retrieval of the coastal surface circulation using Lagrangian data assimilation techniques (Molcard et al., 2003;Taillandier et al., 2006Taillandier et al., , 2008) )  As it has been frequently manifested that the main source of errors in high-resolution coastal simulations is the inadequate wind stress forcing, assimilation of HF radar could improve the realism of the simulations by partially correcting surface wind forcing.However, the amount of available observations (HF radar, along-track altimetry and SST maps from satellites and vertical temperature and salinity profiles from moorings, gliders and profilers) remains sparse compared with the fast, smallscale, nonlinear dynamics characteristic of coastal areas.Moreover, although HF radars can provide information of the surface currents up to 20-70 km from the coast, the actual coverage depends on radio interferences, the time of the day, solar activity, and sea state.
The first work assimilating HF radar surface data into an ocean model was done by Lewis et al. (1998) using a nudging technique to correct the model surface current towards the HF radar estimates.Since then, and driven by the continuous expansion of the network of HF radar systems, different data assimilation approaches have been used to assimilate HF radar currents into non-linear, high-resolution ocean models: nudging (Lewis et al., 1998;Wilkin et al., 2005;Gopalakrishnan and Blumberg, 2012), sequential assimilation (Breivik and Saetra, 2001;Oke et al., 2002;Paduan and Shulman, 2004;Kurapov et al., 2005a;Oke et al., 2009) and four-dimensional variational (4DVAR) assimilation schemes (Hoteit et al., 2009;Zhang et al., 2010;Yu et al., 2012).

Nudging
The first work aiming to assimilate HF radar currents into a regional model of the Monterey Bay (California, US) was published by Lewis et al. (1998).The HF radar observations, u o , were assimilated by adding a fictitious surface wind stress term that nudged the model solution u 1 (uppermost layer) towards the observed values: with ρ being the water density and C D a drag coefficient.The data being assimilated was the 30-minute averaged surface currents, available every two hours and linearly interpolated to the time step of the model.They showed that such a continuous assimilation strategy was able to modify the model currents towards the observed direction.However, significant differences remained in the velocity field after more than 170 hours of assimilation.In particular, the reconstructed velocities remained small compared with the observed ones.The authors pointed out that errors in the Doppler retrieved currents might have been the reason for it and suggested that the HF data should be processed before assimilation.For example, by removing the divergent component from the observation field.The same approach was used by Santoki et al. (2013) to assimilate 1 • × 1 • OSCAR currents (Bonjean and Lagerloef, 2002) in a basin-wide simulation of the Indian Ocean.In this work, the current measurements from three RAMA buoys were used to assess the impact of the assimilation.The authors pointed out, although it is said that OSCAR currents do not provide an accurate representation of the meridional currents at these RAMA locations, the model performs even worse.The assimilation of OSCAR velocities reduced the deficiencies of the model at these locations.
A strategy to simultaneously update the 3D velocity field was used by Wilkin et al. (2005) in thei New Jersey coast (US).
In their application, they estimated the correlation between the surface CODAR data and the measurements provided by a moored Acoustic Doppler Current Profile (ADCP) and used them to The nudging scheme of Gopalakrishnan and Blumberg ( 2012) used a four-dimensional nudging coefficient: where the nudging coefficient, µ, is a function of the distance between the observations and each model grid point.In their work, they propose an analytic form for the nudging coefficient: where ∆r H is the horizontal separation between r o and r, R H is the nudging length-scale, Z d is the depth of influence of the surface observation and T d is a damping time-scale.Each observation may accelerate and decelerate a fraction of the water column, dissemating the corresponding stresses in the four-dimensional neighborhood of each observation.In their application to assimilate HF radar data in the Raritan Bay and the coastal waters of New York and New Jersey, they implemented the limiting case R H → 0, T d → 0, µ o = (1800 s) −1 and Z d = 2 m.The impact of the assimilation was estimated using in situ observations of the ocean currents, temperature and salinity withheld from the assimilation.They found that the verticallyprojected nudging was able to improve both the hindcasting and the 24-hour forecasts of near-surface currents and temperature.

Sequential methods
Breivik and Saetra ( 2001) used what they called a "quasi-ensemble" assimilation scheme derived from the Ensemble Kalman Filter (EnKF) introduced by Evensen (1994) to assimilate HF radar observations into a 1-km, nested, regional model of the Fedje area (Norway).The basic equations of the EnKF are: In equation ( 49), x ∈ R n represents the n-dimensional model state vecor.In an ocean model, the state vecor is usually constructed from the values of sea level, and the three-dimensional fields of temperature, salinity, horizontal currents.The superscripts a and f indicate the analysis and the forecast solutions respectively.The vecor y o ∈ R p represents the set of p observations available at the analysis time.Using a different notation for the model state vecor and for the observations reflects the fact that the observing system may provide information about physical parameters not directly modeled for by the model.As such, assimilation can only be performed if there is an observation operator, H : R n → R p , that faithfully projects The gain matrix K defined at equation ( 50) is said to be optimal (in the sense that it provides the most likely estimate of the system provided the values being observed) if the system is linear and if both forecast and observation errors are Gaussian and unbiased.In such case the covariance matrices, P f ∈ R n×n and R ∈ R p×p are enough to completely define the probability density functions of the forecast and observational error.However, as discussed by Evensen (1994), this is not the case when the dynamical laws followed by the system are non-linear.Indeed, in non-linear systems, the time evolution of Gaussian errors is not longer Gaussian, and the covariance matrix no longer describes the statistical properties of the forecast errors.
For non-linear models, Evensen (1994) proposes equation ( 51) as a Monte-Carlo estimation of the forecast error from the dispersion of an ensemble of plausible estimates of the state of the system.Specifically, let us consider an ensemble of r model states, x i (t), i = 1, . . ., r, evolving according to the non-linear system dynamics and differing because of differences in the initial conditions, external forcing or model parameters.At any time, t, the ensemble mean, x(t) = (1/r) r i=1 x i (t), and the ensemble of anomalies, x i (t) = x i (t) − x(t), can be easily calculated.If we define the matrix X (t) ∈ R n×r as the matrix whose columns correspond to the members of the ensemble of anomalies, the ensemble covariance is given by equation ( 51).The parameter α is introduced to scale the weight of the ensemble versus the observations and/or to to take into account the effect of the model error.
The quasi-ensemble proposed by Breivik and Saetra (2001) consisted of replacing the ensemble of model simulations with an ensemble of model states coming from a unique model simulation: A necessary condition for ensemble (53) to have a meaningful covariance ( 51) is that the collection of states defining the ensemble are taken from a representative model simulation.The advantage of using equation ( 53) is that, once the ensemble has been constructed, the covariance remains constant, reducing the numerical cost of the assimilation algorithm ( 49)-( 51).
The resulting algorithm has been known lately as an Ensemble Optimal Interpolation (EnOI, Evensen, 2003).
In Breivik and Saetra (2001), the radar data is available every 20 minutes, and three data assimilation cycles are used to get the initial conditions for a 6 hour forecast (Figure 7).The low cost of the EnOI made possible that such 6-hour forecasts were available at the Vessel Traffic Service within 45 minutes of the acquisition of the radar measurements.However, although equation (49) allows the correction of the three-dimensional hydrographical fields of the model (temperature and salinity), Breivik and Saetra (2001) found that the model rapidly became unstable.The reason was the nested nature of the simulation.Without correcting the external, coarse simulation, large density gradients built up between the (free) external and the (constrained) internal simulations.Therefore, they had to leave out the cross-updates of temperature and salinity.As such, the information added by the assimilation was lost after 6 hours.Years later, in Zhao et al. (2013), the approach of Breivik and Saetra (2001) was compared with the usual implementation of the EnKF (Evensen, 1994), in an experiment assimilating hourly surface currents over the Qingdao coastal waters (China).In Zhao et al. (2013), the ensemble members corresponded to the difference between successive model outputs every 6 hours during one month.Their results indicated that, although EnKF provides a better fit to independent surface currents, both EnOI and EnKF improve the simulation of the coastal surface currents.
Another seminal implementation of the EnKF to assimilate a subset of observations from an array of CODAR SeaSonde HF radars deployed along the Oregon coast was described by Oke et al. (2002).In their work, they used a stationary version of the Physical-space Statistical Analysis System (PSAS) introduced by Cohn et al. (1998) and a Time-Distributed Averaging Procedure (TDAP).Observations were low-pass filtered to remove the tidal signal, and the average during a full inertial period [0, T ], i.e. approximately 17 hours, was assimilated using an EnOI algorithm to obtain an estimate of the system at time T/2 (Figure 8).The model was then initiated at time T/4 from a true solution of the model and ran until T/2.At each time step, the model solution is corrected as: where k = 1, . . ., N k refers to the time steps of the simulation.One of the advantages of the time distributed strategy is that the model always starts from a pure model output, avoiding initialization shocks.As the assimilation increment is distributed over a quarter of the inertial period, it allows the model dynamics to adjust to the data assimilation increment, better preserving the model dynamical balances.The results were validated using data from a moored acoustic Doppler profiler (ADP) The authors found that, despite the presence of an unexplained bias in the results, the data assimilation increased the magnitude of the fluctuations of the model velocity field increasing the agreement with the observations.The authors pointed out that the assimilation of HF radar data compensated for the unrepresented signal of the wind stress forcing used in their simulation.Paduan and Shulman ( 2004) assimilated low-pass filtered Monterey Bay HF radar measurements using a two-step data assimilation approach: they used an EnOI method to update the velocity field of the first layer of the model, and a second step in which the surface velocity corrections were projected downward using Ekman theory arguments of either energy conservation or momentum transfer.They illustrated the disadvantage of only correcting the surface layer as had been done in Lewis et al. (1998).The simultaneous correction of the 3D velocity field reduced the spurious velocity shear that occurs when only the surface layer of the model is corrected.Kurapov et al. (2005a, b) used an approach similar to Oke et al. (2002) to assimilate velocity profiles measured by a set of moorings in a regional simulation of the Oregon coast.As in Wilkin et al. (2005), only the velocity field was updated and the other variables were allowed to evolve as a result of the dynamical adjustment.Disregarding the ensemble covariance between currents and the hydrography fields was justified by the weak correlation that existed between these variables but also because of the sampling error of the empirical correlations estimated by the EnOI.Their results showed that their EnOI algorithm was able to improve the solution of the model and induces significant dynamical changes.
A slightly different approach was used by Barth et al. (2008) to assimilate 2-day averaged currents in a nested simulation of the West Florida Shelf.Only the radial HF radar component was seen by the data assimilation algorithm, and the background error covariance is used to statistically extrapolate the velocity perpendicular to the radial direction.In their work the background error covariance matrix was built from a set of model simulations differing in the wind-forcing.The reference wind forcing combines the NCEP NAM (North American Mesoscale Model) with in situ wind measurements.The 6-hr wind field during the year 2004 was used to calculate a set of Empirical Orthogonal Functions (EOFs).An ensemble of 100 synthetic wind fields was created by perturbing the reference wind field with a linear combination of these EOFs with Gaussian random coefficients.The analysis step corrected both currents and hydrography.Similar to the findings of Lewis et al. (1998), the authors found that the forecast skill improves if a spatial filter is used to remove spurious barotropic waves from the assimilation increment and if the wind stress is included in the state vecor.This allows the data assimilation to correct both the state of the ocean and the forcing term.In Barth et al. (2011), a similar ensemble approach is implemented with a state vecor 25  2014).In both cases, although some reduction of the error was obtained for surface currents, mixed results were obtained by respect temperature and salinity.
The expected advantage of incorporating HF radar and in situ temperature and salinity observations from glider transects into the operational system used by the Australian Bureau of Meteorology was investigated by Oke et al. (2009).They used the Bluelink Ocean Data Assimilation System (BODAS), an EnOI data assimilation system descendant from the pioneering work of Oke et al. (2002), together with synthetic HF radar and gliders, they checked the added value that these observations would have in their operational system.They found that HF data could reduce the analysis errors by 80%, with improvements reaching 200 km beyond the radar footprint.Moreover, as HF radars are able to detect spatial structures smaller than the ones resolved by the Global Ocean Observing System, they would also help reduce sea level errors.However, glider transects were found to have only a localized impact, probably due to the short spatial scales over the shelf region.It was thus suggested that, if a glider program was to be implemented, transects should be closely spaced (around 100 km) to resolve the mesoscale variability.

4DVAR
which is a weighted average of the model-data misfit and the changes to the control variables.Although not explicitly noted, the observation operator H, the observation error covariance, R and the control variance, B are a function of time.The control vecor u(t) must be defined according to each particular application.It usually contains the initial model state (currents, temperature and salinity), the fields at the open boundaries, atmospheric forcing fields (mass and momentum) or model parameters.
The goal of the 4DVAR is to find the optimal value of the control, u * , for which the cost function (55) reaches its mimimum value.For linear and perfect systems, it has been shown that the solution that minimizes equation ( 55) can be written as ( 49)-( 50).See Lorenc (1986) for a detailed discussion.In the 4DVAR assimilation, the cost function is minimized iteratively.At of heat, mass and momentum.The tidal components of the currents were removed using a least-square fit to four diurnal and four semi-diurnal tidal lines over a 1-year period.Their results showed that the observed surface currents could be fitted by adjusting the wind stress controls and that the resulting surface currents showed skill over persistence for about 20 h.However, they found that without constraining the surface winds, the resulting solution was weakly sensitive to the control o f initial and boundary conditions after about two inertial periods.Moreover, and similarly to the findings of previous works using different data assimilation methods, they concluded that surface current observations alone were not enough to constrain the three-dimensional structure of the system.
The first implementation of a multivariate assimilation of multiple data sources including HF radar currents was done by Zhang et al. (2010) in the New York Bight using the Regional Ocean Modeling System (ROMS) model (Haidvogel et al., 2008) and its adjoint model (Di Lorenzo et al., 2007).Their data assimilation method was an incremental strong-constrain 4DVAR (Powell et al., 2008) that only adjusted the initial conditions using assimilation windows of three days, overlapping the data assimilation windows, advancing the beginning of the data assimilation window by one day.In a series of sensitivity experiments they revealed that the assimilation of HF radar currents in the model increased the current prediction skill of the model by 1-2 days.However, assimilation of surface currents slightly degraded the prediction skill of subsurface temperature.These improvement of prediction skill of surface currents by the multi-data assimilation of all the available observations was also reported by Sperrevik et al. (2015).
The ability of the assimilation of ocean surface currents to correct the position of a SST front in a regional simulation was demonstrated by Yu et al. (2012).In their experiments, they assimilated daily-averaged maps of HF radar derived surface cur-  AR and chlorophyll concentration observations.Another approach that has been investigated during the last ten years is based on the Surface Quasi-Geostrophic (SQG) framework, which only needs a single image to reconstruct the velocity field at very high spatial resolutions (∼ 10 km), if the environmental conditions are appropriate.
A large effort is also being devoted to the direct measurement of ocean currents using remote sensing techniques based on the measurements of the Doppler shift.Two complementary approaches are underway: the use of satellite platforms (e.g.SAR) and the use of land-based systems such as HF coastal radars.Presently, the main constraint of these systems is their limited sampling characteristics, which restrict them to case studies.Nevertheless, they do provide insight about the expected contribution than the assimilation of ocean currents will provide to operational oceanography.Although various approaches have been successfully used to use observations of ocean currents to partially constrain non-linear simulations of various coastal areas, and even improve the geometrical location of the temperature fronts, it has been shown that multiple data sources need to be simultaneously assimilated to better constrain the hydrograpy of the system.In addition, as the main source of errors i n these simulations, advanced multivariate methodologies (ENKF or 4DVAR) need to be used to be able to retrieve wind stress information from ocean currents to further increase the prediction skill of coastal operational systems.

Figure 1 .
Figure 1.Summary of current observations from moorings .Map available at Woods Hole Institution in http://www.whoi.edu/page.do?pid=68916.Colors indicated the availability of data, see a detailed explanation of data compilation in Geophys.Discuss., doi:10.5194/npg-2017-14,2017   Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 2 .
Figure 2. Sea Surface Temperature from MODIS Aqua with Sea Surface Height from AVISO (black lines) obtained from the combination of measurements provided by different altimeters.Lines show the available measurements in the period of ± 12 hours around the time the image was taken provided by Jason-1 (red), Envisat (blue) and GFO (purple).Arrows correspond to the cross-track geostrophic velocities.
have shown that the smallest scale that can be resolved by the new generation 5 Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017-14,2017 Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License. of altimeters is in the range between 40-50 km in high energy areas but it can be as large as 90-100 km.This variability has motivated the development of adaptive approaches to fully exploit the sampling capabilities of current altimeters (Isern-

(
1905) provided a solution to the ageostrophic part of equation 9 by setting ṼS = 0, A v (z) = A 0 and modifying the bottom boundary condition (equation 7) by u → 0 and τ → 0 as z → ∞.

Figure 3 .
Figure 3. Ageostrophic velocity field for the Ekman component (green), the 'Eulerian' Stokes component (blue), the Ekman-Stokes component (red) and the resulting velocity (black).The parameters used are the same as in Polton et al. (2005).Wind and wave propagation is in the x-direction.Arrows in the lower-right plot correspond to the transport a SVP and a CODE drifter would see obtained by integrating velocities for the layers marked with gray bands.
20) assuming the same boundary conditions as in the classical Ekman solution (equation 13).Here, ṼE is the Ekman current at the ocean surface (equation 14), ṼS the stokes velocity and d S = 1/(2k), where k is the wavevector (see equation 12).Notice that the Coriolis-Stokes forcing not only change the direction of ageostrophic current but also has a contribution with a vertical extent similar to the Ekman term.Moreover, the fitting of the heuristic model given by equation 17 to wind measurements and drifter trajectories may mix the purely wind contribution with the wave contribution.As it is evident in figure 3, surface drifters such as CODE or Argo drifters are expected to have different direction in comparison with SVP drifters.Although the determination of upper wind and wave-driven currents provided by the above equation may not be accurate, observations do 10 Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017-14,2017 Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.
Although the slab model has vertically homogeneous velocities, the inclusion of the Coriolis-Stokes induces vertical variations of the velocity field since, in general d S is smaller than the Mixed Layer Depth.In a recent paper Hui and Xu (2016) have included the Stokes-Coriolis force into the model proposed by Lagerloef et al. (1999) showing an improvement of the velocity field observed by SVP drifters to respect the standard OSCAR products, particularly in the Southern Ocean.

Figure 4 .
Figure 4. Singularity exponents derived from the Brightness Temperature of the image shown in figure 2.
we show for the matter of example the map of singularity exponents derived from the SST map shown in 2. As shown in the figure, the singularity exponents provide very detailed information about the streamlines underlying SST, and provide a constant, homogeneous value along streamlines, despite the progressive change in the amplitude of the gradient of SST.Fronts and sharp transitions in general are associated to negative values and so they are shown in white colors in the figure, but also subtler transitions (i.e., smaller amplitude gradients) are associated to negative values, what allows to uncover a more detailed view of the circulation.On the other hand, positive values (represented in dark colors in the figure) are also in correspondence with other streamlines but which have less dynamic relevance.
scalars and abstracting it to a mathematical structure that is common to all of them.The introduction of the fusion techniques mentioned above open the way to provide high-quality remote sensing ocean scalars in which the amplitude of the dynamics and sense of circulation could be provided by altimetry and similar techniques, while the geometry of the circulation could be generated by fusion with the singular structure extracted from a different, better resolved variable.13Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017Discuss., doi:10.5194/npg--14, 2017     Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 5 .
Figure 5. Sea Surface Temperature from MODIS Aqua with Sea Surface Height from AVISO (black lines) obtained from the combination of measurements provided by different altimeters.Lines show the available measurements in the period of ± 12 hours around the time the image was taken provided by Jason-1 (red), Envisat (blue) and GFO (purple).Arrows correspond to the cross-track geostrophic velocities.
Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017Discuss., doi:10.5194/npg--14, 2017     Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.but it is difficult to recover velocities (see also section 2.4).On the other hand, altimeters provide information about ocean velocities but it is difficult to recover the location and geometry of ocean structures.Moreover, within the eSQG framework SSH and SST are in phase and contain the same information.These ideas motivatedIsern-Fontanet et al. (2014)  to reconstruct the surface stream function combining SST and SSH measurements through the definition of an empirical transfer function
project surface CODAR data to to the depth.The authors 21 Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017-14,2017 Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.compared two methodologies to feed their 3D maps into the dynamical model: a continuous nudging and the intermittent melding described by Dombrowsky and De Mey (1992).Their results indicate that the intermittent corrections of the 3D ocean currents better allowed the model to freely adjust and develop than the continuous nudging of the model observations toward observations.

Figure 7 .
Figure 7. Data assimilation cycle inBreivik and Saetra (2001).Surface currents are used to initialize, every hour, a 6-hour prediction.In the initialization procedure, three cycles of EnOI are used to assimilate the current data available every 20 minutes.

Figure 8 .
Figure 8. Data assimilation cycle inOke et al. (2002).The time-distributed averaging procedure approach used to initialize the problem at time T/2 uses all the observations in the period [0,T].In their application, the time T is approximately the inertial period.
Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017Discuss., doi:10.5194/npg--14, 2017     Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.that contained only the wind forcing of the model, i.e. x = (τ x , τ y ).In that case, the implicit observation operator provides the the corresponding upper ocean surface current, i.e.Hx = u 1 .The rationale behind this approach was that too frequent assimilation of observations often produces unrealistic features that, if not dissipated, degrade the model results.They opted for correcting the main source of the model error (the wind stress forcing) rather than the state of the ocean itself.Their results were validated against independent wind and SST observations.Their results indicate that improvements in the amplitude of the wind stress drove the corresponding improvement in the SST.However, in places where the SST was driven by other factors (e.g., open boundary conditions), changes in the forcing wind had no impact.The effort of using HF radar measurements to correct (separately) wind forcing and the open boundary conditions was done byMarmain et al. (

Figure 9 .
Figure 9. Data assimilation cycle in Hoteit et al. (2009).The pair of direct model run and adjoint model run is repeated iteratively until the pre-defined convergence criteria is reached.After convergence, the solution at the center of the assimilation period is used as the restart point for the next assimilation cycle: Overlap of five days.
results indicated either the presence of deficiencies in the background error covariance, B, used by the assimilation algorithm or deficiencies in the dynamical model itself (and its forcing), leading to over-correction of the model initial condition.The 27 Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2017-14,2017 Manuscript under review for journal Nonlin.Processes Geophys.Discussion started: 9 March 2017 c Author(s) 2017.CC-BY 3.0 License.
rents defined in their 6-km grid.The ocean model was nested inside the 9-km grid Navy Coastal Ocean Model of the California Current System (NCOM).Although ROMS was the ocean model used to simulate the circulation, the data assimilation used a stand-alone linear tangent model (LTM) and its exact adjoint code (ALTM).The LTM was dynamically compatible with the non-linear model and its reference ocean state is obtained by the temporal interpolation of the ROMS trajectory, sampled every 4 hours.With the data assimilation strategy shown in Figure10, they control the initial condition.After minimization of the cost function, the initial condition was used to provide a 6-day forecast with ROMS.The model output after three days was used as a first guess for the next assimilation cycle.Although the surface winds were not corrected by the assimilation, it was found that t he assimilation of the HF radar data was able to improve the geometry of the SST front.

Figure 10 .
Figure10.Data assimilation cycle inYu et al. (2012).The data assimilation is done with the help of a linear tangent model (LTM) ans its adjoint code (ALTM).The LTM is an approximation to the linearized dynamics of the ROMS model, used for both the forecast step and to define the reference solution of the LTM model.No overlap between the different assimilation cycles.