*3.1. Theoretical Background*

Under the geostrophic approximation, surface velocities *v*(*x*) are non-divergent and are determined from a stream function *ψ*(*x*)

$$
\vec{\sigma}(\vec{x}) = \vec{e}\_z \times \nabla \psi(\vec{x}),
\tag{1}
$$

which is proportional to sea surface height (*η*(*x*), SSH), i.e.

$$
\Psi(\vec{x}) = \frac{g}{f\_0} \eta(\vec{x}).\tag{2}
$$

Here, *ez* is the vertical unit vector, *x* = (*x*, *y*), *f*<sup>0</sup> is the local Coriolis parameter, and *g* is the gravity constant.

When there is a strong correlation between sea surface buoyancy (*bs*(*x*), SSB) and potential vorticity (PV), it is possible to invert the quasi-geostrophic PV equations and derive the stream-function from SSB. For a semi-infinite ocean with constant stratification *N*(*z*) = *n*<sup>0</sup> *f*0, the stream function at any depth *z* is given by

$$
\hat{\psi}(\vec{k}, z) = \frac{1}{n\_{\ell} f\_0 k} \delta\_s(\vec{k}) \exp(n\_0 kz),
\tag{3}
$$

where ˆ stands for the Fourier transform, *ne* is an effective Prandtl ratio that takes into account the stratification and the interior PV,*<sup>k</sup>* = (*kx*, *ky*) is the wavevector, and *<sup>k</sup>* <sup>=</sup> *k* [22]. This is the so-called effective surface wuasi-geostrophic (eSQG) model.

Temperature and salinity are related to buoyancy through

$$b(\vec{x}) = -\frac{\mathcal{G}}{\rho\_0} \left[ \mathfrak{a} \left( T(\vec{x}) - T\_0 \right) + \mathfrak{f} (\mathcal{S}(\vec{x}) - \mathcal{S}\_0) \right],\tag{4}$$

where *α* < 0 is the thermal expansion coefficient, −*β* > 0 is the haline contraction coefficient, and *T*<sup>0</sup> and *S*<sup>0</sup> are reference temperature and salinity, respectively. Then, assuming that temperature and salinity anomalies are related, Equation (4) can be approximated as

$$b(\vec{x}) \approx -\frac{\mathcal{g}}{\rho\_0} \left[ a'(T(\vec{x}) - T\_0) \right],\tag{5}$$

with *α* being an effective thermal expansion coefficient that takes into account the partial compensation between salinity and temperature [28]. Introducing this approximation to Equation (3), it is possible to derive surface currents from satellite observations of sea surface temperatures (*Ts*(*x*), SST) [28]:

$$
\hat{\psi}(\vec{k}, z) = -\frac{a'}{\rho\_0 n\_\ell f\_0 k} \mathcal{T}\_s(\vec{k}) \exp(n\_0 kz). \tag{6}
$$

Notice that this approach requires settin up the energy level of these velocities because of the presence of a free parameter *α n*−<sup>1</sup> *<sup>e</sup>* to take into account both interior PV and the partial compensation of salinity and temperature [28]. Interestingly, Equation (3) can also be used to derive buoyancy from SSH [28]

$$
\delta(\vec{k}, z) = n\_e \text{g} k \mathfrak{h}(\vec{k}) \exp(n\_0 kz). \tag{7}
$$

Equations (2), (6), and (7) relate key dynamical fields with satellite observations of SST and SSH.

## *3.2. Implementation*

The diagnosis of surface currents from thermal images is straightforward using the SQG approach, if spatial resolutions below 5–10 km are not necessary. Following Isern-Fontanet and Hascoët [33], the surface stream function (*z* = 0) is computed by applying Equation (6) to channel 4 brightness temperatures (BT) instead of SST due to their lower levels of noise. Before this, data gaps in BT images are filled using the approach of Isern-Fontanet et al. [28] and they are interpolated onto a doubly periodic regular grid. Then, geostrophic velocities are obtained using centered finite-differences and the free parameter *α n*−<sup>1</sup> *<sup>e</sup>* in Equation (6) is fixed imposing that the resulting velocity field has the same energy level as the velocity field derived from altimetry [28], i.e.,

$$
\langle E\_T \rangle \equiv \langle E\_\eta \rangle,\tag{8}
$$

where *ET*(*x*) is the kinetic energy of the velocity field derived from thermal images, *Eη*(*x*) is the kinetic energy of the velocity field derived from altimetry, and · stands for the spatial average. Notice, however, that the velocity field derived from thermal images has a much higher spatial resolution than the field derived from altimetry. Therefore, *ET*(*x*) is first low-pass filtered with a cut-off wavelength of 60 km in order to determine the constant *α n*−<sup>1</sup> *<sup>e</sup>* . It is worth mentioning that it is not necessary to have simultaneous observations of BT and SSH to calibrate the velocity field. Moreover, it is possible to use any other independent measurement of the velocities such as those obtained from HF radars, current-meters, or drifting buoys.

Figure 3 shows two examples of the SST representative of the period under study. These images unveil the presence of at least three different water masses: the warm waters associated to the WAG and EAG, a very warm water filling the northeast part of the Alboran Sea as well as the nearby Algerian basin, and a tongue of cold water newly entered from the Strait of Gibraltar that separates the previous ones. For this period, the EAG was smaller than the average situation and the front separating fresher waters from the saltier Mediterranean ones, usually located in the line Almeria–Oran, was displaced allowing Mediterranean waters to penetrate the Alboran basin, see [8] and references therein. The velocity field derived from thermal images correctly reproduce the currents associated to the WAG and partially the currents associated to the EAG. However, it also shows an intense current flowing from the southeast in the middle of the basin, i.e., along the limits of warm waters filling in the northeast of the image, which is not coherent with drifter trajectories (compare Figures 2 and 3). These waters are, very likely, saline enough to dominate density and make this water mass denser than the waters of the Alboran vortices, rather than lighter as it would be expected from observed temperatures alone. Several arguments support this. First, previous in situ measurements showed that resident Mediterranean waters are saltier than the waters that form the Alboran vortices with typical salinities *S* > 38 psu, e.g., [34]. Second, SSS maps derived from SMOS show the penetration of salty waters into the Alboran basing from the east (Figure 4). Third, buoyancy derived from altimetric maps using Equation (7) shows that this water mass, which is represented by the lined area in Figure 4, has less buoyancy and, consequently, it has to be more saline. Consequently, velocities have the opposite sign compared to velocities observed from drifters.

**Figure 3.** Brightness temperature anomaly derived from AVHRR's channel 4 measurements corresponding to: 11 September at 14:12 GMT (**left**); and 12 September at 14:00 GMT (**right**). Arrows show the velocity field derived from the images.

**Figure 4.** Brightness temperature anomaly derived from satellite measurements (**top**); sea surface salinity derived from SMOS (**middle**); and surface buoyancy derived from altimetric maps (**bottom**) corresponding to: 11 September at 14:12 GMT (**left**); and 12 September (**right**) at 14:00 GMT. Lined area identifies resident Mediterranean waters obtained using the method described in Section 3.3.

#### *3.3. Phase Correction*

The first approach to compute surface currents is to assume that temperature is a good proxy of density and, thus, the approach given by Equation (5) is valid. While this is a reasonable approach for the Alboran vortices and surrounding waters [29], results (Section 3.2, Figure 3) show that this is not valid for the warmer waters of the northeast. Consequently, the direct application of the eSQG model to these images produces velocities in the wrong sense in some parts of the image during this period of the year. Moreover, SSS from SMOS is not yet able to provide the necessary spatial resolution to correct full resolution infrared SST data, is not available in real time, and still slightly underestimates the dynamical range of SSS, see [26], for a more detailed discussion.

High resolution SST can still be used to derived currents if the dynamics in the area under study is dominated by the presence of water masses with distinctive salinity. Indeed, suppose there are *M* water masses with constant salinities *Sm*, where *m* = 1, . . . , *M*. Then, Equation (4) becomes

$$b(\vec{x}) = -\frac{g}{\rho\_0} \left[ \alpha (T(\vec{x}) - T\_0) + \beta \sum\_{m=1}^{M} (S\_m - S\_0) \theta\_m(\vec{x}) \right],\tag{9}$$

where *θm*(*x*) is a function such that

$$\theta\_{\mathfrak{M}}(\vec{x}) = \begin{cases} 0 & \text{if } \mathcal{S}(\vec{x}) \neq \mathcal{S}\_{\mathfrak{m}} \\ 1 & \text{if } \mathcal{S}(\vec{x}) = \mathcal{S}\_{\mathfrak{m}} \end{cases} \tag{10}$$

with the additional property that

$$\sum\_{m=1}^{N} \theta\_m(\vec{x}) = 1.\tag{11}$$

This function segments the domain into non-overlapping areas with different constant salinities. Notice that these functions are also a function of time. If, as in the Alboran Sea, only two water masses are present, the above buoyancy can be written as

$$b(\vec{x}) = -\frac{\mathcal{G}}{\rho\_0} \left[ a \left( T(\vec{x}) - T\_0 \right) + \beta \Delta S (1 - 2\theta\_{\text{med}}(\vec{x})) \right],\tag{12}$$

where the salinity of Atlantic waters is *Satl* = *S*<sup>0</sup> − Δ*S*, the salinity of the Mediterranean waters is *Smed* = *S*<sup>0</sup> + Δ*S*, and *θmed*(*x*) gives the extension of Mediterranean waters. The major practical difficulty is to determine the extension of Mediterranean waters, i.e to determine the function *θmed*(*x*). Here, we compute *θmed*(*x*) by computing the temperature anomaly and extracting the different water masses as the connected areas with the same sign of temperature anomaly. This idea is implemented by decomposing thermal images applying the *à trous* wavelets algorithm with B3-splines [33] and then removing the largest and the shortest scales. Finally, the wide, warm, and saltier mass located in the northeast can be identified as the largest connected area with positive temperature anomaly (lined area in Figure 4). Once the resident Mediterranean water mass has been identified, the effect of salinity can be corrected. To prove the concept and avoid numerical issues associated with the sharp transition of SSS, Equation (12) is further simplified to

$$b(\vec{x}) \approx -\frac{\mathcal{G}}{\rho\_0} \left[ a'(T(\vec{x}) - T\_0)(1 - 2\theta\_{med}(\vec{x})) \right].\tag{13}$$

This approach implies that gradients of temperature and salinity tend to be aligned but with a change of sign in the different parts of the area under study (see Figure 2 in Isern-Fontanet et al. [24] and discussions in Isern-Fontanet et al. [24], Isern-Fontanet et al. [26]).

The application of Equation (6) to the corrected temperature anomalies given by Equation (13) provides more realistic results (Figure 5). Indeed, Atlantic waters now flow from west to east, as seen in previous studies [35] and in agreement with the trajectories of the drifters released during the MEDESS-GIB experiment (Figure 2). It is important to take into account that the transfer function between the stream function and temperature (or buoyancy) scales as ∼*k*<sup>−</sup>1, which implies certain dominance of larger scales over smaller scales. This has as a consequence that the incorrect determination of the buoyancy sign for the Mediterranean waters can have a non-local impact. Therefore, ocean velocities around the WAG are wrongly diagnosed if this is not taken into account. Finally, its worth mentioning that the use of this crude approach (Equation (13)) instead of the approximation given by Equation (12) may affect the sign of the small scales in the Mediterranean waters. Nevertheless, we could not investigate this limitation due to the lack of data.

**Figure 5.** Corrected brightness temperature anomaly corresponding to: 11 September at 14:12 GMT (**left**); and 12 September (**right**) at 14:00 GMT. Arrows show the velocity field derived from the images.

#### **4. Comparison between Velocity Estimations**

The eSQG model, together with the phase correction approach described in Section 3.3, were used to derive velocities from all the available thermal images (Figure 6 shows six of them). The resulting fields exhibit more temporal variability of the WAG than the velocity derived from SSH (Figures 6 and 7). Indeed, thermal images shows that the WAG changes its shape significantly at temporal scales of the order of 12 h, but this variability is not observed in altimetric maps, which show an almost stationary vortex. A second major difference in the western part of the basin is the circulation east of the Strait of Gibraltar, the area between the strait and the vortex approximately limited by the meridian located at 4.8◦W. There, altimetric maps show a northwards current while thermal images generate an eastwards current, indicating entrance of Atlantic waters into the Mediterranean Sea. The comparison between velocities derived from temperature and from sea level also show differences in the time evolution of the EAG. While the first shows an evolving vortex, altimetry shows an almost stationary bipolar structure. The differences between altimetric and thermal images are also evident comparing BT and buoyancy derived from SSH (Figure 4). It is worth mentioning that, contrary to SST, buoyancy derived from SSH tends to be more representative of the patterns below the mixed layer rather than the ocean surface [28] and, consequently, fields may not be fully comparable. Nevertheless, the temporal evolution of this vortex and its characteristics points to sampling limitations of the altimetric measurements rather than the difference between density anomalies in the Mixed Layer and below it. Finally, both fields seem to capture the separation of the west–east flow from the WAG located around longitude 3.55◦W. The exception is for the AVHRR image of 10 September at 14:24, which shows currents that do not agree with other satellite images or altimetric maps. The reason for this disagreement relies on the failure of the phase correction approach for this particular image. Indeed, the separation between the EAG and the large extension of warm waters is not good enough to correct for their contribution to the flow on this side of the basin.

The velocity fields derived from altimetry and thermal images were compared to the in situ measurements provided by the drifters of the MEDESS-GIB experiment. A first comparison comes from the juxtaposition of trajectories with the instantaneous velocity fields (Figures 6 and 7), which shows a strong tendency of drifter trajectories to move tangent to the velocity fields. There are, however, some exceptions, particularly in the area east of the Strait of Gibraltar where some drifters have swirling trajectories that do not agree with any of the velocity fields derived from satellite measurements. Although the temporal length of available thermal images is not long enough to analyze the penetration of drifters in the EAG, it is long enough to see the detachment of drifters from the WAG and their path towards the east, which is in agreement with the velocities seen from satellites. Notice that, in this area, swirling trajectories are also present and in disagreement with both velocity fields.

**Figure 6.** Sequence of corrected brightness temperature with the associated geostrophic velocities. Magenta dots correspond to the position of drifters within the the period of ±0.5 days around the map date.

The velocities derived from surface drifters were compared to surface velocities derived from satellite observations by interpolating the latter onto the position and time of drifters (see Figure 8). Since no correction for the presence of inertial oscillations was applied, only those drifter segments that did not show small loops were taken from the whole set of drifters depicted in Figure 2. The global comparison between velocity components suggests a slightly better agreement between the velocities derived from SST than those derived from SSH (Figure 9). This can be seen in the global correlations between reconstructed velocities and drifter velocities, which are *ru* = 0.63 and *rv* = 0.74 for the velocities derived from SST and *ru* = 0.48 and *rv* = 0.66 for the velocities derived from SSH. Since some ageostrophic corrections modify the speed but do not modify the direction of the geostrophic current, e.g., the cyclostrophic flow, it is interesting to compare velocity directions. Here, the correlation between the velocity directions diagnosed from satellite and those derived from drifter trajectories are quite similar: *r<sup>θ</sup>* = 0.74 for SST and *r<sup>θ</sup>* = 0.71 for SSH. Notice, however, that all correlations are well below 0.9. Another difference between the two types of velocities derived from satellite observations is their kinetic energy. Although velocities derived from SST are calibrated with altimetric data, the lack of small scale signals in altimetric maps implies that higher resolution velocities derived from SST have to be low-pass filtered for calibration (Equation (8)) but the resulting field has higher kinetic energies closer to the kinetic energy of surface drifters (Figure 9).

**Figure 7.** Absolute dynamic topography and associated geostrophic velocities. Magenta dots correspond to the position of drifters within the the period of ±0.5 days around the map date.

**Figure 8.** Drifter trajectories for the period under analysis and velocities derived from the IR radiometer (red) and the Radar Altimeter (blue) interpolated onto drifter position and time. The vertical line indicates the longitude used to separate the area east of Gibraltar from the WAG. Black lines correspond to the same trajectories shown in the upper panel.

**Figure 9.** Scatter plot between the drifter velocities (ordinate) and the velocities derived from thermal images (abscissa, **left**) and altimetry (abscissa, **right**). Green corresponds to zonal velocity, orange to meridional velocity, and black to the velocity angle with respect to the east–west direction. Diamonds correspond to the area east of the Strait of Gibraltar, dots to the central part of the WAG, and crosses to the eastern area of the vortex (see Figure 8).

Figure 8 shows a relatively good agreement in the direction of currents, i.e., small angles between velocities and drifter trajectories, in the WAG for both the velocities derived from altimetry and velocities derived from thermal images. However, there is an important difference in the velocity angles observed in the velocities derived from SST with respect to drifter velocities for values around −1.5 rad (Figure 9). This discrepancy is mainly due to the velocities in the eastern edge of the WAG. Moreover, between the Strait of Gibraltar and the WAG, altimetric velocities have angles that can be up to 90◦ with respect to the trajectories of drifters while, velocities derived from thermal images show significantly smaller angles. The above observations and the patterns seen in the velocity fields (Figures 6 and 7) suggest that the agreement between satellite derived and in situ velocities may not be homogeneous and point to define three different areas: the area east of Gibraltar (the area between the Strait of Gibraltar and meridian located at 4.8◦W), the WAG (between the meridians located at 4.8◦W and 3.45◦W), and the bifurcation area (east of the meridian located at 3.45◦W), where part of the flow is trapped by the WAG and part of it flows eastwards towards the EAG (see Figure 8). These areas were used to better quantify the agreement between in situ and satellite velocities. Taylor diagrams confirm that the ability to reconstruct velocities is better for WAG than for any other area (Figure 10), although velocities derived from temperature have slightly higher correlations for the direction. In the area east of the bifurcation point, on the contrary, the performance is better for altimetry, although correlations are quite low for the orientation of the current. The opposite situation is found in the area east of Gibraltar, in which altimetry shows very poor performance. Interestingly, the meridional velocity is better reconstructed than the zonal velocity for both SSH- and SST-derived velocities.

**Figure 10.** Taylor diagram for zonal (orange) and meridional (green) velocity components and velocity direction (black) for the area east of Gibraltar (cross), WAG (circle) and stagnation point (diamond): (**top**) velocities derived from thermal images; and (**bottom**) velocities derived from altimetry. To facilitate the comparison, standard deviations were normalized by the standard deviations of drifters given in Table 1.


**Table 1.** Standard deviation of zonal velocity (*σu*), meridional velocity (*σv*), and velocity direction (*σu*) derived from the selected drifter trajectories shown in Figure 8.

Finally, the topology of the two velocity fields derived from satellite observations was unveiled by using the methodology presented by Jiménez Madrid and Mancho [36], which has been already applied on altimeter datasets over the area of the Kuroshio Current [37]. The idea was to produce a two-dimensional field covering the region, which computes at each point the arc length of the trajectory passing for the point. More precisely, by using synthetic trajectories obtained integrating forwards and backwards for a determined amount of time and getting the distance travelled, the arc lengths of the trajectory are computed for the central instant of the time series to have the maximum possible time period to integrate, i.e., 2.5 days (see Appendix A for technical details). Figure 11 depicts the arc length plots for both velocity fields. One can see both velocity fields are able to capture the two-vortex patterns typical of the Alboran Sea but with significant differences. First, altimetric field is more stable, i.e., vortex cores are clearly delimited and the flow of Atlantic water is very well defined by large arc length values. This indicates, as already stated, that the temporal evolution of the altimetric field is slow. On the contrary, thermal field shows the effect of rapidly evolving WAG producing a vaguer structure but it is still possible to see the vortex core and its contour. It is worth mentioning that the detection of EAG is challenging due to the similar temperature it has with respect to the Mediterranean waters. It is only seen thanks to the flow of fresh waters that clearly separates both structures. Notice the strong differences between them. Second, the thermal field is able to capture the entrance of Atlantic waters while the altimetric field is not. Indeed, the shadowed area in Figure 11 is due to the particles that escape from the domain; consequently, it is not possible to track their evolution for the entire time period considered. It is worth remarking that, in the plot corresponding to the altimetric field, there are big regions of very low arc length values (purple color in the figure), showing the limitation of the altimetric fields when one approaches the coast in comparison to the velocities derived from SST.

**Figure 11.** Arc length of the trajectory for the central time of the period analyzed for the velocity fields derived from: thermal images (**top**); and altimetry (**bottom**). The white area west of the Strait of Gibraltar corresponds to the point at which particles escape through the Strait of Gibraltar when integrating backwards.

#### **5. Discussion**

The circulation in the Alboran Sea is usually dominated by the presence of two large vortices that dominate the upper 100–300 m of the ocean [34,38]. These conditions are favorable to the exploitation of the eSQG approximation to retrieve surface current, as shown by LaCasce and Mahadevan [29] and confirmed by the results shown in Section 4. This opens the door to use SST to retrieve surface currents in this area, which allows resolving both the high frequency evolution of these vortices and the smaller scales structures that surround them [12]. On the contrary, Altimeters have spatial and temporal sampling limitations. Alboran vortices, however, are large enough to be properly sampled by altimeters. Consequently, both velocities derived from infrared radiometers and velocities derived from radar altimeters show similar capabilities when compared with the trajectories of drifting buoys. Nevertheless, there are some relevant exceptions: the area east of Gibraltar. This area is not properly sampled by altimeters and, therefore, velocities diagnosed in this zone are significantly different from those observed by drifting buoys. On the contrary, velocities derived from thermal images give velocities with directions relatively close to the observed by drifters (|*θη* − *θd*| = 40.2◦ while |*θ<sup>T</sup>* − *θd*| = 19.3◦). The comparison between drifters and the geostrophic velocities derived from SST or SSH, however, has an important drawback: drifter trajectories have ageostrophic contributions not taken into account in the derivation of surface currents from SST nor SSH. In particular, drifters are affected by wind [39] and inertial oscillations as well as ageostrophic contributions such as the cyclostrophic flow associated to the curvature of the Alboran vortex. Furthermore, drifter speed in the area east of Gibraltar showed some temporal behavior that suggest a relevant contribution from tides.

Although thermal images and altimeter maps have similar performances when the flow is dominated by structures large enough to be sampled by altimeters, both approaches have different latencies, with that of currents derived from SST being of the order of few hours (typically less than 4 h) (This is the time passed between satellite measurement and the delivery of surface currents derived from them. Notice that the time needed to estimate currents from SST is of the order of the time needed to compute the FFT). Of course, velocities have to be calibrated and independent measurements such as the ones provided by altimeters are necessary; however, when only the energy level has to be fixed, simultaneous data are not needed. On the contrary, it requires some time to get accurate altimetric measurements and near-real time maps can only use past data, which further limits the capability to resolve two-dimensional structures if not enough altimeters are available [18]. Thermal images have, however, two main limitations. First, they can only be used under cloud-free situations, which restricts their applicability (particularly for those methods requiring consecutive cloud free images such optical flow and MCC). It is worth mentioning that the Mediterranean Sea is a quite favorable situation, with a relatively large fraction of cloud-free images (∼40%). For moving clouds, gaps can be filled combining several images. One possible approach would be to compute vorticity from SST for several images, calibrate them, combine vorticities for a short enough period of time, and, finally, recompute the stream function inverting the vorticity.

A second limitation, as discussed above, is the need to correct for the impact of SSS. Presently, SSS measurements are not available in real time, although important advances have been done and it is already possible to detect the signature of the Atlantic waters in the Algerian basin [26] and the key patterns in the Alboran Sea with low resolutions. This implies that other approaches are needed. Here, we explored a simple method that consists of detecting Mediterranean waters from thermal images and modifing y buoyancy to take into account the effect of salinity. Although this approach is very crude, it shows that it is able to correct the role of salinity and provide results similar to those provided by altimeters. This raises a question: To what extent can the approach here proposed be used in other areas of the ocean? As shown in Section 3.3, there are two necessary conditions: the existence of water masses with salinities homogeneous enough to be considered constant and a way to identify such water masses from existing satellite observations. Multiple strategies are possible to identify water masses. In areas dominated by the presence of vortices of distinctive salinity, e.g., the Algerian basin [26], vortex identification techniques can be used [40,41]. In areas where SST univocally identifies water masses, temperature can be used to identify the extension of water masses. In other cases, image processing techniques combined with oceanographic knowledge can be applied to satellite observations of temperature or chlorophyll to segment the image. Future improvements could include the use of new SSS products (if they are available in real time), climatological SSS, or the exploitation of the buoyancy derived from altimetry to give a better estimation of density anomalies. However,

these approaches require the implementation of accurate numerical methods that are beyond the scope of this study.

#### **6. Conclusions**

This study provides further evidence that the SQG approach is able to reconstruct velocities in the Alboran Sea. Moreover, the results show that altimetric measurements do s not capture the temporal evolution of the WAG while infrared measurements do, which can be used to derive surface velocities in this area, although the role of salinity can limit this approach. Nevertheless, the strong relationship between salinity and temperature allows identifying the water masses and correct for the contribution of salinity using image processing techniques. Once corrected, the SQG approach can retrieve ocean currents slightly better to that of near-real-time currents in general but much better in areas badly sampled by altimeters such as the area east of the Strait of Gibraltar. Although this area is far from the geostrophic equilibrium, the above results show that the good sampling capabilities of infrared radiometers allows at least retrieving the direction of ocean currents in this area. These results can be extended to other areas of the ocean, provided that water masses with relatively homogeneous salinity can be identified from remote sensing measurements.

**Author Contributions:** Conceptualization, J.I.-F. and E.G.-L.; Data curation, A.O. and M.G.-S.; Investigation, J.I.-F., E.G.-L., J.A.J.-M., A.O., E.O., M.G.-S. and A.T.; Methodology, J.I.-F., J.A.J.-M. and E.O.; Software, J.I.-F.; Supervision, J.I.-F.; Writing—original draft, J.I.-F.; Writing—review & editing, J.I.-F., E.G.-L. and J.A.J.-M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by European Space Agency grant number 4000109513/13/I-LG and by Spanish Ministry of Economy and Competitiveness and the European Regional Development Fund grants numbers ESP2015-67549-C3-1-R and CTM2016-79474-R. JIF was funded by Fundación General CSIC through Programa ComFuturo.

**Acknowledgments:** This work was funded by the European Space Agency through the GlobCurrent Data User Element project (4000109513/13/I-LG) and by the Spanish Ministry of Economy and Competitiveness, through the projects PROMISES (ESP2015-67549-C3-1-R), COSMO (CTM2016-79474-R) and the ERDF (European Regional Development Fund). Financial support by Fundación General CSIC (Programa ComFuturo) is also acknowledged. In situ measurements are available through the PANGAEA (Data Publisher for Earth and Environmental Science) repository, with the following doi:10.1594/PANGAEA.853701. Altimetric data was produced by the GlobCurrent Data User Element project and is available through its web site: http://www.globcurrent.org. Infrared images are available through the Insitut de Ciències del Mar (CSIC) at https://coo.icm.csic.es/site-page/satellite-data as well as Sea Surface Salinity Maps. We would like to thank the anonymous reviewers who help to improve the present manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Arc Length Of Trajectories**

In this part, we explain the technical details to obtain the figures for the arc length of trajectories. We define the function *M*, which associates to each initial position *x*∗ in the region the arc length of the trajectory that passes through *x*∗ at the selected time *t* ∗ and then the trajectory is advected forward and backward during an integration time *τ*. More precisely, adapting the definition from [36] to the notation used here, we have the dynamical system

$$\frac{d\vec{X}}{dt} = \vec{v}.\tag{A1}$$

Let*j*(*x*∗, *t* ∗, *t*) denote a trajectory of the dynamical system in Equation (A1) which at time *t* ∗ passes through the point *x*∗. Then, for any initial condition *x*∗ in an open set included in the area of study, we define the function *M*(*x*∗)*t*∗,*<sup>τ</sup>* as

$$M(\vec{x}^\*)\_{t^\*,\tau} = \int\_{t^\*-\tau}^{t^\*+\tau} \left[ \left(\frac{dj\_u}{dt}\right)^2 + \left(\frac{dj\_{\bar{v}}}{dt}\right)^2 \right]^{\frac{1}{2}} dt\_\prime \tag{A2}$$

where *ju*(*x*∗, *t* <sup>∗</sup>, *t*) and *jv*(*x*∗, *t* ∗, *t*) are the zonal and meridional components of*j*(*x*∗, *t* ∗, *t*), respectively. *M*(*x*∗)*t*∗,*<sup>τ</sup>* is a function that associates to each initial condition *x*∗ the arc length of the trajectory that at time *t* ∗ passes through *x*∗. The trajectory is computed in the time interval [*t* <sup>∗</sup> − *τ*, *t* ∗ + *τ*].

To perform the numerical integration of Equation (A1) for advecting particles, a fifth-order Runge–Kutta was used. To prevent numerical artifacts getting the velocity values at the untabulated points within the grid cell, we utilized bicubic interpolation in space and third-order Lagrange polynomials in time. This choice of interpolating methods is recommended in [42] to get accurate results at a moderate computational cost.
