The data used in this study originate from three types of drifters: the CARTHE, CODE, and SVP designs.

The Consortium for Advanced Research on Transport of Hydrocarbon in the Environment (CARTHE) drifters measure the surface currents within 60 cm of the sea surface. They are developed to be compact, easy to transport and assemble, and 85% biodegradable (Novelli et al., 2017) so that very large deployments can be attempted in the ocean while being eco-friendly. They have low windage, thanks to the flat toroidal floater, in addition, their motion is not affected by rectification caused by wind waves, thanks to the flexible tether connecting the floater and the drogue. CARTHE drifters slip velocity is estimated to be less than 0.5% for winds up to 10 m/s (Novelli et al., 2017).

The Coastal Ocean Dynamics Experiment (CODE; Davis, 1985) drifters measure the currents within the top meter of the water column. They are composed of a 1-m-long negatively buoyant tube with four drag-producing sails extending radially from the tube and four small spherical surface floats to provide buoyancy (Poulain, 1999). Poulain and Gerin (2019) demonstrated that CODE drifters follow the currents with an accuracy of about 3 cm/s, even under strong wind conditions. The wind and wave- induced slippage was estimated to be 0.1% of the local wind speed. If compared to drogued SVP drifters, CODE drifters are estimated to move downwind (and to the right of the wind) by about 1% of the wind speed, due to slip, Ekman and Stokes drift currents (Poulain et al., 2009).

The Surface Velocity Program (SVP) drifters are the standard drifters of the Global Drifter Program (Niiler, 2001; Centurioni, 2018). They consist of a spherical surface buoy tethered to a weighted drogue that allows it to track the horizontal motion of water at a nominal depth of 15 m. Below the surface buoy, a tether strain gauge measures the tension of the buoy-drogue connection to monitor the drogue presence and a thermistor measures sea surface temperature. Measurements of the water-following capabilities of the SVP have shown that they follow the water near 15 m depth to within 1 cm/s in 10 m/s winds (Niiler et al., 1995) when the drogue is attached.

It is important to note that the CODE/CARTHE, and SVP drifters measure different dynamics, being, respectively, in the first meter and at 15 m below the surface, especially in case of strong wind conditions. Indeed the near-surface drifters are more directly affected by the winds and waves and sample the Ekman wind-driven currents, Stokes drifts and possibly Langmuir cell circulation. The SVP drifters, on the other hand, are designed to measure the currents at 15 m depth which are less influenced by the above-mentioned dynamics forced by the winds and waves, while responding to wind-induced Ekman currents at the specified depth of the drogue.

A number of aspects regarding near surface effects on the sampling of surface drifters have been investigated in the literature. During open-air tank experiments to study surface transport induced by steep waves, Novelli et al. (2020) found that CARTHE drifters (at 60 cm depth) remained stationary, i.e., insensitive to the wave induced surface current. On the other hand, tracers and very surface floaters (1–5 cm depth) showed similar advection along the direction of wave propagation. The drastic difference corroborates the rapidly depth-decaying near-surface shear induced by Stokes drift and wind stress (Novelli et al., 2020). The different drift behavior observed in controlled tank experiments contributed explaining the very low landing rates of CARTHE/CODE drifters compared to very surface floaters (1–5 cm depth), deployed during massive releases in the Gulf of Mexico (D'Asaro et al., 2020). The combination of near-surface current shear, induced by winds and high-frequency wind-waves, together with near-shore wave-induced drift (via Stokes drift and wave-breaking) was identified as a potential major mechanism leading to cross-shelf transport intensification and to the stranding of tracers and very surface floaters (1–5 cm depth) (Novelli et al., 2020). As for the Langmuir circulation regime, Chang et al. (2019) compared CARTHE drifters' dispersion to very surface floaters (1–5 cm depth). At scales (r) below 1 km, the energy spectrum of very surface floaters shows much slower decay r^1/3, i.e., higher energy, compared to CARTHE drifters, which shows similar 2/3 power law turbulence scaling as CODE drifters (Poje et al., 2017). Overall the results highlight the complexity of the surface response and the very high shears involved. In all the tested cases drifters appear to provide consistent dynamical measurements, even though it is possible that in presence of very high winds and breaking waves spurious effects could occur.

For the CALYPSO experiments, the effects of local winds and surface waves (sea state) on the velocities measured by the CARTHE, CODE, and SVP drifters have been analyzed in Poulain et al. (under review)¹. The shear of horizontal currents between the surface drifters and the SVPs is typically downwind with a magnitude of 0.5% of the wind speed. However, maximum velocity differences between CODE/CARTHE and SVP of 15–20 cm/s can occur for all the sea state conditions encountered, even for winds of a few m/s. While CARTHE and CODE drifters have similar water-following capabilities, the difference with the currents measured by SVP is likely due to multiple factors, including the different sampled physics (surface wave motion, shear due to Stokes drift, Ekman currents and any other, mainly ageostrophic, currents in the top layer), direct windage and possible slip biases. Unraveling these aspects is very difficult, but it should be noted that for computation of kinematic properties, that are the focus of this paper and that depend on differential values, the scale of the wind and surface waves, and presumably of the top layer shear, are typically larger than the size of the triplets used for the computation (see section 2.5). Therefore, slip bias and shear problems are not expected to significantly affect the statistical estimates if the scales considered are smaller than 10 km or so. We anticipate (section 3) that statistical differences between surface and 15 m are found to decrease at scales greater than 5 km, further suggesting that bias effects are limited. Nevertheless, further studies to specifically address this point are necessary especially in presence of high winds and wave breaking condition.

We also remark that SVP and CARTHE drifters are prone to drogue loss. For SVPs, drogue loss is automatically detected and flagged (Grodsky et al., 2011; Lumpkin et al., 2013), and during the CALYPSO experiments there was no such event reported. For CARTHE drifters, following the drogue loss during LASER experiment in 2016 (well-documented in Haza et al., 2018), a number of design changes were implemented by replacing the rubber connection between the floater and the drogue by a metal chain, as well as by significantly strengthening the attachment areas by a redesign of the molding. After that, the resilience of the drifters to drogue loss greatly improved. One of the first and most obvious sign of drogue loss is the periodic interruption of the GPS signal, as the floaters start being submerged under the action of waves (as explained in Haza et al., 2018) and this was not seen in CALYPSO data.

In Table 1, top panel, the number of drifters deployed during the two CALYPSO experiments (in 2018 and 2019) is specified for each type of drifter. The datasets include drifter trajectories for the periods: May 28–June 29, 2018 and March 28–April 29, 2019. The strategy of deployments will be discussed in the following.

All drifters transmitted their GPS data (and other ancillary data) via either Iridium (SVP) or GlobalStar (CARTHE and CODE) satellites, every 5 or 10 min. Only the drifter GPS position data were considered for this work. The positions were edited for outliers (manually and using automatic procedures based on speed and acceleration criteria) and were linearly interpolated at 5 min intervals. Velocities were estimated by finite differencing the successive positions (central difference over 10 min). The positions and velocities were also low-pass filtered with a Hamming filter, producing two filtered data sets with cut-off at 30 min and 1 h, respectively. In the following analysis the quality-checked unfiltered data set is primarily used, while the filtered datasets have been considered to test the sensitivity of results to very short-time variability.

Here, we focus on drifter data taken during the two CALYPSO experiments in the periods May 28–June 29 in 2018 and March 28–April 29 in 2019, respectively. In order to provide a framework for the drifter results, we first characterize the two periods in terms of average ADT, mesoscale dynamics, and associated geostrophic velocity (Figures 1, 2) in the area covered by the drifters. We also show the typical stratification during the two periods (Figures 3, 4) as measured by UCTD sections during the CALYPSO cruises.

The Alboran Sea is known to be characterized by vigorous and semi-permanent mesoscale dynamics with anticyclonic eddies and associated strong density fronts (Gascard and Richez, 1985; Tintore et al., 1988; Benzohra and Millot, 1995). In particular, the Western Alboran Gyre (WAG) is typically present in the most western side of the basin, while the more variable Eastern Alboran Gyre (EAG) is located more to the East at the boundary with the Algerian basin.

During the 2018 experiment, centered on the month of June, the average mesoscale field shows a well-defined WAG and EAG (Figure 1), with a connecting boundary current and a smaller anticyclone between them along the African coast at approximately 2°-3°W. Another anticyclone is present further East, along the Algerian current around 3°-4°E. The meteorological conditions during the period of interest are characterized by mostly moderate winds (Figure 1), with average speed ≈ 6.5 m/s over the area of interest, reaching ≈ 10 m/s during three meteorological events. The wind direction is mostly zonal (Macias et al., 2016; Oguz et al., 2017). An example of the stratification during the experiment is shown along a section perpendicular to the western side of the EAG, in terms of Temperature (T), Salinity (S), and potential density (Figure 3). The results are typical also of the other sampled sections during the experiment (not shown). A surface mixed layer can be seen in the first ≈ 50 m, mostly driven by temperature that reaches ≈ 19°C at the surface. There are indications of shallower (around 20–25 m) frontal structures in the density within the mixed layer, that reflect horizontal gradients in both T and S, likely due to filaments of Atlantic and Mediterranean waters (see for example ≈ km 18 and 55 in the transect in Figure 3). We also notice the sharp pycnocline sloping between 100 and 200 m, that is associated to the anticyclonic mesoscale eddy.

During the 2019 experiment, centered on the month of April, the ADT average field (Figure 2), shows a well-defined WAG while the EAG is not developed. Several smaller anticyclones can be seen, one along the African coast around 3° W (as in the 2018 experiment), a weaker one more to the East off the African coast around 1° W, and one in the open sea around 3°W, 36° 30′ N. The wind during the 2019 experiment has ≈ 7 m/s mean speed, similar to the 2018 experiment, but is characterized by higher variability in direction and amplitude, with maximum values of almost 15 m/s. The stratification sampled during the 2019 experiment is different from the 2018 experiment, possibly because of interannual seasonal processes or because of local variability. A typical section across the open sea anticyclone is shown in (Figure 4). Compared to the 2018 section (Figure 3), there is no well-defined surface mixed layer in the upper 50 m, (even though the surface is starting to warm up reaching ≈ 16.5°C), and the main pycnocline tends to outcrop to the surface at the edge of the gyre. On the other hand, there are also some similarity with the 2018 surface stratification, in that density gradients are also present in the upper ocean (≈ 20–30 m) away from the outcrop region (see for instance around km 8 in Figure 4). As for June 2018, these surface fronts are likely to be associated to the presence of filaments of water of different Atlantic and Mediterranean origin (Lodise et al., 2020; Tarry et al., 2021). A similar picture, characterized by surface outcropping of the isopycnals around the edge of the mesoscale gyres and smaller surface density gradients occurring away from the outcrop, is commonly observed also in the other 2019 sections (not shown).

The deployment sites for the three types of drifters are shown separately in the panels of Figures 1, 2, color coded with time. The launching strategy involves both extended transects along targeted features, and cluster releases at different depths.

In this study we focus on the estimate of divergence and vorticity (KPs) as statistical trademark of submesoscale activity in the area of interest, since submesoscale processes are typically characterized by positive vorticity and convergence. A number of methods have been proposed in the literature to estimate KPs from drifter clusters. Most of them are based either on Least Square velocity estimation or on area rate of change estimates (Molinari and Kirwan, 1975; Okubo et al., 1976; Kawai, 1985), using clusters with a minimum configuration of three drifters.

Here we use the area rate of change method (Saucier, 1953; Molinari and Kirwan, 1975; Berta et al., 2016) and the smallest configuration, provided by drifter triplets, chosen because it is more flexible to resolve small scales. The computation is based on the fact that divergence can be expressed as the rate of change of the area of a water parcel, in our case represented by the triangular area delimited by the drifter triplet. The formulation for vorticity can be written in an analogous form, under appropriate rotations of the velocity vectors associated to the vertices of the triangle. Even though estimate precision increases with increasing number of drifters (Essink, 2019), the recent work by Berta et al. (2020b) has shown that statistical estimates from triplets are robust, based on the comparison with estimates from larger drifter clusters and from concurrent independent observations based on Eulerian X-band radar. These findings are also confirmed by modeling results by Sun et al. (2020).

In the following, we show results from estimates of horizontal divergence

and vertical vorticity component

where u, v are the x, y components of drifters' velocity, computed from triplets evolution during a time step Dt. Drifter triplets are re-sampled every 15 min (only triplets sharing at most one drifter are considered) and followed for the time interval Dt over which the KPs are estimated. Using independent chance triplets for a short time window, instead of following original triplets for longer time, allows to maximize the statistics and to minimize the effect of inhomogeneous sampling driven by convergent regions (Pearson et al., 2020; Sun et al., 2020). The KPs results presented in the following are expressed in terms of the local Coriolis parameter, f (8.572 × 10⁻⁵s⁻¹ at the reference latitude 36° N). KPs triplet-based estimates are valid as long as triplets maintain non-collinear configuration, while they deteriorate as triangles elongate. The aspect ratio of the triangles is evaluated at each time step using the deformation metric λ (Berta et al., 2016, 2020b)

where A and P are the area and perimeter of the triangle, respectively. The λ normalization ensures that λ=0 for collinear triangles and λ=1 for equilateral triangles.

For the purpose of investigating scale dependency of KPs, we separate triplets in selected classes of scales, depending on their size—i.e., the average side length (Berta et al., 2020b). The processes of class selection and KP estimation is performed in such a way to minimize possible error sources, as discussed in the following.

A main error source is the tendency toward deformation and alignment of the triplets since KPs estimation becomes unreliable when clusters tend toward collinearity (Berta et al., 2016; Ohlmann et al., 2017). In order to minimize this error, we perform a targeted triplet selection similarly to what was done in Berta et al. (2020b). At each time interval (15 min), we resample the data set choosing chance triplets of various class sizes that are quasi-equilateral, i.e., with values of λ greater than the cutoff value λ=0.7.

The other main potential error source is due to the GPS error on drifters' position, that is of the order of ≈ 10 m. In order to avoid error propagation from this source, we choose a time step Dt for the estimation of δ and ζ such that the increment in triplet size Dr during Dt is much larger than 10 m. Details on uncertainty estimation from GPS error are provided in the Supplementary Material, while the implementation of the procedures are provided in the next subsection.

As mentioned in section 2.1 the drifter data sets are available in three different formats depending on data treatment and filtering. The sensitivity of the KPs estimates to these formats and to the choice of Dt has been tested. As a first step, we have chosen a “benchmark configuration” that consists of the quality-checked unfiltered data set, from which KPs are computed using Dt = 30 min.

Triplets at different scales are selected from the data set as mentioned in section 2.4. The various classes of scales are chosen to maximize spatial resolution still retaining a statistically significant number of triplets, leading to a set of classes centered around values of 1, 3, 5, 9, and 16 km. Higher resolution classes, centered at less than 1 km, cannot be achieved because of the restricted number of triplets. The number of triplets in each class for each drifter type is shown in Table 1, bottom panel. KPs for each type of drifter (CARTHE, CODE, and SVP) are then computed for each size class and statistics are evaluated when the number of triplets is larger than a cutoff value of 85.

A few outliers have been removed (approximately 10 points for each type of drifter over thousands of samples), corresponding to isolated points outside the KP distributions. Typical values of Dr corresponding to Dt=30 min are evaluated as averages for each class, and in all cases they are found to be greater than 100 m, which is significantly bigger than the GPS error of 10 m (see also the Supplementary Material).

The robustness of the benchmark configuration is then tested using as metric the root mean square (rms) of the divergence δ (that is one of the main statistics considered in this study), and considering the sensitivity to the formats of the data sets and to different values of the time step Dt. The following tests have been considered:

Overall, the results show that the benchmark configuration is robust, and that minimizes the impact of errors on the estimates. We therefore confirm the use of the benchmark configuration, i.e, the unfiltered dataset and the time step Dt=30 min to compute KPs, and from now on we refer to these results only.

We start by providing a general description of the distribution of the trajectories, that are shown in Figure 5 color coded according to drifters' speed (the superimposed black dots indicate triplet distribution, that will be discussed in the next sections).

For both experiments, drifter trajectories show a good qualitative agreement with the main features of ADT in Figures 1, 2. During June 2018, the trajectories from the three drifter types sample the EAG and the eastern anticyclone in the Algerian current, while SVPs also partially sample the WAG. During April 2019, all drifter types sample the WAG and the African boundary current with the two smaller eastern anticyclones close to the coast. They also sample the small open ocean anticyclone at around 3°W, 36° 30' N, especially for the CARTHE drifters.

The energetics of the two experiments are similar, but stronger for 2018 probably because of the prevalence of the strong WAG and EAG gyres as shown also in the altimeter maps. Mean speeds are similar for the three types of drifters (CARTHE, CODE, and SVP) and have average values of ≈ 45 cm/s in 2018 and ≈ 40 cm/s in 2019.

In this section, we use the methodology described in section 2.4 and 2.5 to compute divergence and vorticity using triplets in the various classes of scales shown in Table 1, bottom panel. We notice that the number of triplets, especially for small scales, depends not only on the number of drifters (Table 1, top panel), but also on the details of the deployment (see Figures 1, 2). A typical spatial distribution of triplets is shown in Figure 5, considering as an example the positions (black dots) of triplets in the 5 km class characterized by higher divergence values, exceeding 0.5f. The 5 km class, corresponding to the range (4–6) km, has been chosen because is the smallest range with enough data to provide significant statistics for the three drifter types in both years. In the 2018 experiment, the triplets are mostly found along the front of the EAG (left panels) especially for SVPs. CARTHE and CODE triplets are also seen along the African boundary current reaching the eastern anticyclone in the Algerian Current. In the 2019 experiment (right panels) triplets are found in the WAG, in the African boundary current and in part of the open ocean smaller anticyclone around 3° W,36° 30' N.

As a first step in the analysis, we investigate the possibility of computing the statistics separately for the three drifter data sets. Specifically, the question that we address is whether or not there are enough data to significantly differentiate between the very surface drifters, i.e., CARTHE at 60 cm and CODE at 1 m. Previous works have shown the presence of significant shears in the first meter of water (Laxague et al., 2018; Lodise et al., 2019), suggesting that it would be interesting to keep them separate, but the number of data is limited, so that a preliminary investigation is needed.

We consider two classes of triplets as examples, one at the smallest resolved scale of 1 km (range <2 km), and the other at the larger 9 km scale (range 6–12 km), and we investigate the KP statistics, comparing the two surface drifters and the 15 m SVPs. Probability density functions (Pdfs) of divergence and vorticity are shown in Figures 6, 7, respectively, for the two classes and the 2 years. The corresponding statistics are shown in Table 2. As shown in Table 1, bottom panel, in both years for the smallest class (<2 km) only one type of surface drifters has sufficient data for analysis (i.e., more than 85 values) : CODE for 2018 and CARTHE for 2019.

For the smallest class (<2 km) the pdfs of divergence (Figure 6 upper panels) show a clear difference between the SVP and the surface drifters considered for each experiment. The SVP distribution is more peaked, while the surface drifters are more spread-out and characterized by higher standard deviation (std) and root mean square (rms) (Table 2, top panel), with individual values reaching 3f. For the vorticity (Figure 7 upper panels), the surface pdfs show prevalent positive values, with significant means (up to 0.8f for CODE 2018), and heavy positive tails especially for CARTHE 2019, with individual values exceeding 6f. The SVP pdfs have qualitatively similar properties but reduced values of std and rms (Table 2, bottom panel). For the larger class (6–12 km) the situation is different (Figures 6, 7 lower panels). For all drifter types and for both KPs, the pdfs are more peaked than for the smaller class. The std and rms values are smaller (Table 2), and there is no clear difference between the distribution of SVPs and the more surface drifters. The vorticity distribution is positively skewed for both classes and for both surface drifters, as it can be expected in presence of submesoscale processes in an area with strong lateral gradients and frontogenesis. For SVPs the skewness is positive for the smaller class, and slightly negative at the bigger scales. The divergence skewness is smaller than for vorticity and less well-defined, even though mostly negative.

Overall, the results indicate that the there is no detectable systematic difference between the statistics of the two surface drifter sets at the larger scales of 9 km, while at the smaller scales of 1 km the coverage is insufficient to differentiate among them. Differences with SVPs are comparable for the two surface sets. As a consequence in the following we analyze the two surface drifter data sets together. This allows increasing the number of triplets in each class (Table 1, bottom panel), giving a minimum of 123 surface triplets for the smallest class (<2 km) in 2018 and 290 in 2019.

The results presented above provide an interesting and articulated picture of the upper ocean between the surface and 15 m at horizontal scales of 1–16 km, bridging the meso and submesoscale range. Here we summarize the main findings and point out some interesting open issues.

At the surface, CARTHE and CODE drifters sampling the first meter of water indicate that the flow is highly ageostrophic at scales of 1 km, with typical rms values of order f for both KPs (higher than 1f for ζ and slightly less than 1f for δ) and individual values exceeding 3-6f. Positive vorticity skewness indicates the presence of submesoscale dynamics, in agreement with previous results (Shcherbina et al., 2013). Both surface KP values decrease with increasing horizontal scales, as expected, but there is a significant difference between divergence, that decreases faster with variance slope ≈ −4/3 at least up to scales 5–9 km, and vorticity variance that decreases more slowly, with a slope ≈ −4/5 or smaller, and tends to flatten for scales larger than 5 km, suggesting an interaction with mesoscales. The dynamical reasons behind these different decays will be further discussed in section 4, comparing our results with other literature results.

At 15 m, results from SVP drifters show that the flow has similar but attenuated characteristics with respect to the surface. The attenuation is especially evident for divergence, that has rms values less than 0.5f at 1 km and significantly slower decay with scales with respect to the surface, suggesting that the flow is less dominated by small scales. The vorticity shows similar but less marked results, with rms below 1f at 1 km and reduced scale dependence. The skewness values are positive at the smallest scales, but there is more variability with scales, with the occurrence of some negative values. This indicates that the flow is still ageostrophic but the influence of submesoscale processes is attenuated with respect to the surface.

From the dynamical point of view, the attenuation with depth points to the relevance of boundary layer processes (D'Asaro, 2014) that are expected to energize surface dynamics, increasing divergence and vorticity and influencing the smaller scales. This is indeed in agreement with TTW and TTTW results, that indicate that ageostrophic overturning is intensified in presence of strong vertical turbulence and lateral buoyancy gradients (Dauhajre et al., 2017). Also the presence of Langmuir cells, that can reach scales of 1 km (Malarkey and Thorpe, 2016; Suzuki et al., 2016) is expected to increase surface kinematic properties, as well the interaction between fronts and waves (McWilliams and Fox-Kemper, 2013; Suzuki et al., 2016). The wave interaction occurring through Stokes shear forces, can inject energy and momentum at the scales of the front, i.e., well beyond the actual scales of Stokes drift. Finally, ageostrophic overturning in shallow mixed layer fronts have also been shown to evolve in surface gravity currents (Pham and Sarkar, 2018), with strong across front velocities. All these results provide a useful general framework that goes in the direction of explaining the present results. There are some specific questions, though, that are still not resolved at least at the authors knowledge, and need further investigation. For instance it is not clear whether high wind stress directly correlates with enhanced surface values of KPs and with depth attenuation in the upper 15 m. Some mechanisms, such as the dependence of the overturning circulation on vertical mixing suggests that high wind stress induces higher surface KPs, but other mechanisms such as Ekman pumping (Thomas et al., 2013) and wave interactions (Suzuki et al., 2016) can act both ways, since they depend on the specific wind direction with respect to the front. Also other mechanisms related to atmosphere feedback to submesoscale fronts can act both ways (Renault et al., 2018). A statistical answer based on the present data set is not feasible, because there are not enough information on triplet sampling with respect to frontal structures nor enough co-located triplets. Further studies would be necessary to address these questions either considering specific events (Johnson et al., 2020a,b), or using model results (Oguz et al., 2017). Another interesting question arises from the difference in the attenuation of ζ and δ. A possible hypothesis is that higher values of ζ at larger scales and at depth are sustained by mesosocale interactions. Again, further studies are necessary to quantitatively address these points.

The results on scale and depth dependence are generally valid for both experiments, but it is interesting to verify whether there are some quantitative differences between them related to the different environmental characteristics, especially in terms of stratification.

Comparing the 2018 and 2019 results, differences can indeed be seen in the surface δ and ζ joint probability density (Figure 11). The flow in the 2019 experiment presents more extreme values of vorticity and divergence, suggestive of the occurrence of stronger submesoscale events, in agreement with the different stratification characteristics and possibly also with the higher wind effects. The bulk statistics, though, are similar for both experiments, and the 15 m attenuation of submesoscale signatures is seen in both cases, even though is slightly more evident in the 2018 experiment, as indicated by the more variable skewness.

A possible hypothesis to explain the limited statistical differences, despite the stratification differences shown in Figures 3, 4. is that submesoscale instabilities might occur mostly within the shallower filamentations with depth of ≈ 25–30 m, that are indicative of different water masses and that are present during both experiments (Lodise et al., 2020; Tarry et al., 2021). Even though an in depth investigation of the effects of stratification is out of the scope of this paper, we have quantitatively tested this hypothesis characterizing the stability properties (Thomas, 2008) of the two sections in Figures 3, 4 in the upper 50 m, computing the vertical and lateral buoyancy gradients, N² and M², and the balanced Richardson number RiB=f2N2/M4 (see Supplementary Figures 5, 6). Values of Ri_B smaller or close to 1 are typical of submesoscale regimes, and indicate regions prone to frontal instabilities such as symmetric instabilities (SI) (Taylor and Ferrari, 2010), forced symmetric instabilities (FSI) (Thomas et al., 2013), and ageostrophic baroclinic instabilities (ABI) (Boccaletti et al., 2007). The results in Supplementary Figures 5, 6 show values of N² and M² reaching 10⁻⁴ and 10⁻⁶ s⁻², respectively, for both experiments, in the range of previously observed values in submesoscale regimes (Thompson et al., 2016; Ramachandran et al., 2018; Johnson et al., 2020a,b). Also, while regions with Ri_B < 1 are restricted especially for the 2019 section, values close to 1 can be seen in both years. Actually, even though M² reaches the highest values in the outcrop region in 2019, values of Ri_B close to 1 are indeed more prominent in 2018, close to the surface filaments shown in Figure 3. These results therefore reinforce the fact that water masses gradients support submesoscale processes in both experiments.

Also, another factor to be considered is that the statistical description at scales of order of 1 km, as considered here, might not be sufficient to effectively diagnose the complete development of submesoscale frontal phenomena. Recent studies have shown that frontal structures are characterized by scales of the order of hundred to tens of meters and values of KPs of the order of tens of f (Gonçalves et al., 2019; Lodise et al., 2020). We will come back on this point in section 4 in the framework of literature comparison.

We start by comparing our statistical results on depth dependence over the first 15 m with results obtained by Tarry et al. (2021) using the same CALYPSO data set, but focusing on a specific event in June 2018, sampled also by a neutrally buoyant Lagrangian float (D'Asaro, 2003). Tarry et al. (2021) analyze a submesoscale instability at a surface front, where the mixed layer was about 15 m deep in the dense part, shallowing to less than 6 m at the front and then reaching 20 m in the lighter part. KPs estimates during the event show vorticity values reaching approximately 2f at the surface and barely 1f at 15 m. A strong convergence signal of approximately 2f is observed at the surface, while at 15 m divergence reaches 0.5f and does not show any obvious correlation with the surface. Indeed, measurements of vertical velocity from the float indicate that the divergence is not a linear function of depth between the surface and 15 m. A very high spatial variability in KP values is found at the surface, while at 15 m the values are smoother. Overall, the description of this specific event is in good agreement with our statistical findings, that indicate a decrease of small scale KP values at 15 m, especially evident for divergence, and increased variance at scales of the order of 10 km, indicating smoother divergent patterns.

Surface intensification of submesoscale KPs in the upper ocean has also been observed during a 1-day survey with a ship-towed Triaxus profiler in the California current (Johnson et al., 2020a). Vorticity and divergence were inferred at several depths in the first 20 m and in the layer 100–140 m, showing strong fluctuations in time and approaching values of order f near the surface.

Other experimental results in the literature regarding KPs variability do not resolve the upper ocean in the first 15 m. Statistical submesoscale depth dependence has been studied by Shcherbina et al. (2013), Buckingham et al. (2016) in the North Atlantic comparing currentmeter observations in the range of 50 m with deeper results at 350–500 m depth. In winter, when the mixed layer depth is deeper than 50 m, they report attenuation of KPs std of approximately a factor two and a marked difference in the shape of vorticity pdfs, with skewness values exceeding 1 at 50 m and decreasing to values close to zero at depth.

Our statistical results on the KPs dependence on horizontal scales are first compared with previous results by Berta et al. (2020b), obtained in a different part of the world ocean, i.e., in the surface Northern Gulf of Mexico (GoM), using a similar drifter analysis and covering a (partially) compatible range of scale. We notice, though, that the drifter sampling strategy in Berta et al. (2020b) is different from the one used here. While in the present analysis chance triplets over an extended geographical region and during 1 month period are used, providing an “average” picture of the dynamics in the Alboran Sea region, Berta et al. (2020b) consider short analysis periods of 1–3 days, providing a more synoptic description of specific environmental situations. The analysis is based on massive cluster deployments (order of 100–300 drifters) synoptically launched during the CARTHE experiments in 2012 and 2016 (D'Asaro et al., 2020), that specifically targeted submesoscale dynamics in a GoM area influenced by the Mississippi outflow. Four main releases are analyzed in Berta et al. (2020b), two during summer (S1, S2) and two during winter (W1, W2), characterized by different conditions. In particular, the winter W1 deployment targeted an energetic cyclonic eddy of approximately 20 km in diameter, characterized by strong surface density gradients associated with filaments of fresh Mississippi River water, while the W2 deployment occurred in an area with nearly homogeneous horizontal density during the onset of a storm with high winds.

In Figure 12, the scale dependence of the variance of δ and ζ are shown in a range between 0.2 and 16 km, covering both the Berta et al. (2020b) and our results. The scale dependence of divergence variance (Figure 12 upper panel) is quite consistent in all cases with an approximately −4/3 slope up to 5–9 km scales. The GoM summer values are the lowest, as it can be expected given the very shallow mixed layer in this area (5–10 m), while the two GoM winter cases are in the same range as the two Alboran Sea cases considered here. Overall, the consistency of results despite the differences in geographical and dynamical regions is noteworthy. In particular, the consistency with the W2 case, obtained in environmental conditions indicating weak submesoscale instabilities and strong wind effects (Berta et al., 2020b), suggests the importance of boundary layer dynamics in shaping divergence scale dependence.

The vorticity variance (Figure 12 lower panel) shows a more complex situation. Summer GoM and the winter W2 cases have lower vorticity variance and show a quick decay with scales, compatible with the same −4/3 slope as divergence. The GoM W1 and the two Alboran Sea cases, on the other hand, are characterized by higher variance values and slower decay (−4/5). This suggests an interaction with larger scale vorticity, and possible mesoscale transfer fueling vorticity values at intermediate scales (Lodise et al., 2020).

Scale dependence of KPs have also been considered by Ohlmann et al. (2017) using surface drifters in the coastal California currents, based on repeated small cluster deployments targeting coastal eddies and frontal features. Ohlmann et al. (2017) consider the scale dependency of mean square (ms), rather than variance, in a range covering approximately 400 m–3.5 km (their Figure 5). The comparison between our ms results and the Ohlmann et al. (2017) ones is therefore limited to 2 data points, corresponding to approximately 1 and 3 km scales. A qualitative comparison of ms values shows that δ and ζ values are in the same range, but the Ohlmann et al. (2017) values are consistently higher. The values are 1.2f and 2.5f for δ and ζ at 1 km (vs. 0.7f and 1.6f in the Alboran Sea), and 0.3f and 1f at 3 km (vs. 0.25f and 0.5f). The difference is at least partially due to the fact that the California Current deployments were targeted to eddies and fronts, that are characterized by high KP values. The ratio between 1 km and 3 km values are in the range of 2.5–3 in all cases, except for the divergence in Ohlmann et al. (2017), that has a 4 ratio, suggesting a faster decay at those scales. Indeed Figure 5 of Ohlmann et al. (2017) shows a very fast decrease of δ for scales greater than 1 km, while ζ oscillates and decreases slowly. The difference between the δ and ζ decay is consistent with our findings, but even more marked. We also mention that scale dependence has been considered by Johnson et al. (2020a) using Triaxus and SEASOR data at scales of 5 and 12 km during the frontal survey in the California current. Even though the methodology is quite different, they find consistent values (of order 0.7f at 5 km and 0.1f or less at 10 km) with more pronounced scale decay for δ.

The statistical properties presented here provide a description of KP properties at a given scale between 1 and 16 km. As discussed in section 3.3, though, the characteristics of specific frontal events might need a more detailed description, likely including smaller scales. Lodise (2020) has recently proposed that values of maximum vorticity and divergence could provide useful metrics to characterize frontal events, and has provided a literature overview (his Figures 3.20, 3.21) of such values for scales smaller than 1 km, considering various types of instrumentations and methodologies. Here we report the Lodise (2020) results in our Figure 13, integrated with results from the present study in the Alboran Sea at 1 km and from Berta et al. (2020b) in the GoM at the available scale of 500 m (their Figures 10, 11).

The results in Figure 13 have been taken in various world oceans, including the Gulf of Mexico, North Atlantic, California Current, and Mediterranean Sea (see legend). They are mostly obtained from surface drifters, aside from Shcherbina et al. (2013) that are taken from ADCP, Rascle et al. (2017) from visible camera on airborne surveys, and Lund et al. (2018) from shipboard X-band radar (see also Table 5 in Lodise et al., 2020). Also, various methodologies have been used to analyze drifter data, following two main approaches. The first approach, (used in Ohlmann et al., 2017; D'Asaro et al., 2018; Berta et al., 2020b) includes techniques that estimate KPs directly from drifter clusters, provided that the cluster aspect ratio is beyond a certain threshold value. The second approach (used in Gonçalves et al., 2019; Lodise et al., 2020) is based on the reconstruction of the Eulerian velocity field using a Gaussian Process Regression (GPR), from which KPs are then evaluated. While also this technique is limited by drifter availability and cluster deformation, the flow parameter estimation allows a greater flexibility reaching smaller scale estimates.

Figure 13 shows the relationship between maximum δ and ζ values and the horizontal scale of the estimates. For the GPR based studies (Gonçalves et al., 2019; Lodise et al., 2020), the horizontal scale corresponds to the smallest correlation scale evaluated for each analyzed event. The horizontal scales associated with the visible camera measurements made by Rascle et al. (2017) are defined by the two different estimates of the width of front sampled (30 and 50 m). For statistical analyses at relatively coarse resolution (Shcherbina et al., 2013; Ohlmann et al., 2017; Lund et al., 2018; Berta et al., 2020b and the present study), for each class of scale, the maximum values from the distributions of δ and ζ are chosen and plotted against the resolution of the other datasets. Assigning a horizontal scale to the maximum kinematic values found by D'Asaro et al. (2018) is less trivial due to the fact that the ellipses method used becomes less accurate as drifters converge and become less isotropically spaced. Lodise et al. (2020), as well as this study, choose a corresponding length scale of 300 m. Varying this scale between 200 and 500 m still results in good agreement with the trends seen in Figure 13.

Most of the values in Figure 13 are clustered in a range between 300 m and 1 km, where most of the observations are taken. In this range, the maximum divergence values are approximately three times larger than the corresponding variance values in Figure 12 (upper panel), with the summer values in the GoM being distinctively smaller than the others. Despite the fact that the results are obtained from various instruments and in various locations of the world ocean, the divergence trend in this range shows a well-defined quasi-linear decrease at increasing scales, qualitatively consistent also with what shown in the variance plot (Figure 12 upper panel). We notice that the values include also the Alboran Sea estimates by Lodise et al. (2020) with GPR, corresponding to δ approximately 7f at 550 m. A second cluster of values is found for smaller scales in the range 70-30 m, with maximum δ values reaching almost 100f at the smallest scales. The highest values are obtained by Rascle et al. (2017) using visible camera on airborne surveys over the GoM which sampled a very narrow front, 30–50 m in width, followed by GPR estimates at 70 and 30 m also in the GoM. The vorticity plot (Figure 13 lower panel) shows similar general features, even though in the range 300m–1km the values are more scattered, in keeping also with the variance plot (Figure 12 lower panel). The Lodise et al. (2020) estimate for the Alboran Sea reaches 11f at 550 m. The highest values of order 100f are confirmed at 30 m for the camera measurements, while GPR estimates at those scales reach 30f.

Overall the results indicate that there is a good consistency between different types of instruments and methodology at various scales. They also show that frontal features can reach small scales of the order of tens of meters with KPs values of the order of 10–100f, and that targeted techniques are necessary to reach such resolution.