High-Resolution Ocean Currents from Sea Surface Temperature Observations: The Catalan Sea (Western Mediterranean)

: Current observations of ocean currents are mainly based on altimetric measurements of Sea Surface Heights (SSH), however the characteristics of the present-day constellation of altimeters are only capable to retrieve surface currents at scales larger than 50–70 km. By contrast, infrared and visible radiometers reach spatial resolutions thirty times higher than altimeters under cloud-free conditions. During the last years, it has been shown how the Surface Quasi-Geostrophic (SQG) approximation is able to reconstruct surface currents from measured Sea Surface Temperature (SST), but it has not been yet used to retrieve velocities at scales shorter than those provided by altimeters. In this study, the velocity ﬁeld of ocean structures with characteristic lengths between 10 and 20 km has been derived from infrared SST using the SQG approach and compared to the velocities derived from the trajectories of Lagrangian drifters. Results show that the SQG approach is able to reconstruct the direction of the velocity ﬁeld with observed RMS errors between 8 and 15 degrees and linear correlations between 0.85 and 0.99. The reconstruction of the modulus of the velocity is more problematic due to two limitations of the SQG approach: the need to calibrate the level of energy and the ageostrophic contributions. If drifter trajectories are used to calibrate velocities and the analysis is restricted to small Rossby numbers, the RMS error in the range of 10 to 16 cm/s and linear correlations can be as high as 0.97.


Introduction
The observation of currents is of key importance for understanding the ocean and managing human activities at sea. However, observations are currently limited by sampling difficulties and the lack of satellite-based direct measurements. Indeed, the cost and difficulty of deploying in situ instruments, such as current meters or drifting buoys, lead to the need of using satellite measurements. At present, some direct satellite measurement do exist based on the Doppler signal of Synthetic Aperture Radars (SAR), but they are very limited, implying that currents have to be estimated by indirect approaches (see in [1], and references therein). To remediate such a limitation, the European Space Agency has recently started a pre-feasibility study of a new mission concept that would provide ocean surface current vectors at 1 km resolution for all the coastal ocean and shelf seas: the SeaStar concept [2]. Nevertheless, beyond the technological challenges such a mission may face, the limited availability of independent measurements of surface velocity vectors at an equivalent resolution is a major difficulty.
The most common approach to indirectly retrieve surface velocities relies on alongtrack measurements of Sea Surface Height (SSH) by altimeters, whose current sampling characteristics and noise level restricts its spatial resolution to scales larger than 50 km [1]. However, a better understanding of the dynamics in the upper layers of the ocean has shown that, under the appropriate conditions [3][4][5], their dynamics can be modeled using the Surface Quasi-Geostrophic (SQG) approximation [6]. This approach has been observed to be consistent with a wide range of ocean observations [7][8][9] and it has been successfully used to reconstruct velocities from Sea Surface Temperature (SST) measurements [10][11][12] and Sea Surface Salinity (SSS) [13]. The combination of the SQG framework with infrared measurements of SST has the potential to provide estimations of surface velocities comparable to those that will be provided by the SeaStar mission [14]. Nevertheless, its capability to reconstruct the velocity field at scales below altimeter capabilities has not been yet tested.
The Mediterranean sea is well suited for testing methods for deriving high-resolution velocities from infrared observations due to the large fraction of cloud-free images (∼40% in the Western Mediterranean, see Figure 1). Moreover, the Rossby radius in the Mediterranean is small (∼20 km) implying that the geostrophic approximation is valid at quite small scales. At the same time, the resolution attained by altimeters precludes to resolve a substantially fraction of the mesoscale dynamics in the Mediterranean, while SST observations may still remain adequate. The surface circulation in the Mediterranean is characterized by the entrance of fresh waters through the Gibraltar strait that then propagate anti-clockwise along the coast, particularly in the Western Mediterranean sea, e.g., [15]. The destabilization of these waters of Atlantic origin leads to the generation of vortices with sizes in the range of 50 to 100 km approximately, which have been studied with SSH and SST, e.g., [16][17][18][19].  Here, we explore the capability of the SQG framework to reconstruct the velocity field in its lower limits exploiting both the optimal conditions for using SST measurements in the Mediterranean and the validity of the geostrophic approach at scales significantly smaller than those captured by current altimeters. In particular, we focus on the reconstruction of the velocity field associated to coherent structures smaller than 20 km that were sampled by drifting buoys.

Theoretical Framework
Assuming that upper ocean velocities are non-divergent, a stream function can be defined such that the velocity v(x) is given by where e z is the vertical unit vector, x = (x, y), ∇ z = (∂ x , ∂ y , 0), and ψ( x, z) is the stream function. Invoking the principle of invertibility of Potential Vorticity (PV) [21], the stream function of a balanced flow can be diagnosed from the knowledge of PV in the ocean interior and buoyancy on the boundaries. In particular, the Quasi-Geostrophic (QG) PV anomaly q( x, z) is given by with the hydrostatic equation providing the boundary conditions at the surface, and at the bottom, at depth H, ∂ψ Here , and n(z) is the Prandtl ratio, which is defined as the quotient between the Brunt-Väisälä frequency N(z) and the Coriolis parameter f 0 Then, the stream function ψ( x, z) can be recovered if the PV distribution q( x, z) and the surface buoyancy b s ( x) are known. In general, the PV is not known, unless in situ measurements are available, contrary to surface buoyancy, which can be estimated from satellite observations. Lapeyre and Klein [22] noted that, under the appropriate conditions, the PV is separable and in phase with surface buoyancy, implying that it can be written as which allows to solve the above differential equation, once the vertical variation of the Brunt-Väisälä frequency or the Prandtl ratio are known. The solution to Equation (1) can be split into two contributions: a surface solution ψ sr f ( x, z) obtained taking q( x, z) = 0 and b s ( x) = 0, and an interior solution ψ int ( x, z) obtained taking q( x, z) = 0 and b s ( x) = 0 [22]. For a constant stratification n(z) = n 0 , the solutions to PV inversion problem with the PV given by Equation (5) arê for the surface solution, see, e.g., in [23], and for the interior one, see, e.g., in [24]. Here, the caretˆstands for the 2D Fourier transform, k = (k x , k y ) is the wavevector and k = | k| its modulus. Then, the total solution is given bŷ Notice that the classical SQG solution [25] ψ sr f ( k, z) =b s ( k) f 0 n 0 k exp(n 0 kz) is recovered from Equation (6) in the limit H → −∞. Beyond the above, solutions for exponential stratification [26] and constant stratification with a mixed layer [27] are also available. Lapeyre and Klein [22] and LaCasce and Mahadevan [28] proposed that the classical SQG solution Equation (8) can be used to derive the total stream function from buoyancy thanks to the non-orthogonality between the interior and surface solutions. In that case, the total solution would bê where c is a constant to be fixed with independent observations. This theoretical framework can be applied to satellite observation of SST and/or SSS assuming that buoyancy at the surface is given as where T s ( x) and S s ( x) are the corresponding sea surface observed values, α is the thermal expansion coefficient, and β the haline contraction coefficent. If only SST are used, the constant c in Equation (9) not only accounts for the interior contribution [22], but also for the partial compensation between SST and SSS gradients [3,29]. The SQG approach can be generalized using the transfer function formalism, as proposed by Isern-Fontanet et al. [29] and González-Haro et al. [30] for surface currents. In that case, the surface stream function would be written asψ where F(k) can be determined theoretically using Equations (6) and (7) or Equation (9) together with Equation (10) or, alternatively, estimated empirically using Sea Surface Height (SSH) and SST observations [29,30]. As a final remark, Isern-Fontanet et al. [3] proposed to use this framework to reconstruct the subsurface dynamics from high resolution SSH, which has been widely investigated during the recent years [27,[31][32][33], and improved through the combination of SSH and SST measurements [34,35]. In those cases, high resolution SSH are needed, which will be available when the SWOT mission will be launched.

Materials and Methods
In this study, we focus on the reconstruction of surface velocities from high-resolution SST. Then, in order to explore the capability of the framework described in Section 2 to reconstruct the velocity field in its lower limits (<20 km), we have searched for small structures in the Mediterranean simultaneously sampled by surface drifters and infrared SST images. We selected those structures that could be tracked during a time period of at least 2 days and for which drifter trajectories and cloud-free satellite measurements were sampled with a maximum time difference of less than 12 h.
The database of drifting objects trajectories with vertical extensions of less than 1 m available at our institution was searched to identify looping patterns with diameters of the order of 20 km, or less ( Figure 2). Then, those trajectories were compared to infrared SST satellite measurements. We found three coherent structures that could be followed during at least two to three days with SST images and simultaneously sampled with drifting buoys. This consisted of two coastal structures-a northwards propagating elliptical vortex with semi-axes of 5 and 10 km in front of Barcelona and a southwards propagating vortex or meander with a radius of 5 km in front of Valencia-as well as an open sea vortex with a radius between 5 and 10 km (Figure 2). The scarcity of available observations fulfilling the above requirements implied that we had to use not only standard devices such CODE drifters but also other drifting objects such as the dummies deployed during the exercises done by the Spanish Search and Rescue Service (SASEMAR). In particular, we used a dummy in horizontal position with a PK-2 buoy attached to the neck (see Figure 3). The PK-2 Satellite tracked Lagrangian drifting buoy has been developed at the ICM. It is composed of a low density polyethylene cylinder with 110 mm of diameter and 430 mm height. It hosts a high-performance GPS, with an error smaller than 6 m and a sampling frequency of 30 min [36]. The dummy in a horizontal position has almost no drag, which implies that it is very sensitive to wind, consequently, we focused on those situations with no wind. The CODE drifters were deployed during different experiments such as the CONECTA cruise [37].  The trajectories of the drifting objects were compared to the available sequence of thermal images and we found that the temporal evolution of these structures was very fast and barely captured by SST images. Two types of satellite data were used: Level 1 MODIS (Aqua and Terra) images downloaded from the NASA Goddard Space flight Centre (oceancolor.gsfc.nasa.gov (accessed on 9 September 2021)) and AVHRR images from NOAA-19 and MetOp platforms available at the ICM Coastal Ocean (coo.icm.csic.es (accessed on 9 September 2021)) in real time. MODIS data were processed using SeaDAS 7.0 to obtain SST and Brightness Temperature (BT), while AVHRR data already provided both SST and BT and no additional processing was applied. BT was preferred in front of SST due to its lower levels of noise [14]. In particular, we used the channel centered at 11 µm for MODIS and the channel between 10.3 µm and 11.3 µm (channel 4) for AVHRR.

Framework Applicability
The applicability of the framework outlined in Section 2 depends on a combination of environmental [3,29] and dynamical [4,5] conditions. Among these conditions, there is the assumption that the flow is close to the geostrophic equilibrium. Moreover, as some of the ocean structures selected are located over the continental slope (see Figure 2), it is unclear whether that the classical SQG assumption may be suitable for retrieving ocean velocities (H ∼ 200 m).
The Rossby radius of deformation R provides information about the scale above which the geostrophic approximation is valid. For a continuous stratification, it is given by [38]. Taking it as the depth of the continental shelf (H ∼ 200 m) and using stratification values of n 0 = 50 and n 0 = 200 typical of winter and summer respectively, we obtain Rossby radii of R = 10 km and R = 40 km. A more accurate determination of the Rossby radius requires to solve a Sturm-Liouville eigenvalue problem. According to Chelton et al. [39], the solution to this problem can be approximated as Its worth mentioning that this approximation is known to overestimate the Rossby radius providing radii systematically higher by about 6.5% [39]. The Brunt-Väisälä frequency N(z) was estimated from the in situ temperature and practical salinity provided by the the World Ocean Atlas (2018) [40] using the TEOS-10 equation of state for sea water [41].
Then, Equation (12) was integrated at each available point (see Figure 4). Results show that the Rossby radius in this area is shorten than 10 km, which, a priory, implies that the geostrophic approximation can be used to derive the dynamics of the small structures selected for this study. The next question is to define the transfer function that will be used to reconstruct velocities. The transfer function corresponding to the surface solution for a finite depth is given by Equation (6): The limiting behavior of this transfer function is F sr f (k) ∼ k −2 for large scales (small k) and F sr f (k) ∼ k −1 for the small scales (large k). Using representative values of summer stratification and depth above the platform and in the blue sea, Figure 5 shows that the departure of the classical SQG behavior [25], i.e., the departure from F sr f (k) ∼ k −1 , is important at scales larger than the scales here investigated pointing to the use of the classical SQG solution even on the continental platform.
The evaluation of the interior contribution given by the transfer function requires to evaluate the function ξ(z). According to Lapeyre and Klein [22], it can be approximated as where q y (y, z) and b y (y, z) are the meridional averages of PV and buoyancy, respectively. In particular, the meridionally averaged PV is approximated by taking advantage of the dominance of the stretching term. These quantities were also estimated using the World Ocean Atlas (2018) [40]. In particular, density ρ( x, z) was computed from the in situ temperature and practical salinity using the TEOS-10 equation of state for sea water [41] and, then, used to estimate buoyancy using the mean density as the reference density ρ 0 . PV was obtained from the buoyancy and the Brunt-Väisälä frequency using the approximation given by Equation (16). Once these quantities were computed, the meridional averages were obtained and the function ξ(z) was estimated by least-squares fitting At the ocean surface, ξ s ≡ ξ(0) ∼ 5 × 10 −6 m −1 in the area under study ( Figure 2) although certain seasonal variability can be observed. Using this value for ξ, the relative importance between the interior and the surface solution can be explored ( Figure 5). On one side, it is also evident that the departures from the F int (k) ∼ k −2 are observable at wavelengths longer than the ones of interest. On the other side, the relative importance of the interior solution to respect the surface solution has a seasonal dependence through its functional dependence on the stratification n 0 . Despite the interior contribution at wavelengths longer than 5 km in summer, we applied the same transfer function during the whole year:

Reconstruction of Velocities
Surface velocities were estimated from BTs using Equation (19). First, they were calibrated to the same dynamical range of the corresponding SSTs and gaps corresponding to clouds were filled using the same approach as in Isern-Fontanet et al. [3], while gaps due to land were set to zero. Surface buoyancy was calculated using the Equation of SeaWater TEOS-10 and, then, interpolated to a regular local mercator grid. Finally, a Fast Fourier Transform algorithm was applied to the interpolated surface buoyancy and Equation (9) was used to estimate the surface stream function with z = 0 and n 0 = 100. The constant c was initially set to 1, although it was then modified when compared to surface drifters. Finally, the surface stream function in the real space was recovered applying the inverse Fourier transform with a high-pass Lanczos filter with a cut-off wavelength of λ max = 70 km.
The reconstruction of velocities from thermal measurements was performed applying a high-pass filter with a cut-off wavelength larger than the size the observed coherent structures. As a consequence, we assume that the contribution from the large-scale flow discarded by the high-pass filter can be approximated by a constant value vector v ls = (u ls , v ls ). Then, the reconstructed velocity v is given by where v sqg is the velocity derived using the SQG approach and c a constant value that takes into account the contribution of the interior PV at these scales as well as the partial compensation due salinity [3,29]. The calibration of both c and v ls was obtained through least-squares fitting of the above model (Equation (20)) against velocities derived from drifter trajectories v d . The goodness-of-fit was then assessed using the Pearson correlation coefficient and the Root Mean Squared Error (RMSE) of velocities and a propagation direction, i.e., The comparison between drifters was initially restricted to observations within a time period of ±24 h around the date of the SST image. Figure 6 (upper left) shows the SQG velocity field for the coastal vortex propagating along the coast and the trajectory of the drifter, released about one day before the image capture and retrieved at 24 h later. It is clear that SQG velocities tend to be tangent to the drifter trajectory. SQG velocities were calibrated with Equation (20) using only drifter velocities with speeds v d < 50 cm/s. This is motivated by the saturation in the velocity speed shown by the SQG velocities, which is evident in the scatter plot ( Figure 6). This is typical from the ageostrophic corrections due to the increase of importance of the advective term for intense and curved trajectories. Indeed, the gradient flow approximation exhibits more intense velocities than its geostrophic flow counterpart. Fitting an ellipse to the drifter trajectory provides a major radius of 9.26 km and a minor radius of 4.53 km, suggesting Rossby numbers that can be close to 1, which implies significant ageostrophic contributions. Notice that correlations for the two components are bigger than 0.9 if only velocities such that v d < 50 cm/s are considered (Table 1). This is also true for the velocity RMSE, which drops from 32 cm/s to 10 cm/s, if the most intense velocities are not considered. Table 1. Comparison between SQG-derived velocities and velocities derived from drifter trajectories: start (∆t − ) and time (∆t + ) of the drifter data used in hours to respect the SST image date; maximum wavelength (λ max ), correlations for zonal (r u ) and meridional (r v ) velocities, and correlation for velocity direction (r θ ); and RMSE of velocity (ε v ) and direction (ε θ ).   Table 1. (Right) Scatter plot between velocities derived from drifters and those obtained from the BT images. White areas correspond mainly to land areas. Figure 6 (lower left) also shows the SQG velocity field associated to the second coastal coherent structure. The comparison between the SQG velocity field and drifter trajectories also shows a good agreement. However, its temporal evolution (Figure 7) shows that the coastal vortex moved by 0.05 • southwards during a time period of 36 h approximately while the time delay between the dummy's deployment and the available SST was of 9 h. Assuming a constant propagation velocity, we estimate that the structure has propagated ∼1.4 km between this period of time (0.01245 • ). To face this problem, we moved the Lagrangian trajectory in the opposite direction of propagation of the meander taking into account its propagation velocity and the time between the images and the deployment of the dummy. If the drifter trajectory is displaced in the opposite direction to take this into account, we observe an increase on the correlations (Table 1)  The third example (Figure 8) corresponds to an open ocean situation outside the continental shelf. Its rapid evolution could be captured by four consecutive images. A CODE drifter was trapped in the core of the vortex during one day approximately and then escaped. We focused on its evolution after leaving the vortex core because there were no clouds hiding any part of the vortex. As in the previous examples, SQG velocities tend to be tangent to the trajectory of the drifters, and the scatter plot between the reconstructed velocities and the drifter velocities is reasonably good. Pearson's correlation coefficient shows values between 0.80 and 0.97 (Table 1) and velocity RMSE of the order of 10 cm/s. The comparison was extended to their directions defied as θ( x) ≡ arctan( v). It has the advantage that it is not affected by ageostropic corrections such as the centrifugal acceleration. Then, in addition to the Pearson correlation coefficient, the RMSE for the directions is given by Figure 9 shows the scatter plot between the angles of the reconstructed velocities and those from drifter trajectories for all the examples. Results show a striking correlation of 0.9. Notice that the two vortices in front of Barcelona have correlations above 0.95 and only the vortex/meander in the southernmost part of the domain has lower correlations, i.e., ∼0.8, ( Table 1). As discussed above, Lagrangian measurements started 9 h later than the SST image was taken, which may affect the comparison for that case.

Summary and Conclusions
The validation of surface current vectors at 1 km provided by future missions such as the SeaStar mission concept will require developing and exploring complementary approaches to provided currents at similar spatial resolutions, at least in enough places and/or periods of time to carry on the validation. Here, we have investigated the potential of the the SQG approach applied in the north-western Mediterranean Sea to retrieve the velocity field of coherent structures with diameters less than 20 km from high resolution SST. Our results have confirmed that the Rossby radius is small enough to allow the use of the geostrophic approximation, although for small curvature radii the gradient flow approximation should be used. Our results have also confirmed that the SQG approximation can be used, even on the continental platform, to retrieve surface currents of coherent structures, provided that they have a strong signature in SST. These results open the door to use this area for validating future missions devoted to measuring high resolution velocity fields. Data Availability Statement: Infrared images are available through the Insitut de Ciències del Mar (CSIC) at https://coo.icm.csic.es/site-page/satellite-data (accessed on 9 Septmebre 2021).