Investigating the Formation of Submesoscale Structures along Mesoscale Fronts and Estimating Kinematic Quantities Using Lagrangian Drifters

Much of the vertical transport near the surface of the ocean, which plays a critical role in the transport of dissolved nutrients and gases, is thought to be associated with ageostrophic submesoscale phenomena. Vertical velocities are challenging not only to model accurately, but also to measure because of how difficult they are to locate in the surface waters of the ocean. Using unique massive drifter releases during the Lagrangian Submesoscale Experiment (LASER) campaign in the Gulf of Mexico and the Coherent Lagrangian Pathways from the Surface Ocean to the Interior (CALYPSO) experiment in the Mediterranean Sea, we investigate the generation of submesoscale structures along two different mesoscale fronts. We use a novel method to project Lagrangian trajectories to Eulerian velocity fields, in order to calculate horizontal velocity gradients at the surface, which are used as a proxy for vertical transport. The velocity reconstruction uses a squared-exponential covariance function, which characterizes velocity correlations in horizontal space and time, and determines the scales of variation using the data itself. SST and towed CTD measurements support the findings revealed by the drifter data. Due to the production of a submesoscale instability eddy in the Gulf of Mexico, convergence magnitudes of up to ∼20 times the planetary vorticity, f , are observed, the value of which is almost 3 times larger than that found in the mesoscale dominated Western Mediterranean Sea.


Introduction
In recent decades, submesoscale surface currents have received added attention due to their importance pertaining to the energy cascade from large to small scales in the ocean, the transport and dispersion of oil spills and plastics, and the strong vertical velocities associated with these smaller scale ageostrophic motions [1][2][3][4]. The Lagrangian Submesoscale Experiment (LASER), which took place in the Northern Gulf of Mexico in 2016, and the ongoing Coherent Lagrangian Pathways from the Surface Ocean to the Interior (CALYPSO) program, being performed in the Western Mediterranean Sea are among the most recent and extensive observational campaigns aimed at sampling and understanding the three-dimensional dynamics of submesoscale flows.
The LASER experiment was conducted in the aftermath of the Deepwater Horizon event, the largest accidental oil spill in history, in order to explore the material transport pathways at the surface of the ocean, mainly because the socio-economic damage from oil spills is associated with the beaching of floating oil. The main underlying concept of LASER was to use a very large number of drifters (in excess of 1000) to capture the spatio-temporal variability of submesoscales in the northern Gulf of Mexico. During the experiment, shortly after the deployment of a large array of drifters, a convergent submesoscale instability caused the area covered by the drifters to contract by one hundred thousand times over just a few days [5]. The strong convergence exemplified the capacity of submesoscale features to facilitate vertical transport into the interior ocean. The follow up CALYPSO program is aimed at exploring the full 3D structure of ocean eddying, and in particular identifying the role of submesoscale processes with subduction. The CALYPSO program is being conducted in the Western Mediterranean Sea, specifically centered on the Alboran Gyre systems, which are characterized by large persistent mesoscale eddies in the confluence region of Atlantic and Mediterranean water masses. Here, we compare two mesoscale frontal systems, one observed during the LASER experiment in 2016 and one observed through the CALYPSO effort in 2019, to study the causality of submesoscale formation along mesoscale fronts, as well as the dynamics within the observed submesoscale surface flows. These novel measurements of submesoscale kinematics in the proximity of larger mesoscale fronts can serve as a much needed baseline comparison to help improve operational ocean model performance.
Submesoscales in the ocean are characterized as flows with Rossby numbers (Ro = ζ/ f ) of O(1) or greater, meaning that the relative vorticity, ζ, is equal to or greater in magnitude than the planetary vorticity, f . A consequence of such a flow regime is the breakdown of geostrophic and hydrostatic assumptions that can be applied to larger scale flow, but fail to describe the 3 dimensional flows typical at smaller scales [6,7]. The 3 dimensionality of these small scale flows allows for significant vertical velocities at convergent or divergent zones, having magnitudes of divergence on the order of δ ≥ f [8,9]. These convergent or divergent associated vertical velocities act as major pathways in the exchange between the atmosphere and the interior ocean, which has large implications for the transport of nutrients and carbon storage in the ocean [10,11]. Convergence zones, specifically, collect large amounts of buoyant material, like plastic or oil, as well as plant and animal life which can have large implications for ecosystem dynamics [5]. Submesoscale flows have also been shown to have significant control over the horizontal transport and mixing at the ocean surface, thus playing a critical role in the dispersion of floating debris, spilled oil, and Lagrangian drifters [3,4]. Modeling studies have shown the importance of submesoscales as a key player in the forward cascade of energy from large scale currents towards microscale turbulent energy dissipation [12,13]. These submesoscale instabilities form along regions of high strain in mesoscale frontal features, extracting energy from the mesoscale field in the process [13]. Using a neutrally buoyant float, in Ref. [14] authors measured strong vertical motions and enhanced turbulence at a sharp submesoscale front, highlighting the role submesoscales play in ocean energy dissipation.
Among the first observations of submesoscale features were sun glitter images captured during the space shuttle Challenger 41-G mission in 1984, which showed an abundance of submesoscale eddies within fields of strong horizontal velocity shears [15]. Although images, like the ones taken aboard the Challenger space shuttle, lead to promising theories to explain submesoscale dynamics [16], they lacked the quantitative, in-situ measurements necessary to fully describe the structure and dynamics of these small scale features. One of the first studies to estimate values of surface velocity gradients using in-situ measurements, Ref. [17] used clusters of drifters and ship drift to infer surface velocities and calculated values of horizontal shear greater than several times f and cross-frontal convergence larger than 0.8 f . In Ref. [18], authors similarly utilized two vessels to take parallel ADCP measurements about 1 km apart and found cyclonic vorticity on the order of 3 f amidst a weak, anticyclonically dominated background velocity field.
Ref. [19] was able to estimate large values of horizontal velocity shear (∼10 f ) from an aerial image which captured the refraction of swell waves across a sharp submesoscale front. More recent remote sensing techniques have also showed promise in capturing detailed observations of strong submesoscale fronts. During the LASER campaign, Ref. [20] used a ship based marine X-band radar to observe sea surface roughness and calculate surface currents based on surface waves captured in the radar backscatter. Similarly, during the LASER campaign, Ref. [21] used sun-glitter images taken from an airplane to observe fronts in sea surface roughness measurements. The studies by [20,21] show very good agreement, measuring similar values of across-front convergence and velocity shear of up to 5 f at a horizontal resolution of ∼500 m. However, the high resolution sea surface roughness derived from images of sun-glitter, revealed a very narrow front measuring 30-50 m, which indicated the velocity gradients present to be from ∼45-80 f . These measurements were based on the amplitude of surface roughness contrast across the narrow front [21].
Other recent studies have developed several other methods to quantify surface gradients from Lagrangian drifter measurements. Ref. [22] investigated the scale dependency of kinematics from a large drifter data set (∼180) using drifter triad metrics, revealing that horizontal divergence and strain rates are inherently larger at smaller scales. Ref. [5] calculated values of relative vorticity and convergence on the order of 5 f within a heavily sampled submesoscale cyclonic eddy observed during the LASER campaign. These calculations were made by fitting ellipses to drifter clusters based on their center of mass and fitting a plane to the velocities of the drifters. Ref. [8] also used a drifter clustering technique to compute horizontal velocity gradients based on a linear least squares approach and found magnitudes of divergence and relative vorticity on the order of 1-10 f .
Using drifter clustering techniques to compute horizontal velocity gradients becomes challenging when drifter clusters become less isotropically spread, and begin to collapse onto a line, as drifters often do at convergent frontal zones. Ref. [8] explains how the aspect ratio of drifter clusters is the most significant source of error when calculating horizontal velocity gradients in this manner, and how large separations between drifters make the calculated gradients difficult to interpret, as they may miss energetic small scale features present between drifter locations. As a result that surface drifters tend to elongate and become one-dimensional in the presence of a front, methods reliant on drifter clustering or triads may not be well suited to describe such flow regimes [5,8,22].
Ref. [23] utilized a dense deployment of 326 drifters during LASER to develop a velocity reconstruction method, based in Gaussian Process Regression (GPR), to transform Lagrangian velocities from drifters into Eulerian velocity fields, on which gradients can easily be calculated. The GPR method is unique in that it uses a chosen covariance function, which characterizes velocity correlations in horizontal space and time, and determines the scales of variation using the data itself through an optimization process. Applying this method over a 12 h period to 326 drifters caught inside a frontal cyclonic eddy, revealed maximum values of horizontal convergence and cyclonic vorticity of 8 f and 13 f , respectively.
In executing this study we further develop, apply, and validate different GPR frameworks, the core method of which was first developed for use with Lagrangian velocities by [23]. GPR is used extensively in supervised machine-learning methods, and has a large range of applications in geostatistics and uncertainty propagation in numerical modeling [23][24][25][26][27][28]. Here, we expand upon the work done thus far by [23], which focused solely on drifters within a cyclonic convergent eddy, to include drifters which align on mesoscale fronts, as well as their subsequent behavior. To do so, we explore several new technical variations, around the core GPR method, in an effort to improve the performance of the velocity reconstruction under this new dynamical regime. Taking advantage of the coinciding drifter data and marine X-band radar measurements made by [20] during LASER, we assess the performance of each GPR framework by comparing the resulting velocity fields to simultaneously derived X-band radar velocity fields. We then utilize the best performing frameworks to carry out the quantitative comparison between the mesoscales fronts sampled in the Gulf of Mexico and in the Mediterranean Sea.
The main objective and novelty of this work is the identification, quantification, and comparison of newly observed instances of submesoscale motions along mesoscale fronts captured during the LASER and CALYPSO experiments. Most in-situ measurements of submesoscale velocity field gradients have been made near coastal outflows, where mesoscale structures are largely absent [5,21,23]. Thus, far less is known about how mesoscales and submesoscales interact. The drifter data in the Gulf of Mexico presented in this study, captures the quantitative surface signature during the formation of a submesoscale swirling eddy, created through the extraction of energy from a mesoscale front. Sea Surface Temperature (SST) data reveals the presence of multiple swirling eddies, very similar in size, forming along the same front. In the Mediterranean Sea, while we observe some submesoscale motions within the mesoscale front sampled, there is no clear sign of submesoscale formation as result of flow instability or energy extraction from the mesoscale. We compare the reconstructed surface velocities to CTD measurements in order to view the 3-dimensional structure of the flow features sampled in both regions. We then investigate the mechanisms responsible for such submesoscale formation observed along the mesoscale front sampled in the Gulf of Mexico, and the lack of such submesoscale dynamics along the mesoscale front in the Western Mediterranean Sea.
The paper is constructed as follows: Section 2.1 describes the drifter, X-band radar, SST, and CTD data sets used for the analysis. Section 2.2 describes the core GPR method, the variations to the GPR framework employed here, and the method of validation using X-band radar derived velocities. Results are presented in Section 3, followed by a discussion in Section 4. Concluding remarks can be found in Section 5.

CARTHE Drifter
Both field campaigns, in the Northern Gulf of Mexico and in the Western Mediterranean Sea, featured large deployments of biodegradable CARTHE (Consortium for Advanced Research on Transport of Hydrocarbon in the Environment) drifters at adequate spatial resolution to capture the dynamics of submesoscale features. The GPS position of each drifter was initially recorded on 5-min intervals, and were then quality controlled for missing transmissions and linearly interpolated to regular 15-min intervals. The drifters are composed of a torus float containing a Spot Trace GPS unit from Global Star and two interlocking panels that form the drogue. The drogue occupies the vertical layer from 0.2-0.6 m, with the drifter having an overall draft depth of 0.6 m. At the time of the The Lagrangian Submesoscale Experiment (LASER) campaign, a flexible rubber tube was used to connect the drogue to the circular float, however due to drogue loss during the experiment, the CARTHE drifter was updated to utilize an equivalent length of steel chain to attach float and drogue before the Coherent Lagrangian Pathways from the Surface Ocean to the Interior (CALYPSO) experiment.
Drogue loss during LASER was determined using differential velocities of neighboring drifters and the sporadic transmissions received by undrogued drifters flipped by waves. The drogue detection algorithm used by [29] was validated on a subset of 50 drifters and was able to differentiate drogued from undrogued drifters with an accuracy of 94-100%. The precision of the time of drogue loss was 0.5-3 h for 85% of the drifters [29]. Only known drogued drifters were used for the analysis in either region. The trajectories used in the Gulf of Mexico can be seen in Figures 1 and 2, while the trajectories used in the Mediterranean Sea can be seen in Figure 3. An animation of the drifters captured in the Gulf of Mexico submesoscale eddy is also shown in S1 of the Supplementary Materials. The drift properties of the CARTHE drifter were extensively tested in laboratory settings before the LASER campaign [30]. The drifters were found to follow the mean Eulerian current averaged over their draft depth within 0.01 m s −1 in the absence of wind and waves, and have a maximum total slip velocity of 0.3 cm s −1 , which decreased with increasing wind speed in the lab. Drogued drifters were also found to feel a reduced wave-induced acceleration due to the flexible tether holding together the float and drogue [30].   Figure 12 are also plotted as black or blue crosses. (b) Shows the same drifters positions as plotted in a (a), along with the trajectories and end points, ending in the southwest corner of the domain, of the drifters captured in the submesoscale eddy. Ship track of the MVP CTD transects shown in Figure 16 are also plotted. Velocity data from green and magenta drifters are both used during the velocity reconstruction over the mesoscale front, however only magenta drifters are used for the reconstrucions corresponding to the smaller scale eddy.

Marine X-Band Radar Velocities
In the Gulf of Mexico, this study utilizes corresponding Marine X-band Radar velocity measurements made during the The Lagrangian Submesoscale Experiment (LASER) campaign to validate the different GPR frameworks used here. The 9.4 GHz Doppler marine radar, developed at Helmholtz-Zentrum Geesthacht (HZG), Germany, was mounted on a mast on top of the wheelhouse of the R/V F. G. Walton Smith at about 12.5 m above mean sea level. Ref. [20] calculated, validated, and compiled all the velocity measurements made using the marine radar during LASER, a small portion of which serve as a velocity comparison for the Gaussian Process Regression (GPR) velocity fields calculated here. The marine radar has a 7.5 ft horizontal transmit and receive (HH)-polarized antenna, a rotation period of 2 s, and a 12 kW peak power output. The marine radar has a maximum range of about 3.1 km, with a range resolution of 7.5 m.
The surface current velocities were calculated using only measurements of radar backscatter intensity. Based on the wavenumber range of the signal used to calculate these velocities (0.1-0.3 rad m −1 ) the approximate effective depth of these measurements range from 1-5 m, sampling larger depths in the presence of longer waves. The X-band velocity fields are processed over a time window of ∼30 min, with individual measurements of velocity corresponding to an area of ∼0.7 km 2 and ∼4 min of backscatter intensity images. The final velocity fields have a horizontal resolution of ∼475 m with time increasing from one side of the image to the other, as the research vessel travels along its heading [20]. The complete details of the marine radar velocity calculations can be found in [20].
Alongside the X-band radar measurements, wind data were also collected via an anemometer also mounted aboard the Walton Smith, relevant data of which are plotted in Figure 4. No radar velocity measurements or anemometer data are used in the analysis of the Mediterranean Sea data.

SST Data
In the Gulf of Mexico, Sea Surface Temperature (SST) products produced by Collecte Localisation Satellites (CLS) from AVISO satellite data, provide additional observations of the water masses which shape the mesoscale front and the subsequent submesoscale features sampled by the drifters. Each SST product is a 24 h average of measurements made using four inter-calibrated satellite infrared radiometers. Each SST map has a resolution of 0.02 • (∼2 km) and each 24 h averaging period is centered at 00:00 UTC of each day. Figure 1 shows all the drifter data from the Gulf of Mexico along with two SST images which correspond roughly in time with the drifter data shown.
In the Mediterranean Sea, SST data from the NASA satellite, Aqua, measured with the Moderate Resolution Imaging Spectroradiometer (MODIS), is used to complement the drifter data in the region. SST images are derived from the Aqua/MODIS swath data for the given day, which has 1 km resolution, across and along swath. The SST image used for the analysis in the Mediterranean can be seen in Figure 3.

CTD Data
During the The Lagrangian Submesoscale Experiment (LASER) campaign in the Gulf of Mexico, towed Moving Vessel Profiler (MVP) CTD casts were made to sample through the depth of the mixed layer to observe the 3-dimensional structure of submesoscale features. A full description of this instrumentation can be found in [31]. Temperature, conductivity, and pressure were measured and gridded onto 1 m vertical bins, along with the derived variables, salinity and potential density anomaly (σ θ ). Salinity and potential density anomaly were calculated using the CSIRO (EOS-80) Seawater toolbox. Ship tracks of the transects of CTD data used in the Gulf of Mexico can be seen in Figure 2.
Similarly in the Mediterranean sea, CTD transects were also collected across the flow feature of interest. Salinity and temperature were bin-averaged to 1 dbar vertical resolution with pressure being defined at the bin centers. Density was calculated using the Gibbs-Seawater toolbox and reported here as σ (ρ − 1000). Ship tracks of the CTD transects used for the analysis in the Mediterranean are shown in Figure 3.

Velocity Reconstruction with Gaussian Process Regression
The core method of reconstructing surface current velocities measured with CARTHE drifters using Gaussian Process Regression is the same used by [23], with the novelty of this study being the varying treatments of the data before and after the reconstruction is performed. These details will be described in Section 2.2.3. Following [23], we assume the drogued CARTHE drifters being used for this study, accurately follow the water velocity over their draft depth of 60 cm. In order to perform the reconstruction, we treat the zonal u and meridional v components of the observed velocities, as separate, scalar quantities, with each individual velocity datum being associated with a unique coordinate in horizontal space and time, such that (u i , v i ) = [u (p i ) , v (p i )] , with p = p (x, y, t) , i = 1, ..., n, and n being the total number of observations. Since the method of reconstruction for u and v is identical, the procedure will only be presented for the u component.
A Gaussian process consists of a collection of random values, such that any finite number of the collection has a joint Gaussian distribution. A Gaussian process is entirely defined by its mean function and covariance function [25]. In practice, the mean function is also referred to as the prior mean, because it is used as a best guess of what a given quantity could be in the absence of data. Assuming the quantities of surface velocity we're interested in follow a Gaussian process, G P, we can define the velocity at all points in a continuum as: where u is the scalar component of the velocity in the zonal x direction, u is the mean function, and K (p, p ) is the covariance function, which specifies how velocities at any two points p and p , jointly deviate about the mean. For practical use with a finite number of points, we must write Equation (1) as a joint, or multivariate, Gaussian distribution for the velocity, defining the prior mean as a vector of velocities and the covariance function as a matrix of covariances. In doing so, we specify between our drifter measurements of velocity, denoted here as the vector u d , and the target velocities we wish to estimate, denoted as vector u t . We will also now define points p to refer to the coordinates of the observed velocities, and points p' to refer to the coordinates of the target velocities. In order to account for some degree of measurement error, we also assume that u d are noisy observations of u, the noise of which is independent, normally distributed with zero mean and variance σ 2 N . With u d and u t being defined as two finite collections of u values, each of which following a Gaussian distribution, the multivariate Gaussian distribution for u becomes: where u d and u t are the prior mean estimates of u at observed and target points, respectively, N denotes a normal or Gaussian distribution, K dd is the covariance matrix between pairs of observation points, K dt is the covariance matrix between observation and target points, (K dt ) T is the Transpose of K dt , K tt is the covariance matrix between pairs of target points, and I is the n x n identity matrix. Expanding the prior mean vectors and covariance matrix entries we obtain: Then, by conditioning the Gaussian prior distribution (Equation (2)) on the observations, u d , we arrive at the joint Gaussian posterior distribution: whereũ t refers to the predictive, or Gaussian posterior mean, and cov(u t ) refers to the Gaussian process posterior covariance. Following [25], Equation (4) serves as the predictive equations for GPR, as values of u t can be sampled from the joint posterior distribution at target points, p , by evaluating the prior means u d and u t and the covariance matrices K dd , K dt , (K dt ) T and K tt .
For the core GPR method used here, both u d and u t are effectively set to zero for all variations of the method, some caveats of which will be described in Section 2.2.3. This allows all the reconstructed velocity fields to be contingent only on the covariance function and the regression itself. To define our covariance function, K , we use the sum of 2 squared exponential functions, each having one defined scale of motion in the zonal, meridional, and temporal dimension: where r t ,r y , and r x are the correlation time scale, meridional length scale, and zonal length scale, respectively, σ 2 is the signal variance, and M is the number of scales per dimension. Using a sum of two, three-dimensional covariance functions allows the regression to capture two distinct correlation scales per dimension. As a result that the flow sampled by the drifters in the Gulf of Mexico seems to be driven by a combination of a larger, mesoscale front, and a subsequent submesoscale spiraling eddy, the use of two scales per dimension seems appropriate to capture the two dominant scales of motion present. The parameters that define the correlation scales in the squared exponential functions r ti , r yi , and r xi , along with the velocity variance, σ 2 i , and the Gaussian noise variance σ 2 N , are known as the hyperparameters of the covariance function. A distinct feature of GPR is the ability to optimize all the correlation scales, and the signal to noise ratio, using the observations of velocities themselves [23].

Optimization of Hyperparameters
The set of hyperparameters, θ = r ti , r yi , r xi , σ i , σ N | i = 1, 2 , which define the covariance function, K (θ), is determined through an optimization procedure which maximizes the probability of the hyperparameters conditioned on the observed velocities p θ | u d . We carry out the optimization through the application of the marginal likelihood method, explained in Section 5.4.1 of [25]. As [23] states, this is equivalent to maximizing the log marginal likelihood: where, B = K dd + σ 2 N I and n is the number of observations. To optimize the hyperparameters by maximizing the log marginal likelihood, we take the partial derivatives of the log marginal likelihood function (Equation 6), with respect to each hyperparameter θ j ∈ θ: where β = B −1 u d and Tr is the trace of a square matrix. We then utilize the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) gradient based optimization algorithm [32] to find the values of θ that maximize the log marginal likelihood. Any gradient-based optimization algorithm can be applied to Equation (7) to perform this step of the optimization [23].

Variations of the GPR Framework
Building upon the work by [23], we explore three different variations to the framework around the core Gaussian Process Regression (GPR) method, mainly involving differences in the preprocessing of the velocity data before the optimization of the hyperparameters and reconstruction is performed. We also present the results using only the core GPR method, as used by [23], with no variation (described here in Sections 2.2.1 and 2.2.2) to serve as a control for the experiment. The different GPR frameworks tested, using the LASER data in the Gulf of Mexico, are as followed: (1) Rotational: The coordinate system and drifter velocities are rotated such that the x-axis aligns with the mesoscale front making the spatial correlation scales parallel (x-direction) and perpendicular (y-direction) to the front. This is executed as followed: where u d r and v d r are the vectors of the observed u and v velocity components u d , v d that have been rotated by angle φ. x r and y r are the vectors of rotated x,y coordinates corresponding to u d r and v d r . The coordinate system and reconstructed velocities are rotated back to the original axes, after the optimization of hyperparameters and reconstruction are completed. φ, the angle of rotation, was determined by fitting a linear trend to the final position of all the drifters at the end of each reconstruction time window, such that the x-axis aligns with the front.
(2) Center of Mass: At every time step in the observations, the average of the drifter velocities and the center of mass of the drifters are subtracted from each drifter's velocity and position, respectively, before the optimization and reconstruction are performed. The velocity components and horizontal coordinates become: where u d i,j and v d i,j are the velocity component deviations from the mean drifter velocity u d j , v d j at each time step, j. In this framework, the mean velocities at every time step are acting as the prior mean of the observations u d in Equation (4) , however the subtraction occurs during the preprocessing stage before the regression is performed. x i,j and y i,j are the spatial deviations from the center of mass of the drifters x j , y j at every time step, n is the number of drifters at each time step, and m is the number of time steps used in the reconstruction. The optimization of the hyperparameters and reconstruction are then performed using u d i,j , v d i,j , x i,j , and y i,j . After the reconstruction, the same average drifter velocities and drifter positions are added back to the reconstructed velocity fields, and the coordinates there of, at each corresponding time step. In removing the center of mass, we aimed to reduce the effect of the background flow on the drifters' displacement, thus removing this signal from the regression.
(3) Low-pass Filtered GPR Prior Mean: The goal here was to separate the flow field into one large and one small scale component, and estimate each one separately. We attempt the optimization and reconstruction in two stages by first using low-pass filtered drifter velocities, meant to act as a theoretical larger scale component of the flow. For this, we use GPR with a 3-D squared exponential covariance function (Equation (5) with M=1). Afterwards, we estimate the smaller scale using the GPR-derived large scale component as the prior mean. Again this is done by subtracting the prior mean from the drifter velocities in the preprocessing stage. As done with the large scale component, we use GPR considering a 3-D squared exponential covariance function to calculate the small scale component of the flow. We then add these fields together to arrive at the final velocity field. We perform this method 3 times testing cutoff frequencies of the low-pass filter corresponding to periods of 12 h , 24 h, and 36 h, results of which did not vary significantly.
For all variations of the method's framework, the time window over which each optimization and reconstruction were performed was defined as an eight hour window, the middle six of which were used for the analysis. This was motivated by the initial comparison of the X-band derived and reconstructed GPR velocities, which showed relatively shorter time windows produced better agreement between the two. A buffer hour was added to the beginning and end of each 6 h window to ensure continuity between reconstructed velocity fields across adjacent time windows. All the time periods used for the analysis of the mesoscale front in the Gulf of Mexico, the submesoscale eddy in the Gulf of Mexico, and the mesoscale front in the Mediterranean Sea can be seen in Tables 1-3, respectively.

Velocity Validation with X-Band Velocity Fields
From the LASER experiment in the Gulf of Mexico, we utilize two X-band radar velocity fields which overlap with the drifter data in Figure 2, to validate the four different techniques of Gaussian Process Regression (GPR). The gridded, 475 m resolution X-band velocities are first interpolated, using a biharmonic spline interpolation, to the same horizontal resolution as the desired GPR grid. For each of the four GPR frameworks, we evaluate Equation (4) using the grid points of the X-band velocity field as the target velocities, u t , and the previously optimized hyperparameters for the given 8-hr period in which the X-band velocity field falls. We also keep the observed velocity vector, u d , unchanged, as to only rely on the drifter data for the evaluation of the target velocities. As a result that the drifters, and subsequent GPR-reconstructed velocity fields, have a temporal resolution of 15 min, we evaluate the X-band radar target velocities at the closest 15 min time step which corresponds roughly to the middle of each X-band velocity field seen in Figure 5. This procedure is performed twice, once during the convergence of the drifters onto the mesoscale front, corresponding to the X-band velocities in Figure 5a, and once after the formation of the submesoscale eddy occurs, corresponding to the X-band velocities in Figure 5b. This allows a one to one comparison of the X-band-derived and GPR-reconstructed velocities for each framework and flow regime. Based on the findings of the GPR velocity field validation performed using the Gulf of Mexico data, we select the best performing framework to apply to the drifter data along the mesoscale front sampled in the Mediterranean during CALYPSO. No data from the CALYPSO experiment is used in the velocity validation stage of this study.  Figures 6 and 7 show the velocity fields created using each of the Gaussian Process Regression (GPR) framework variations, evaluated at the X-band velocity field coordinates for the case of the mesoscale front and submesoscale eddy in the Gulf of Mexico, respectively. Following [23], we use the square root of the diagonal of the posterior covariance matrix, cov(u t ), as an estimate for the error in our reconstructed fields. Ref. [23] found the error estimated from the posterior covariance matrix to overestimate their calculated absolute errors by 0.01 m s −1 , showing one can make meaningful estimates of errors using this metric. The velocity fields shown in Figures 6 and 7, along with all other GPR velocity fields presented in this study, are masked to exclude any values associated with a posterior covariance error estimate greater than 0.04 m s −1 , for either the u or v velocity component.   Figure 6 shows that for the case of the converging mesoscale front, rotating the coordinates and velocities, such that the x-axis aligns with the front, results in the best qualitative agreement with the overlapping X-band velocities in Figure 5a. In addition, the rotational framework displays the largest area coverage for reliable velocity reconstruction whose error estimate is below the threshold. Results from the center of mass framework (Figure 6c) display poor coverage as estimated errors increase rapidly away from the location of the drifters. In comparison to the control method, where no variation is applied (Figure 6a), the utilization of optimized low-pass filtered velocities as a prior mean shows practically no change in the resulting field (Figure 6d). Figure 7 illustrates that for the case of the submesoscale eddy, some variations of the GPR framework can produce quite different results compared to that of the mesoscale front. The rotational framework (Figure 7b) no longer seems to have quite a drastic effect on the resulting velocity field compared to that of the mesoscale front case (Figure 6b). The center of mass framework shows some distinct effects, however it is difficult to say qualitatively if the resulting velocities are in better or worse agreement with the corresponding X-band velocity field (Figure 5b). One notable difference of the resulting velocity field is the less drastic spike in velocities in the southwest region of the eddy. Again the framework of using optimized low-pass filtered velocities as a prior mean (Figure 7d) shows almost no detectable deviation from the control (Figure 7a). Figures 8 and 9 show histograms of the velocity differences between the GPR-reconstructed and X-band-derived velocity fields (method described in Section 2.2.4). Figure 8 shows the velocity comparison for the case of the mesoscale front. Figure 8b highlights the extent of agreement between the X-band velocities and the velocities resulting from the rotational GPR framework, having RMSD values of 0.019 m s −1 for the u component and 0.047 m s −1 for the v component. Compared to the control method, this shows a good improvement in the estimate of the v component, while the estimates of the u component perform similarly. The effect of higher error estimates is again evident for the center of mass framework (Figure 8c), which still shows elevated RMSD for the data under the estimated error threshold. Changing the error mask threshold in order to utilize a similar amount of data points for the center of mass case, did not improve the performance of this framework. The use of low-pass filtered optimized velocities as a prior mean for the observation (Figure 8d) show little to no improvement upon the control method. In addition, the use of the rotational framework in unison with either the center of mass or low-pass filtered frameworks did not improve upon the results for using solely the rotational framework shown here.   Figure 7. Additionally listed are the Root Mean Square Differences (RMSDs) calculated between GPR and X-band velocity components. Only data associated with an error estimate less than 0.04 m s −1 are used for the calculations of RMSD and shown in the histograms. Figure 9 shows the quantitative velocity comparison between X-band and GPR-derived velocity fields in the presence of the submesoscale eddy. Compared to the control method, the rotation of axes and velocities results in an improvement in RMSD for the u velocity component, but an increased RMSD in the v component. The histograms for the control method (Figure 9a), rotational framework ( Figure 9b) and optimized low-pass filtered prior mean framework (Figure 9d), all exhibit long tails displaying absolute errors in both velocity components on the order of 0.5 m s −1 , which can be linked to the increased velocities seen in the southwestern region of the velocity fields in Figure 8a,b,d. The framework of removing the center of mass and average velocities results in significantly smaller deviations from the X-band velocity field, although some evidence of the same tail of increasing differences persists (Figure 8c). This framework however, still results in more agreement with the X-band velocity field, reducing the RMSD for the u component by more than half that of the control method, while not diminishing that of the v component. Applying the rotational framework, in addition to the center of mass or low-pass filtered prior mean framework for the submesoscale eddy case, does not result in improvement upon applying the center of mass framework alone, with regards to the RMSDs.

Comparison of GPR and Marine Radar Velocities in the Gulf of Mexico
Due to the performance of the rotational framework of GPR in the presence of drifters aligning on the mesoscale front, all further results corresponding to this dynamic regime, in the Gulf of Mexico and the Mediterranean Sea, will only be presented for this framework. Similarly, for the case of the submesoscale eddy in the Gulf of Mexico, only the results for the center of mass framework will be presented. Figure 10 shows the relative vorticity calculated from the interpolated X-band radar velocities and corresponding GPR velocity reconstructions used for the GPR velocity validation. The GPR-derived relative vorticity fields show more intense instances of relative vorticity, which is most evident across the mesoscale front. The GPR field shows the mesoscale front to be much more narrow with large values of positive vorticity concentrated along the front. The large values of relative vorticity shown in the GPR-derived fields are a direct result of the drifter's ability to sample the mesoscale front and submesoscale eddy with higher spatial resolution than that of the X-band radar which samples velocities at horizontal resolution of ∼475 m. Higher resolution measurements of submesoscale fronts have been shown to drastically increase horizontal velocity gradients and have shown these fronts to have widths as small as 50 m [21]. Large differences in velocity gradient-based measurements are expected between the X-band and GPR velocities due to differences in sampling strategy, any small phase shifts of the flow features sampled, and the natural tendency for gradient-based calculations to have larger deviations than between the original velocity fields. For these reasons, we do not consider the absolute differences of velocity gradient-based quantities between the X-band and GPR velocities, but rely only on the velocity statistics for comparison. Figure 10. Relative vorticity normalized by the planetary vorticity f , calculated as ζ = dv dx − du dy / f , using (a) X-band interpolated velocities and (b) GPR-reconstructed velocities using the rotational framework, associated with the velocity validation for the case of the mesoscale front. Additionally shown are fields of relative vorticity calculated using (c) X-band interpolated velocities and (d) GPR-reconstructed velocities using the center of mass framework, associated with the velocity validation for the case of the submesoscale eddy. Only data points associated with a posterior covariance error estimate less than 0.04 m s −1 are shown. Origins of (a,b) are that of Figure 6 and origins of (c,d) are that of Figure 7.    Figure 11 shows two snapshots of the velocity fields created using the rotational framework of GPR (i.e., the optimized hyperparameters form Table 1) for the mesoscale front in the Gulf of Mexico. Additionally shown are corresponding fields of relative vorticity and divergence, calculated as ζ = dv dx − du dy and δ = du dx + dv dy , respectively. The convergence of the drifters onto the front is clearly observed between the time steps shown. From Figure 11a,b, it is evident that as the drifters converge, the velocity shear across the front quickly intensifies. This is also reflected in the fields of relative vorticity, shown in Figure 11e,f, which depict a narrow band of strengthening positive vorticity along the front. Accompanied by this increase in vorticity, is the presence of strong convergence (Figure 11c,d), which is also evident given the collection of the drifters along the front. Animations of all divergence and relative vorticity fields used for the analysis of the Gulf of Mexico mesoscale front can be found in Supplementary Materials S2 and S3, respectively.
Transects of potential density anomaly, σ θ , salinity, and temperature along the two ship tracks shown in Figure 2a, which cut across the mesoscale front during the time period of the corresponding GPR reconstruction, are plotted in Figure 12. Figure 12 shows a strong front that is dividing cold, fresh Mississippi outflow water from the warm, salty waters of the interior Gulf of Mexico. During these hours, when the drifters converge on the mesoscale front, we see evidence of vertical mixing, predominately on the cold (west) side of the front. The cold, fresh surface water can be seen mixing with the warmer, saltier water directly below it down to about 75 m, where it even seems to be mixing with the waters just below the mixed layer, mostly evident in the transects of temperature (Figure 12e,f) and salinity (Figure 12c,d). There also seems to be some degree of downwelling occurring below the surface signature of the front, evident in the depression of the deeper pycnocline and deep mixing directly under the drifters captured in the front. In contrast, the warm (east) side of the front seems to be already well mixed at the time of the measurement, showing almost uniform temperature and salinity from the surface to the bottom of the mixed layer. From the first transect which begins at 07:25 UTC to the next which begins on 09:23, the mesoscale front intensifies, becoming very sharp at the surface and spanning the entire depth from the surface down to the pycnocline at ∼80 m.

Submesoscale Eddy: Correlation Scales, Kinematics, and CTD Transects
All the optimized hyperparameters for the five time intervals of reconstruction associated with the submesoscale eddy in the Gulf of Mexico are shown in Table 2. Every optimization procedure, across the five time windows results in optimized noise, σ N , with values of 0.002-0.007 m s −1 . Each optimization also found 2 distinct scales of motion, one smaller and one larger, having a smaller scale standard deviation of velocity ranging from 0.001-0.067 m s −1 and a larger scale standard deviation ranging from 0.063-0.187 m s −1 . Although every optimization period shown in Table 2, produced distinct scales in the zonal, meridional, and temporal dimensions, not all the correlation scales seem to contribute significantly to the velocity fields produced. Four of the larger time correlation scales (r t2 ) listed are very large numbers, on the order of thousands of hours or greater. There are also two values of the larger meridional correlation scales (r y2 ) for the v velocity component that are on the order of thousands of km. These very large correlation scales, being in the denominator of the terms in the squared exponential covariance function, cause these terms to have little influence on the reconstructed velocities. Physically, this can be interpreted as representing the larger scale motions of the flow that will not change over the short distances and time periods being investigated here. Table 2. Optimized hyperparameters corresponding to the drifter velocities within the submesoscale eddy in the Gulf of Mexico for the five time windows chosen. Optimization was performed after average velocities and center of mass of the drifters were subtracted from the velocity and position of each drifter, respectively, at each time step. ** The final time window was shifted by one hour to correct for unrealistic velocities present in the final reconstructed velocity fields for the original time window. Excluding these outliers, the larger time correlation scale range is ∼3-10 h, while the smaller time correlation scale range is ∼0.5-1.5 h. The larger meridional and zonal correlation scale range is ∼200 m-3 km, with the meridional scales having a notable intermediate outlier of ∼12.5 km for the optimization period beginning at 18:00 UTC on 29 January. The smaller spatial correlation scales for u and v almost mirror one another within each time window, ranging from ∼30-800 m. These values of r x1 and r y1 also decrease for both the u and v velocity components from the second time window to the last, which seems appropriate given the drifters convergence upon the center of the submesoscale eddy. The smaller time correlation scales (r t1 ) for both velocity components, generally trend towards smaller values (r t1 < 1 h) in the later time windows as well.
Six snapshots of the GPR-reconstructed velocity fields of the Eddy are shown in Figure 13. Moving from Figure 13a-f, we follow the initial creation of the submesoscale eddy in the form of an instability along the mesoscale front. As the eddy begins to spiral it is also being advected by the strong background flow in the region, which is why the velocity vectors seen in Figure 13 never show the closed circulation of the eddy. The animation of the drifter motion, found in S1 of the Supplementary Materials, shows the formation and advection of the spiraling eddy much more clearly. The signature of the eddy, however, is evident in the fields of velocity magnitude (Figure 13) showing the portion of the eddy's rotation which aligns with the larger scale background flow to be much greater than that opposing the direction of advection.
Fields of divergence and relative vorticity of the same panels in Figure 13 are plotted in Figures 14  and 15. From the fields of divergence, it is evident that strong episodes of surface convergence is measured as the drifters spiral around the eddy. Beginning in a line along the mesoscale front, stretched over ∼6 km, the drifters themselves converge to fit within a 10 m by 18 m box over a time period of ∼29 h (see S1, S4, or S5), exemplifying the strong convergence present in the eddy. Associated with these high levels of convergence in the eddy, are even larger values of relative vorticity, seen in Figure 15, increases of which seem to coincide with the episodic convergence events.
Animations of all the divergence and relatively vorticity fields for the submesoscale eddy can be found in Supplementary Materials S4 and S5, respectively.    Figure 16 shows the MVP data measured in the vicinity of the eddy. The first transect in time (Jan 30th 01:30-02:24 UTC) is recorded a few hours before the drifters cross the transect. Although not perfectly symmetric, the sloping lines of constant σ θ on either side of the drifter positions located in the center of the eddy in Figure 16a, seem to suggest that the eddy is penetrating the plane of the transect. This can also be seen very clearly in the salinity data in Figure 16c. On the western side of the salinity transect, there is a large blob of high salinity water, which has been subducted underneath the fresh, cool Mississippi outflow waters. Underneath this blob there is also evidence of strong mixing between the two water masses, likely a result of the spiraling eddy forcing the subduction of the heavy saltier water, generating turbulent mixing. Similar structures can also be seen in the temperature transect (Figure 16e). To a lesser extent, horizontal advection around the core of the eddy throughout the depth of the mixed layer, likely associated with some mixing, can also be seen in these transects.   Figure 16b shows the core of the eddy extending from the surface down to ∼80 m. Similar indicators of subduction and turbulent mixing can be seen in both western sections of the salinity and temperature transects, along with an enhanced signature of horizontal transport and mixing between the two water masses as they circulate around the eddy core. As this instability eddy forms, it is not only forcing surface water to converge but also acting to stir and mix water at depth around the core of the eddy. In addition, there appears to be some degree of upwelling occurring at the base of the eddy core, sucking water up from below mixed layer, visible in the transects of σ θ , salinity, and temperature. The presence of upwelling is also supported by the surface divergence observed at the core of the eddy, several hours later, in Figure 14f.

Comparison to Mesoscale Front in the Mediterranean
The drifter trajectories used for the GPR reconstruction in the Western Mediterranean Sea are shown in Figure 3, along with an SST image from the Aqua/MODIS satellite. From the SST data, a mesoscale feature, known as the Western Alboran Gyre (WAG), can be seen driving the drifter trajectories towards the southeast. Within the gyre, circulation appears to be some streaks of higher SST, which could be evidence of submesoscale filaments within the larger scale flow. Looking at the trajectories, there does appear to be an alignment of the drifters along the frontal edge of the WAG, similar to the alignment of the drifters in the presence of the mesoscale front in the Gulf of Mexico. For this reason, the rotational framework of GPR is employed to reconstruct the velocity field around the drifters over two 8-h Figure 17 shows two snapshots of the reconstructed velocity fields in the Mediterranean, plotted over fields of velocity magnitude, divergence, and relative vorticity. From these two snapshots, there appears to be some convergence of the drifters onto the front. The velocity magnitudes are generally larger in the WAG than for the mesoscale front in the Gulf of Mexico (Figures 11 and 17), and the velocities in the WAG do not show a sharp velocity shear across a defined frontal feature. The fields of divergence (Figure 17c,d) and relative vorticity (Figure 17e,f) also display less organized patterns of these calculated quantities, which are also considerably smaller in magnitude. The convergence of the drifters themselves and local maxima within the relative vorticity and divergence fields suggest the presence of some submesoscale variability within the flow, possibly being driven by the formation of filaments as water masses are stirred around the gyre. Animations of divergence and relative vorticity fields over all the time steps used for the analysis of the Mediterranean front can be seen in Supplementary Materials S6 and S7, respectively. Two cross sections of density (σ), salinity, and temperature across the mesoscale front in the Mediterranean are shown in Figure 18. In all cross sections evidence of a sharp, slanted pycnocline can be seen at the bottom of the mixed layer varying in ∼100-150 m depths. In Figure 18a,c,e a surface front can be seen extending down to about 50 m where it seems cooler, saltier water on the north side of the front is being subducted underneath the warmer waters to the south. This movement and subduction of cool, salty water into the WAG is consistent with turnstile lobes produced in a model by [33]. Fluid transported in lobes of this type have been shown to elongate into filaments as they wrap around edges of coherent structures [34]. Figure 18b,d,f show that the convergence along this front may be the result of alternating filaments of different salinity and temperature within the mixed layer of the gyre circulation. The SST image shown in Figure 3 shows a similar pattern, depicting small areas of warmer surface waters just to the south of the drifter trajectories. Evidence of subduction can be seen in these cross sections, with enhanced mixing taking place in the waters below the filament, but above the deeper pycnocline. The normalized histograms of divergence and relative vorticity for all three of the investigated flow features are shown in Figure 19. The histograms are created using all the velocity fields within the 6 h analysis windows, accumulated over the entirety of each event. Only velocities associated with posterior error estimates, under the previously stated threshold, are used. The histograms relating to the mesoscale front in the Gulf of Mexico (Figure 19a,b) show considerably skewed distributions of divergence and relative vorticity, the former being strongly skewed towards negative values, and the latter towards positive values. The maximum value of positive vorticity along this front is ∼34 f , while the minimum divergence, or convergence, has a magnitude of ∼10 f . The histograms of the submesoscale eddy divergence (Figure19c) and relative vorticity (Figure 19d) show similar patterns of skewness, but with more robust quantities of negative divergence and positive relative vorticity. While having a larger maximum magnitude of negative divergence (∼20 f ), compared to that of the Gulf of Mexico mesoscale front, the submesoscale eddy has a slightly smaller value of maximum relative vorticity ∼30 f . Compared to the features sampled in the Gulf of Mexico, the front in the Mediterranean shows little to no skewness in either the histograms for divergence (Figure 19e) or relative vorticity (Figure 19f). The magnitudes of the minimum and maximum values are also considerably smaller.

Vertical Velocities
Through the conservation of mass we can apply the continuity equation: where dw dz is the vertical derivative of vertical velocity, w, to estimate order of magnitude vertical velocities from our calculated values of divergence. By integrating over depth, z, we can write: where w z is the vertical velocity at the surface and w z+∆z is the vertical velocity at depth z + ∆z. From the available CTD transects for all three flow features sampled, we can estimate a depth at which w z+∆z ≈ 0. We choose to define this depth of approximately zero vertical velocity at the deepest extent of continuous, uniform salinity and temperature observed from the surface down to a depth ∆z. Doing so, we estimate depths of zero vertical velocity equal to ∼40 m, ∼40 m, and ∼50 m for the case of the mesoscale front in the Gulf of Mexico, the submesoscale eddy, and mesoscale front in the Mediterranean, respectively. Estimated values of w z , based on these depths of zero vertical velocity and maximum convergence for each case, are shown in Table 4. These rough estimates of vertical velocities associated with the sampled features are meant to provide order of magnitude estimates of w z , given the maximum observed values of surface convergence. From the values of w z calculated, we can conclude that maximum vertical velocities on the order 10 −2 m s −1 are present in both mesoscale fronts in the Gulf of Mexico and the Mediterranean. These estimates also suggest that the large convergence observed within the submesoscale eddy in the Gulf of Mexico are creating maximum vertical velocities up to approximately twice that observed in either frontal feature. Table 4.
Vertical velocities estimated using Equation (11), prescribed depth of zero velocity, and corresponding maximum convergence observed in each flow feature.

X-Band and GPR Velocities in the Gulf of Mexico
The values of RSMD between the velocities produced with X-band radar and the velocities produced with the best performing GPR frameworks are comparable to those calculated by [20]. Utilizing all coinciding pairs of drifter velocity and X-band velocity measurements, Ref. [20] found values of RSMD equal to 0.035 m s −1 for the u component and 0.04 m s −1 for the v component, calculated over all the available data from the LASER campaign. In the case of the rotational GPR framework applied during the convergence of the drifters onto the Gulf of Mexico mesoscale front, the RSMD for the u and v velocity components, listed in Figure 8b, are in good agreement with the findings of [20]. A large portion of the difference between drifter and X-band velocities is known to stem from the difference in the effective depth of each measurement [20]. Previous studies have shown wind-and wave-induced currents produce considerable vertical shear in the upper few meters of surface currents in the open ocean [35,36]. The histograms of absolute velocity difference (Figure 8b) shows the GPR-derived v component to be skewed negatively, and the u component slightly skewed positively, compared to the X-band velocities. Considering the wind magnitude (7-8 m s −1 ) and wind direction (traveling SSE 150 • ), shown in Figure 4, at the time of this GPR validation, one would expect a wind-driven velocity shear in the vertical profile of the surface current, which could explain the skewness of absolute difference in both velocity components.
The center of mass GPR framework performed best for the case of the submesoscale eddy in the Gulf of Mexico with relation to the RMSD between GPR and X-band velocities, but still resulted in higher RMSD values of up to 0.07 m s −1 . The wind during this period of validation was ∼4 m s −1 , traveling ENE at ∼75 • , and again the histograms of absolute velocity difference (Figure 9c) show skewness that can be at least partially explained through wind-induced vertical shear of the surface current, however the relatively low wind magnitudes may suggest another mechanism is affecting the vertical structure of the surface current in the upper few meters. The radar backscatter image shown in Figure 5b shows the outline of multiple surface fronts swirling around what appears to be the center of the submesoscale eddy. The strong submesoscale eddy is likely responsible for the generation of complex turbulent flow conditions with variability that could not only be affecting the vertical profile of the surface currents, but also could be under resolved in the X-band radar velocities. As previously stated, the drifters can also sample the flow at higher spatial resolution depending on the spatial distribution of the drifters. The GPR reconstruction also takes into account drifter data before and after the X-band velocity analysis time period, which could also cause increase the velocity difference between the methods. In addition to having the best agreement with the X-band-derived velocities, the center of mass GPR framework seems to limit the area of influence each observation has on the reconstructed velocity field, which makes this framework the most conservative estimate out of the GPR variations tested.
The center of mass approach would be ideal if the center of mass coincided with the center of the eddy, so that when removed, the drifters would be circling around the origin, and the effect of the larger flow would be completely removed from the regression. As the drifters are often not perfectly arranged around the eddy center, removing the center of mass will not always perfectly remove the effects of the larger scale, which could also add to the relatively larger RMSD for this case, compared to that of the rotational framework applied to the mesoscale front.

Gulf of Mexico vs. Mediterranean Flow Features
The SST images in Figure 2 show the development of the mesoscale front in the Gulf of Mexico as a large discharge event from the Mississippi River sends water traveling offshore, which elongates to the southwest and northeast. Figure 2b shows the SST image produced for 00:00 UTC on 31 January. Even though this 24 h average of SST is centered several hours after the drifter observations of the submesoscale eddy in the Gulf of Mexico, it is clear to see that the final drifter positions in Figure 2b show the drifters at the center of a spiraling submesoscale instability which has formed along the mesoscale front. Additionally seen in Figure 2b, along the mesoscale front are what appears to be the signature of other submesoscale instabilities, potentially similar in structure to the drifter sampled eddy. It appears that multiple instabilities developed simultaneously along this mesoscale front, resulting in the rolling wave-like pattern seen in the SST. In comparison, Figure 3 does not reveal the same type of submesoscale signature in the Mediterranean SST data. The movement of the drifters in the Mediterranean seems to be dominated by the larger scale movement associated with the northern front of the WAG. The WAG is known to be the strongest dynamical feature in the region, reaching velocities of up to 1.5 m s −1 , spanning 100 km across and reaching a depth of 200 m [37]. After the analysis period used here for the Mediterranean data, the drifters begin to disperse as they travel around the eastern curve of the WAG, showing little signs of submesoscale activity.
The difference in the scales of motion captured by the Gulf of Mexico and Mediterranean drifter data is also apparent in the optimized hyperparameters. This is perhaps most pronounced in the hyperparameters associate with larger scale motions found in the Mediterranean data. The optimization found both the u and v velocity components to have larger horizontal correlation scales on the order of a couple 10's of km, corresponding to correlation times scales of ∼4-7 h. This set of larger correlation scales shows no tendency of the flow to be correlated differently in the along-front and cross-front directions, suggesting an absence of strong frontal velocity gradients in the large scale motion. Both sets of the horizontal correlation scales for the Gulf of Mexico front show significantly larger correlation scales in the along-front, than in the cross-front directions for both velocity components, which can be linked to the strong velocity shear in the cross-front direction seen in Figure 11. To some extent the smaller horizontal correlation scales in the Mediterranean do show similar cross-front and along-front tendencies, which may be another indication of convergence occurring due to submesoscale filaments aligned with the larger frontal structure. However the front in the Gulf of Mexico, being associated with much smaller cross-front correlation scales (as small as 70 m), displays a more extreme case of cross-front velocity shear.
In the case of the submesoscale eddy, there are four instances where the larger time scale for either the u or v components, and two instances where the meridional scale for the v component, are exceptionally large. As this submesoscale eddy is being transported SSW by a strong larger scale current and the eddy itself is very well organized, it is rational that the individual velocity components could be described without utilizing all 6 correlation scales during each time period, especially those pertaining to the larger scales of time and meridional distance. In other words, the translational velocity of the eddy, along with the velocities which make up the eddy's circulation, could be superimposed in a way that only one time correlation scale, and/or one meridional correlation scale in the case of the v component, dominates the flow. As previously mentioned, the center of mass framework does not always perfectly remove the affects of the larger scale flow due to the irregular distribution of the drifters around the eddy center, which could also explain the asymmetry found in the larger correlation length scales of the v component (large r y2 and small r x2 ). It seems the larger correlation scale often captures the larger background flow while the smaller scales are more consistent with the eddy itself.
The optimization of the hyperparameters allows the velocity measurements themselves to determine the scales of motion captured by the drifters. When comparing GPR based optimized correlation scales to those derived empirically from the same drifter velocity data, Ref. [23] found reasonable agreement between the two, with the spatial and temporal correlations of each showing similar correlation decay. There is the fact that as the drifters align on either front or converge towards the center of the eddy the dynamics sampled can also become limited. Since the optimized correlation scales are dependent on the drifter velocity data, they are also limited by the spatial coverage of the drifters. This could partially explain why the larger scale correlations for the mesoscale front in the Gulf of Mexico are much smaller than one would expect when analyzing a mesoscale feature. Nonetheless, the correlation scales revealed by the optimization seem congruous to what is observed in the raw drifter data. From the tight alignment of the drifters on the mesoscale front in the Gulf of Mexico, along with the signature of this front seen in the radar backscatter data shown in Figure 5a, we expect to see short cross-front, and larger along-front, correlation scales like those found in Table 1. The correlation scales found from the Mediterranean drifter data is analogous to what is observed in the drifter motion itself, with the flow being characterized by two very distinct scales of motion: the dominant mesoscale flow carrying most of the drifters in an almost uniform fashion along the front and localized smaller scale movements causing a few of the drifters to converge. The interpretation of the correlation scales related to the submesoscale eddy has the added complexity that the optimization takes place after the center of mass and average velocity is removed at each time step, but still produces reasonable spatial and time scales representative of the movement of the drifters through the analysis periods.
The histograms of divergence and relative vorticity shown in Figure 19 support the findings of the optimized hyperparameters for the different data sets. The histograms of the front in the Gulf of Mexico (Figure 19a,b) and the front in the Mediterranean (Figure 19e,f) show the different regimes of flow captured by the respective drifter data sets. Most notable is the large tail of positive relative vorticity seen in the histogram for the Gulf of Mexico front, which reaches a maximum of ∼34 f . This large maximum in relative vorticity is three times that found in the Mediterranean front, highlighting the strong cross frontal velocity shear present in the Gulf of Mexico case. The divergence histograms for the two fronts also reveal large differences in the dynamics and water masses at play for each case. The Gulf of Mexico front, being heavily skewed toward negative divergence, is dominated by a convergent flow field, which is a direct result of the subduction of the heavier interior gulf waters underneath the Mississippi River outflow waters. The uniform histogram of divergence for the Mediterranean front suggests that the flow is relatively non-divergent, suggesting either a lack of strong frontal behavior between two water masses of significantly differing densities and/or the presence of dominating mesoscale dynamics which overpowers submesoscale instabilities. Similar dynamic regimes have been tested using multi scale LES by [4], which showed deeper baroclinic instabilities tend to dominate over surface submesoscale flows within a shallow surface layer. The sharp slanted pycnocline present at the bottom of the mixed layer under the Mediterranean front, shown in Figure 18, reveals the presence of larger scale forcing. While not necessarily baroclinic in nature, this deeper instability is likely dominating over the smaller mixed layer features observed. The transects showing the vertical structure of the Gulf of Mexico front ( Figure 12) show a relatively uniform pycnocline around ∼80 m depth, underneath the front which is contained in the mixed layer above, possibly revealing a lack of larger scale forcing that would otherwise dominate the flow.
The histograms pertaining to the submesoscale eddy (Figure 19c,d) are heavily skewed, being dominated by convergence and positive vorticity. The magnitudes of these quantities are among the highest recorded using drifters to date, reaching ∼20 f in convergence and ∼30 f in relative vorticity. In comparison to previous studies, results of which are shown in Table 5, which have analyzed submesoscale structures using drifter data, the magnitudes of divergence and relative vorticity found here are ∼2-3 times and ∼3-4 times greater, respectively, than previously reported measurements [5,8,23]. These previous studies analyzed submesoscale structures that formed along coastal outflows, away from mesoscale influence. The extraction of energy from the mesoscale feature observed here could possibly be the cause of the increased magnitudes of convergence and relative vorticity. We also found horizontal correlation scales an order of magnitude smaller than [23], on the order of 10s of meters, for both the mesoscale front and submesoscale eddy, indicating increased ageostrophic motion associated with smaller scales. Ref. [21] which sampled a front 30-50 m in width from high resolution sea surface roughness amplitude data, found current gradients ∼2-4 times greater than the maximum values of divergence and relative vorticity reported here. This highlights the impact that sampling velocity fields at higher spatial resolutions has on these calculations. Table 5. Synthesis of recent studies on submesoscale kinematics. Note: All maximum absolute values listed are a result of positive relative vorticity and negative divergence, with one exception denoted by (*) in which the maximum absolute value of divergence was derived from a positive quantity. ** Reported only across front current gradients ( du dx ).

Reference
Max The vertical velocities estimated here, using Equation (11), are comparable to those measured by [5] using an ADCP equipped Lagrangian float, which followed waters subducted within a cyclonic frontal eddy. The submesoscale eddy sampled by [5] during the LASER campaign was in close vicinity to, and observed about one week after, the Gulf of Mexico flow features shown here. Both the CTD cross sections shown in Figures 12 and 16 and the cross section of potential density shown by [5] show very similar mixed layer depths and cross frontal structure. Furthermore, the Lagrangian float used by [5] returned to the surface, in part, due to weak upwelling after reaching a depth of about 35 m, which supports the 40 m depth of zero vertical velocity prescribed here for the Gulf of Mexico features. Ref. [5] reported negative vertical velocities at the surface convergence zone of 1-2 cm s −1 , with maximum convergence magnitudes reaching 6 f . Given the similar values of convergence observed along the mesoscale fronts in the Gulf of Mexico and Mediterranean Sea in this study, the corresponding downward vertical velocity estimates of up to almost 3 cm s −1 are congruous with the previously reported findings. The downward vertical velocities found in the submesoscale eddy reach up to 5.6 cm s −1 , showing the effects strong convergence has on vertical velocities in the upper ocean. Vertical velocities of this magnitude are quite large compared to previous estimates on the order of 10 −4 m s −1 [38]. The large estimated vertical velocities and skewness of divergence histograms shown in Figure 19 reflect the large quantity of vertical exchange observed in the Gulf of Mexico. It is noteworthy that we find similar magnitudes of vertical velocities along the mesoscale front in the Mediterranean, however the observed vertical exchange associated with this front is considerably smaller, being limited to the localized regions of convergence captured by the drifter data.
Examining the snapshots in Figures 14 and 15, along with S4 and S5 of the Supplementary Materials, there seems to exist intermittent episodes of large convergence, which coincide with the convergence of surface fronts being swirled around the eddy. Examples of these surface fronts can be seen in the backscatter image shown in Figure 5b, where several bright lines can be seen bending around, and meeting in middle of the eddy center. In the time steps following this image of radar backscatter, some of the drifters can be seen traveling along these fronts, causing them to seemingly collide with other drifters as they also converge onto the same spiraling line. These swirling structures appear similar to those observed by [5], in which these structures were referred to as zippers.
Another interesting feature of the eddy is highlighted by Figures 14f and 16b. Although the observations of the drifters in Figure 14f are plotted several hours after the density anomaly transect in Figure 16b, they both show the presence of a defined eddy core where some small amount of divergence seems to be occurring. The surface divergence at the eddy center appears to be connected, without interruption, all the way to the upwelling seen at the bottom of the mixed layer. This would imply waters from below the deeper pycnocline are being sucked well into the mixed layer, possibly even reaching the surface over the lifetime of this eddy. This shift from convergent to divergent regimes felt by the drifters as they approach the eddy center during the final time period in Table 2, could potentially be the cause of the difficultly to reconstruct realistic velocities over the original time window. The supplementary animation of the eddy divergence (S4) shows that the drifters briefly disperse as this surface divergence occurs, but then quickly converge again either because the drifters have moved away from the divergent center of the eddy or because the surface divergence and upwelling may be intermittent as well. It is not clear from the observations the consistency of the upwelling and surface divergence, but from the time difference between the MVP transect and surface signature of this divergence it seems to be a reasonably stable feature of the eddy.
The strong velocity shear across the front in the Gulf of Mexico seems to be a very important forcing mechanism in the development of the submesoscale eddy sampled by the drifters, and analogously in the other submesoscale instability eddies seen in the SST image (Figure 1b). The highest values of positive relative vorticity are also calculated before the development of the submesoscale eddy, which could be an indication of energy dissipation by the submesoscale eddy. Just before the formation of the submesoscale instability, sloping lines of constant σ θ can be seen on the western side of the front in Figure 12b. However, to the east of the front lies a uniform mass of water having the same density almost through out the entire mixed layer. The sloping lines of constant potential density seen here are clear signatures of baroclinic instability, as shown by [1,4], which most likely acts a forcing mechanism in the creation of the submesoscale eddy. The straight vertical lines of constant σ θ , salinity and temperature which define the edge of the eastern water mass seen in Figure 12b,d,f may indicate that the velocity shear seen across the front at the surface is present to some extent through out the mixed layer directly along the front. This would mean the across front velocity shear is playing a major role in forming the instability we observe. The submesoscale eddy also extends to the very bottom of the mixed layer, likely being driven by forcing mechanisms which extend across the same depths. The similarity in the instabilities observed in the SST in Figure 1b implies that these too are forced by strong cross frontal velocity shear.
The overlying limitation to the GPR method used here is the available drifter data. As drifters begin to converge upon one another, either onto a front or towards the center of an eddy. The coverage of the velocity fields can become very limited, making it necessary to assess the increase in error away from the observations. In the case of the Gulf of Mexico front, due to its increased convergence, once the drifters align on a single line, it is possible that the reconstruction is missing some of the final development stages of the front before it becomes unstable and begins shedding the submesoscale eddy. More work, with increased observations may be needed to further explore this phenomenon. In the case of the eddy, the same problem exists when trying to reconstruct velocity away from the eddy center. Figure 7 illustrates this difficulty, as each GPR variation results in the highest velocities away from any available drifter observations, making it challenging to report on the validity of these values. Figure 14e also highlights the limitation of the data, as it depicts an area of strong convergence in a region where there is not a strong presence of drifters sampling the flow. This could be an artifact of the reconstruction, due to the large values for the optimized time and meridional length scales and/or the aliasing of finer scale motions in between individual drifters.
Comparing the maximum magnitudes of divergence and relative vorticity during the time period corresponding to the submesoscale eddy utilizing each of the different GPR frameworks tested, we find a wide range of values. The rotational GPR framework produced maximums of relative vorticity ∼50 f and convergence ∼30 f , while the control and low-passed GPR prior mean frameworks produced maximums of ∼35 f for vorticity and ∼30 f for convergence. Although the maximums differ by a significant margin in some cases, the overall skewness of the resulting histograms from the different frameworks are very similar to that of the center of mass framework presented here. As a result that the removal of the center of mass from the drifter observations inside the eddy seems to limit the area of influence around each drifter, this method produced the most confident results, but could however be missing some key information due to a lack of observations away from the eddy center.

Conclusions
Through the utilization of available marine X-band radar velocities, we have developed and tested new frameworks of GPR velocity reconstruction and have found certain frameworks to excel for different flow regimes, with none performing best for all cases. Further development and validation of GPR methods will be needed in the case of future studies. Given the current state of available data and methodology, the results shown here are among the most robust estimates of kinematic statistics for the flow regimes focused on during this study.
We have analyzed and quantified 3 different regimes of flow: The unstable mesoscale front in the Gulf of Mexico, the subsequent submesoscale eddy which forms from said instability, and the mesoscale front in the Mediterranean which does not become unstable, but still contains submesoscale filaments responsible for some convergent phenomena. The skewness, or lack there of, in the relative vorticity and divergence distributions calculated for the different flows, helps in the classification of the dominant dynamics in each case, as mesoscale dominated flows are expected to be mostly non-divergent (i.e., δ/ f histogram with little skewness). The front sampled in the Gulf of Mexico, having values of relative vorticity almost 3 times that of the Mediterranean observed front, could give key insights to the unstable flow conditions necessary for the generation of submesoscale eddies along mesoscale fronts. In addition, the production of the submesoscale eddy itself significantly impacted the quantity and magnitude of convergence measured in the region, again reaching almost 3 times that caused by the ageostrophic submesoscale filament associated motions observed in the Mediterranean frontal gyre circulation. As a result of this strong convergence we find vertical velocities associated with the submesoscale eddy, twice that found in either mesoscale front. The submesoscale instability eddy showcases how strong instances of relative vorticity can cause considerable surface convergence and subsequent vertical exchanges in the water column through the lateral mixing of water masses across mesoscale fronts.
Supplementary Materials: The following are available online at http://www.mdpi.com/2311-5521/5/3/159/s1, Video S1: Drifter trajectories corresponding to the submesoscale eddy velocity reconstruction in the Gulf of Mexico, Video S2: Divergence and GPR-reconstructed velocity fields, along with drifter positions, for the mesoscale front in the Gulf of Mexico, Video S3: Relative vorticity and GPR-reconstructed velocity fields, along with drifter positions, for the mesoscale front in the Gulf of Mexico, Video S4: Divergence and GPR-reconstructed velocity fields, along with drifter positions, for the submesoscale eddy in the Gulf of Mexico, Video S5: Relative vorticity and GPR-reconstructed velocity fields, along with drifter positions, for the submesoscale eddy in the Gulf of Mexico, Video S6: Divergence and GPR-reconstructed velocity fields, along with drifter positions, for the mesoscale front in the Western Mediterranean Sea, Video S7: Relative vorticity and GPR-reconstructed velocity fields, along with drifter positions, for the mesoscale front in the Western Mediterranean Sea.