Surface evolution and wind effects during a cyclonic eddy splitting event in the Balearic Sea

During the period of 23–28 February 2022, a cyclonic eddy in the Balearic Sea was observed to split into smaller eddies. Serendipitously, a wealth of data was collected of the event, including satellite chlorophyll maps, Lagrangian drifters at several depths, hydrographic sections intersecting the splitting eddy, and wind speed and direction. Sufficiently many drifters were in the area to estimate kinematic properties (divergence, vorticity, and strain rate) from clusters along the edge of the eddy before and during the elongation period that led to the splitting. The vertical velocity w can be computed from colocated divergence values from surface drifters (CARTHE and CODE, within the top meter) and near-surface drifters (SVP, at 15 m depth). Together, the observations delineate the process of eddy elongation, leading to vorticity and strain-rate intensification on 25 February, followed by the collapse of the ridge in the center of the eddy and the emergence of smaller eddies on 26 February, terminating with the splitting into submesoscale cyclones on 28 February. Dense drifter observations and daily hydrography supplement the remotely sensed descriptive view of the eddy splitting process. In particular, they confirm dominant internal dynamics, consistent with isopycnal doming, but also point to a role played by the winds, which shifted from predominantly southwesterly to predominantly northeasterly and strengthened significantly before weakening again in the area of interest during the eddy splitting period. Nonlinear Ekman pumping w_Eknl is estimated from the wind data and drifter-derived vorticities to capture the contribution of the wind effects to the patterns of up- and downwelling accompanying the eddy splitting. The w_Eknl patterns are consistent with the drifter-based w and divergence estimates. Moreover, the nonlinear Ekman pumping is found to be of the same order of magnitude as (though generally less than) w, suggesting that the wind likely influences the observed surface processes.

1 Introduction

Ocean eddies are ubiquitous and are subject to complex dynamics (Fu et al., 2015; Petersen et al., 2013; McWilliams, 1985). Much of what is known about them has been informed by remote sensing. Global satellite observations have permitted quantitative characterizations of eddies, at least at the mesoscale and the upper end of the submesoscale, compatible with the resolution of the observations (Martínez-Moreno et al., 2021; Fu et al., 2015). Both globally (Chelton et al., 2011) and in regional seas like the Mediterranean Sea (Escudier et al., 2016; Mkhinini et al., 2014), large cyclonic eddies have been found to be on average shorter-lived than anticyclones. This is consistent with theoretical results that cyclones are more prone to weakening and splitting, due to several physical mechanisms. Sufficiently large angular momentum can lead to splitting of isolated rotating vortices, but this condition is only met by intense cyclones (Nof, 1990, 1991). The impact of free surface dynamics has been studied by Arai and Yamagata (1994) and Polvani et al. (1994), confirming the asymmetry between cyclones and anticyclones and pointing to the appearance of a saddle point inside the cyclone core as a necessary splitting condition (Arai and Yamagata, 1994). In addition, strain deformation has been found to play a key role in inducing enhanced cyclone weakening (Graves et al., 2006).

Numerical model studies (Amores et al., 2018) have suggested that at the present spatial resolution of gridded altimetry products such as AVISO/CMEMS, which are frequently used as the basis for eddy counts and analyses, only approximately 16 % of eddies in the Mediterranean Sea can be detected and that the unresolved features can significantly alias the larger structures. Moreover, Stegner et al. (2021) have shown that cyclone detection in the Mediterranean is significantly less reliable than anticyclone detection, because large cyclones tend to split into small submesoscale cyclones characterized by rapid dynamics that cannot be resolved by the altimetry data. These smaller cyclones are, thus, likely to be generally undercounted and yet important to the ocean dynamics. Indeed, sunglitter reflections observed from space are consistent with their pervasiveness (Munk et al., 2000). The process of cyclone deformation and splitting typically involves various scales and could play an important role at the submesoscale, coinciding with the forward energy cascade which dissipates energy from larger to smaller scales in the ocean (McWilliams, 2019). Moreover, cyclones are associated with nutrient upwelling due to isopycnal doming, potentially triggering biological blooms (Cao et al., 2024). Cyclone splitting, on the other hand, is characterized by enhanced vertical exchanges along the water column, influencing the subduction of surface waters rich in carbon to the interior ocean, a fundamental process in the biological carbon pump (Resplandy et al., 2019).

To properly capture the distribution and dynamics of submesoscale eddies, satellite observations must be supplemented with in situ data. During a recent experiment in the Balearic Sea (Northwest Mediterranean Sea), as part of the CALYPSO project (Mahadevan et al., 2020), a comprehensive suite of observations was collected during a cyclonic eddy splitting event, including satellite imagery, ship-based Eulerian measurements of the hydrography, and Lagrangian data from drifter releases at multiple depths within the top 15 m. The eddy was originally identified in remotely sensed chlorophyll-A maps (see Animation S1 in the Supplement), then targeted for in situ measurements. Its initial radius was ≈10 km (Fig. 1), with a shallow mixed layer of ≈20 m (Middleton et al., 2025). Over the course of five days, it split into two smaller cyclones with radii less than ≈5 km, i.e., much smaller than the typical Rossby radius in the area of ≈10 km (Escudier et al., 2016).

Figure 1 Top: Maps of the Balearic Sea with chlorophyll concentration from satellites (background color), wind (vectors), and in situ observations, with underway CTD sections in red and CaC and SVP drifter trajectories in cyan and blue, respectively, for (from left to right) 23, 26, and 28 February 2022. Bottom: Overview of the Mediterranean Sea, indicating the position of the top maps (white box). The black box in all frames shows the position of the geographic panels in Fig. 5 and later.

Middleton et al. (2025) compared the subsurface characteristics of this eddy breakup, derived from a variational mapping of subsurface hydrography observations to an idealized model without external shear or wind stress and determined that the eddy shape, size, and intensity alone could have precipitated the splitting. Meteorological observations, however, confirm that the eddy splitting was accompanied by significant winds, which almost certainly impacted events. In particular, the wind can drive strong localized vertical transport through nonlinear Ekman pumping, enhancing subduction associated with unforced eddy splitting, as suggested by Middleton et al. (2025). Our aim is to quantify this impact. The plentiful drifter observations facilitate this analysis, allowing us to infer fine-scale divergence, vorticity, and strain rate, as well as localized vertical transport. Wind intensification effects due to non-linear Ekman pumping on divergence and vertical velocity in the sea surface layer can thus be captured. The focus here is on the evolution of the eddy within the surface layer and its interaction with the wind, as illuminated by surface (drogued at 0.4 and 0.75 m) and near-surface (drogued at 15 m) drifter observations.

An overview of the observed eddy splitting event and the physical processes likely associated with it is given in Sect.  2. Sections 3 and 4 detail the observations and analysis methods, respectively. Section 5 provides a detailed description of the eddy evolution based on drifter and hydrography data, followed by the analysis of wind effects and their contribution to the observed dynamics. The paper wraps up with a discussion and conclusions in Sect. 6.

2 Background

2.1 The near-surface dynamics of eddy splitting

The eddy splitting event examined here was observed northwest of the island of Mallorca in the Balearic Sea, in the Northwestern Mediterranean Sea (Fig. 1, bottom panel). The area is characterized by intense mesoscale and submesoscale eddy activity, fed by instabilities of the Northern Boundary Current and the Balearic Front (Pinot et al., 1994; Ruiz et al., 2009). The mixed layer depth at the time was relatively shallow, approximately 20 m, as deduced from the Brunt-Vaisala frequency calculated from in situ hydrography (Fig. S1; see also Middleton et al., 2025). The eddy could be identified in satellite chlorophyll concentration maps as early as 17 February 2022 (Animation S1) and was targeted for intense sampling starting 23 February, shortly before it began elongating (Fig. 1, left panel). The wind in the area was approximately from the southwest, roughly parallel to the eddy's major axis, with an intensification lasting about one day (Fig. 2). Subsurface hydrography for the area shows an elliptic cyclone with doming isopycnals, characterized by subsurface vorticity reaching values of the order of the inertial frequency f (Middleton et al., 2025), indicative of a Rossby number R of order 1. The smoothed divergence field, averaged over the top 60 m and 37 h, exhibits a quadrupole pattern (Middleton et al., 2025), consistent with cyclogeostrophic effects that in highly curved flows can alter the geostrophic balance. In essence, as the flow approaches the narrow end of an elliptical eddy, it accelerates due to centrifugal forces creating convergence, while it diverges after passing the tip.

Figure 2 Top: Wind maps at 10 m from ERA5 The red (cyan) path shows the R/V Pelagia (R/V Pourquois Pas?) ship track on the corresponding dates (23, 25, 26 February). The black square indicates the area over which the wind has been averaged to obtain the wind time series plotted in panel (D) (wind speed) and (E) (wind stress) together with in situ wind observations from the two oceanographic vessels involved in the CALYPSO 2022 experiment (R/V Pourquois Pas? in blue and R/V Pelagia in red).

Over the next few days, the elongation is observed to increase, as the eddy strains (Fig. 1 top, central panel. Refer also to the Animation S1). Idealized modeling, suggests that the convergence areas of the quadrupole trigger a positive feedback that intensifies positive vorticity in two distinct patches, leading to the splitting (Middleton et al., 2025). By 26 February, the center of the ridge collapses, while the smaller eddies intensify and become fully separated by 28 February (Fig. 1 top, right panel). Meanwhile, the wind, which was weak throughout 24 February, strengthens swiftly, likely due to a Tramontane event, starting on 25 February and lasting until 27 February, with prevalent direction from the northwest (Fig. 2).

Thus, the eddy experiences both internal instabilities triggered by the eddy shape and its high intensity (high Rossby number) and significant wind forcing, especially during the period 25–27 February. These processes interact nonlinearly. The vertical vorticity ζ in the eddy modifies the `effective' Coriolis parameter $f_{\text{eff}}=f+\mathit{\zeta }/\mathrm{2}$ and thereby influences Ekman transport, potentially leading to nonlinear Ekman pumping (Stern, 1965; Niiler, 1969; Hart, 2000; Wenegrat and Thomas, 2017; McGillicuddy et al., 2007; Mahadevan et al., 2008). This, in turn, can impact the pressure field and therefore the eddy structure. The vorticity gradients associated with the edges of eddies can also modify the local nonlinear inertial oscillations (NIO), leading to surface convergence and divergence in the surface layer with associated pumping effects (Perkins, 1976; Kunze, 1985; Elipot et al., 2010; Asselin and Young, 2020; Esposito et al., 2023). Moreover, Ekman and NIO effects can interact during periods of transient forcing, through energy transfer from high to low frequencies (Chen et al., 2021).

Ekman transport requires a period of at least two to three days under quasi-steady forcing to be fully established, although the time scale depends on several parameters including the strength of the wind forcing and eddy viscosity (Kirincich and Barth, 2009; Berta et al., 2018). On shorter time scales, only a fraction of the full Ekman transport is found in the surface layer. The event investigated here is characterized by high variability in the wind forcing (Fig. 2). The evolution and splitting of the eddy occurs over five days, while the wind variability has a scale of one to three days. This complicates the analysis. We anticipate that the interaction between eddy and wind effects is less straightforward than in the classical quasi-steady picture (Wenegrat and Thomas, 2017).

In order to identify the wind effects, we propose the following strategy. The densely distributed drifter data, from repeated launches in the first meter and at 15 m depth, are used to compute the kinematic properties (KPs), namely vorticity, strain rate, and divergence. Where drifters at the two depths are co-located, it is also possible to estimate the vertical velocity w (Esposito et al., 2023). This allows us to describe the surface layer evolution for 23–28 February, with a temporal resolution of 1 d. Wind effects are expected to be particularly evident in the divergence field, since it is directly influenced by wind (Chen et al., 2021), while vorticity and strain rate are mostly dominated by mesoscale and submesoscale dynamics (Esposito et al., 2021). While the high temporal variability in the wind inhibits the establishment of a fully developed Ekman transport and associated pumping, the tendency toward that state is expected to influence the surface layer. We therefore estimate the theoretical Ekman transport and pumping during the main wind episodes, following the theory briefly reviewed in Sect. 2.2, and relate these estimates to the observed divergence.

As mentioned above, NIOs can also impact the divergence field in the vicinity of the eddy. The daily averaging will dampen their signal but cannot fully eliminate it, in part because the inertial period (approximately 19 h) differs from the 24 h averaging time, but more fundamentally because of the nature of Lagrangian sampling, which is by definition not fixed in space. The quickly evolving eddy field necessarily makes the interpretation of the NIO effects on divergence more qualitative than is possible in the context of a quasi-steady jet (Esposito et al., 2023). We have to rely on the detection of smaller scale patterns and deviations with respect to the Ekman tendency.

2.2 Ekman transport and pumping

According to linear Ekman theory (Cushman-Roisin and Beckers, 2011), non-zero wind stress curl induces positive or negative Ekman pumping, referring to the upwelling and downwelling responses, respectively, of the Ekman layer. Assuming an eddy viscosity of 0.05 m² s⁻¹, a typical Ekman depth at the relevant latitudes is roughly 30 m. Recall that for weak surface ocean currents, the wind stress exerted on the ocean surface is given by

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\mathrm{wind}}={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{D}}{\mathit{u}}_{\mathrm{wind}}\left|{\mathit{u}}_{\mathrm{wind}}\right|\end{array}$$

where C_D is a drag coefficient, ρ_a is the air density, $\left|{\mathit{u}}_{\mathrm{wind}}\right|$ is the wind speed (at 10 m). Non-zero wind stress induces Ekman transport ${\mathbf{M}}_{\mathrm{wind}}=\langle U_{\mathrm{wind}},V_{\mathrm{wind}}\rangle$, integrated over the depth of the surface Ekman layer, that is directed perpendicular to the stress and rotated to the right in the northern hemisphere, with magnitude

$$\begin{array}{}\text{(2)}& \left|{\mathbf{M}}_{\mathrm{wind}}\right|={\displaystyle \frac{\left|{\mathit{\tau }}_{\mathrm{wind}}\right|}{{\mathit{\rho }}_{\mathrm{0}}f}}\end{array}$$

where ρ₀ is the mean water density in the Ekman layer. Non-zero divergence in the Ekman transport gives rise to Ekman pumping, with a vertical velocity proportional to the curl of the wind stress:

$$\begin{array}{}\text{(3)}& w_{\mathrm{wind}}=\mathrm{\nabla }\cdot {\mathbf{M}}_{\mathrm{wind}}={\displaystyle \frac{\mathrm{\nabla }\cdot ({\mathit{\tau }}_{\mathrm{wind}}\times \mathit{k})}{{\mathit{\rho }}_{\mathrm{0}}f}}={\displaystyle \frac{\mathit{k}\cdot \mathrm{\nabla }\times {\mathit{\tau }}_{\mathrm{wind}}}{{\mathit{\rho }}_{\mathrm{0}}f}}.\end{array}$$

Here k is the upward unit vector.

Nonlinear effects modify these formulas (Stern, 1965; Wenegrat and Thomas, 2017). When ocean currents $\mathit{u}=\langle u,v\rangle$ have significant magnitude relative to the wind, the stress τ acting on the ocean surface depends on the relative velocity ${\mathit{u}}_{\mathrm{rel}}={\mathit{u}}_{\mathrm{wind}}-\mathit{u}$:

$$\begin{array}{}\text{(4)}& \mathit{\tau }={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{D}}{\mathit{u}}_{\mathrm{rel}}\left|{\mathit{u}}_{\mathrm{rel}}\right|\end{array}$$

and the associated Ekman transport and pumping velocity are modified by the vorticity ζ in the Ekman layer:

$$\begin{array}{}\text{(5)}& {\displaystyle }\left|{\mathbf{M}}_{\mathrm{Eknl}}\right|={\displaystyle \frac{\left|\mathit{\tau }\right|}{{\mathit{\rho }}_{\mathrm{0}}(f+\mathit{\zeta })}},\text{(6)}& {\displaystyle }{w}_{\mathrm{Eknl}}=\mathrm{\nabla }\cdot {\displaystyle \frac{\mathit{\tau }\times \mathit{k}}{{\mathit{\rho }}_{\mathrm{0}}(f+\mathit{\zeta })}}.\end{array}$$

The total Ekman pumping velocity w_Eknl can be decomposed into two parts, one depending only on the vorticity and the other also involving the vorticity gradient (Song et al., 2020):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{w}_{\mathrm{Eknl}}& =w_{\mathit{\zeta }}+w_{\mathrm{\nabla }\mathit{\zeta }}={\displaystyle \frac{\mathit{k}\cdot \mathrm{\nabla }\times \mathit{\tau }}{{\mathit{\rho }}_{\mathrm{0}}(f+\mathit{\zeta })}}+{\displaystyle \frac{\mathit{k}\cdot (\mathit{\tau }\times \mathrm{\nabla }\mathit{\zeta })}{{\mathit{\rho }}_{\mathrm{0}}(f+\mathit{\zeta }{)}^{\mathrm{2}}}}\\ & ={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}(f+\mathit{\zeta })}}\left({\displaystyle \frac{\partial {\mathit{\tau }}^y}{\partial x}}-{\displaystyle \frac{\partial {\mathit{\tau }}^x}{\partial y}}\right)+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}(f+\mathit{\zeta }{)}^{\mathrm{2}}}}\\ & \left({\mathit{\tau }}^x{\displaystyle \frac{\partial \mathit{\zeta }}{\partial y}}-{\mathit{\tau }}^y{\displaystyle \frac{\partial \mathit{\zeta }}{\partial x}}\right),\end{array}\end{array}$$

where $\mathit{\tau }=\langle {\mathit{\tau }}^x,{\mathit{\tau }}^y\rangle$ and we have assumed an f-plane approximation. In order to concentrate on the contribution to w_Eknl due to ocean currents (Gaube et al., 2015), and in agreement with the observations, we assume a clear scale separation between the wind velocity, that provides a smooth large scale background, and the ocean currents, characterized by mesoscale and submesoscale eddies with significant local vorticity ζ, while neglecting the nonlinear dependence of the stress on the relative velocity.

The two terms in w_Eknl describe two different mechanisms. The first one, w_ζ, (Dewar and Flierl, 1987; McGillicuddy et al., 2007) results from the curl of the relative surface stress τ that, according to the assumptions stated above, is dominated by the wind velocity, scaled by the `effective' Coriolis parameter f+ζ, which is modified by the local surface vorticity ζ. Within a mesoscale eddy, the resulting Ekman pumping is expected to have opposite sign with respect to the local vorticity, resulting in upwelling in an anticyclonic eddy and downwelling in a cyclonic one. Gaube et al. (2015) evaluated w_ζ for an isolated idealized cyclone, finding downwelling in the core, with elongation in the direction of the wind.

The second term, w_∇ζ, (McGillicuddy et al., 2007) results from the vorticity gradients, divided by the square of the “effective” Coriolis parameter. Conceptually, within our assumptions, this term is due to the fact that Ekman transport M_Eknl (Eq. 5) decreases with increasing magnitude of positive ζ, inducing downwelling (and upwelling for decreasing positive ζ), while the opposite is true for negative ζ. This mechanism leads to the formation of dipoles within the core of an isolated eddy (Gaube et al., 2015), as shown in the cartoon in Fig. 3. For a cyclone, downwelling is expected to occur on the left side of the eddy core with respect to the wind (in the northern hemisphere) and upwelling on the right side. This is because the transport M_Eknl is directed to the right of the wind, and therefore enters the eddy from its left side with respect to the wind, leading to downwelling in order to adjust for increasing ζ and decreasing transport. The opposite occurs on the other side of the eddy, where ζ decreases and transport increases. For anticyclones the dipole structure is inverted with respect to the wind. Surface wind forcing has also been found to induce a vertical tilt in the cyclone (Stern, 1965). Eddy tilt is also associated with stronger vertical velocities (Li et al., 2022). Unfortunately, our dataset is insufficient to study such a tilt conclusively for the present eddy.

Figure 3 Schematic of wind-induced 3D transport on a marine cyclone.

3 Observations

3.1 Drifter Data

The analysis here is based on the surface and near-surface observations of drifter trajectories, complemented by ancillary data from satellites and ship-board instruments. Several significant drifter deployments were carried out during the CALYPSO cruise in February/March 2022 to sample the near-surface currents in specific submesoscale features with high horizontal resolution (order of 1 km). Three types of drifters were released in or near the cyclonic eddy considered in this paper: CARTHE (Consortium for Advanced Research on the Transport of Hydrocarbon in the Environment) drifters, drogued at 0.4 m (Novelli et al., 2017; D'Asaro et al., 2018), CODE (Coastal Ocean Dynamics Experiment) drifters, drogued at 0.75 m (Davis, 1985; Poulain, 1999), and SVP (Surface Velocity Program) drifters, drogued at 15 m (Niiler, 2001; Lumpkin and Pazos, 2007; Centurioni, 2018). CARTHE and CODE drifters have been found to measure currents representative of the top meter of the water column with accuracy of a few cm s⁻¹, assuming limited wind or wave induced slippage and Stokes drift (Poulain et al., 2022), and will be treated collectively as “surface” drifters (following, e.g., Esposito et al., 2021; Tarry et al., 2021) with a reference depth of 0 m, referred to below as CaC drifters. The SVP drifters have been shown to follow the sub-surface currents within 1 cm s⁻¹ in winds up to 10 m s⁻¹ (Niiler et al., 1995). GPS positions were sampled every 10 min via GlobalStar satellites for CaC drifters and every 5 min via satellite Iridium telemetry for SVP drifters. The drifter positions were processed with standard methods for quality control (Menna et al., 2017) and interpolated to uniform 5 min intervals starting at the top of the hour. Velocities were estimated by central finite differencing of these positions. Drifter data in the area of interest on three days during the analysis is shown as cyan (CaC) and blue (SVP) curves in Fig. 1. The Supplemental Materials contain an hourly animation (S2) for the entire 6 da period. Of the 140 CARTHE drifters released during the cruise, 53 were in the area of interest (black box in Fig. 1) during this time. Similarly, 32 of 70 released CODE drifters and 22 of 68 deployed SVP drifters entered the area and time period of the analysis.

3.2 Wind Data

To evaluate the wind effects on the eddy, we use the ERA5 reanalysis 10 m winds from the European Centre for Medium-Range Weather Forecasts distributed by the Copernicus Climate Change Service (Hersbach et al., 2023). The data are available hourly on a 0.25°×0.25° grid. Here we use daily averages derived from the hourly archive; see vector field in Fig. 1.

For local validation of the reanalysis winds, we compare wind speed and wind stress to in situ observations obtained from the standard shipboard instruments during the CALYPSO 2022 cruise. Two ships participated in the experiment, the R/V Pourquois Pas? and the R/V Pelagia. Both were equipped with anemometers recording data at a nominal frequency of 1 Hz. These were mounted at 30 and 27 m, respectively. Resulting wind measurements were converted to 10 m winds using the power law from Johnson (1999). Wind stress magnitude is computed from the measured wind speeds using the bulk aerodynamic formula assuming negligible surface ocean currents compared to the wind $\left|\mathit{\tau }\right|=C_{\mathrm{D}}{\mathit{\rho }}_{\mathrm{a}}|{\mathit{u}}_{\mathrm{wind}}{|}^{\mathrm{2}}$, where C_D=0.0015 is the drag coefficient, ρ_a=1.2 kg m⁻³ is an average air density, and u_wind is the wind velocity at 10 m. Ship tracks and time series of wind speed and stress magnitude are shown in Fig. 2.

3.3 Ancillary Data

Context is provided by ancillary data collected during the cruise. The eddy was first identified for intense sampling in satellite ocean color imagery, as illustrated in Fig. 1. The phytoplankton chlorophyll concentrations shown are a product of the Copernicus Marine Environment Monitoring Service (CMEMS) that combines data from multiple satellites (SeaWIFS, MODIS, MERIS, VIIRS-SNPP, VIIRS-JPSS1, OLCI-S3A, and OLCI-S3B) (CMEMS, 2022; Volpe et al., 2019). Temporal resolution is daily, and spatial resolution is 1 km.

Hydrographic information is obtained from EcoCTD profiles taken along multiple sections during the experiment (Mahadevan and D'Asaro, 2025). The EcoCTD includes sensors for conductivity, temperature, and density (RBR Concerto CTD), a chlorophyll fluorometer and backscatter meter (WETLabs BB2F), and an oxygen sensor (JFE Advantech Rinko), all mounted in a weighted aluminum housing and tow-yoed manually from the stern of the R/V Pourquois Pas? with an Oceanscience UCTD winch (Dever et al., 2020). Typical operation of the EcoCTD was to make one profile to 250 m every 5–6 min at a vessel speed of 6 knots. Its data was calibrated against the CTD rosette. While multiple sections were sampled that intersected with the eddy, here we focus on the single temperature and salinity section for each day that is most nearly colocated with the drifter observations in time and space, shown as red lines in Fig. 1.

4 Method

4.1 Kinematic Properties Estimates from Drifters

A variety of methods exist for extracting KPs from drifter observations (Molinari and Kirwan Jr., 1975; Okubo et al., 1976; Kawai, 1985; Berta et al., 2015; Gonçalves et al., 2019). Here we obtain area-averaged KP estimates from drifter triplets, following Method I from Molinari and Kirwan Jr. (1975). Divergence δ, averaged over time and over the area A of a triangle defined by a drifter triplet, is computed from the change in that area over 15-minute intervals:

$$\begin{array}{}\text{(8)}& \mathit{\delta }={\displaystyle \frac{\log \left(A\right(t+\mathrm{\Delta }t\left)\right)-\log \left(A\right(t\left)\right)}{\mathrm{\Delta }t}}.\end{array}$$

Vorticity ζ, normal strain rate N, and shear strain rate S are computed similarly, from rotated velocity components, following Saucier (1955). The latter two quantities are combined in a total strain rate $\mathit{\alpha }=\sqrt{N^{\mathrm{2}}+S^{\mathrm{2}}}$. To eliminate those estimates that are subject to large errors due to the geometry of the drifter triplets, we filter out groups that do not satisfy the threshold Λ>0.2 on the scaled triangle aspect ratio $\mathrm{\Lambda }=\mathrm{12}\sqrt{\mathrm{3}}A/P^{\mathrm{2}}$, where A is the area and P the perimeter of the triangle (Huntley et al., 2022). Area-averaged KPs are known to be scale-dependent (Berta et al., 2020). So, only triplets with the distance between each pair of drifters (side length of the resulting triangle) between 1 and 10 km are used. To ensure independent sampling, triplets sharing more than one drifter with an already listed triplet are eliminated. The resulting estimates are associated with the corresponding triplet centers of mass, then binned in 2 km × 2 km boxes and averaged over bins and time for daily gridded maps.

These KP estimates are subject to errors due GPS position uncertainty, as well as due to representation errors arising since the advected triangles are generally not the triangles defined by their advected vertices and only partially addressed by the geometric subsampling (Huntley et al., 2022; Aravind et al., 2024). These errors are somewhat mitigated by the use of spatio-temporal means over each 2 km × 2 km box and each day.

GPS positioning uncertainties (order of 10 m) result in larger KP errors for small-scale triplets and/or small-time intervals due to their larger magnitude relative to the actual triangle area and/or area change. For triads of a scale of 1–16 km, these errors were found to impact 4 %–14 % of samples from an experiment in the Alboran Sea by Esposito et al. (2021), resulting in an error estimated to be <10 %. The dynamics here are expected to be comparable. However, the pre-processing of the drifter data (see Sect. 3.1) is designed to mitigate position error, so that the anticipated errors are actually smaller than these theoretical results.

While choosing larger triplet scales minimizes the impact of position errors (Berta et al., 2016; Esposito et al., 2021), the representation error grows with triplet size (Aravind et al., 2024). The latter is of particular concern for triangles that are far from equilateral (Huntley et al., 2022; Ohlmann et al., 2017). By applying the threshold of Λ>0.2, we remove those estimates most likely to be effected by these errors. Model results from Huntley et al. (2022) showed that with this threshold, more than 99 % of samples had an error less than 0.25f. Therefore, 0.25f could be considered as an upper bound for any individual estimate of δ. Since we are averaging both spatially (2 × 2 km bins) and temporally (24 h), the error associated with our results should be even less and a function of the number of samples in each spatio-temporal bin. Typically, 10 or more triplets were used for each individual bin in our results (see maps in Fig. S2), suggesting an uncertainty on our KP estimates of no more than 0.1f.

Drifter clusters can also be used to derive pointwise KP estimates by using least-squares fits (Molinari and Kirwan Jr., 1975; Okubo et al., 1976). Pointwise-estimates tend to have larger extremes than averages. So, a direct one-to-one comparison with area-averaged estimates is not advisable. However, as a form of validation, we carried out a qualitative comparison between the triplet-based area-averaged and cluster-based pointwise estimates of vorticity from CaC drifters with spatial scales of 2–4 km. (This is a subset of the estimates entering the bins used in the rest of the paper). For the least-squares fits, clusters were formed consisting of 6 or more drifters, all at most 4 km apart from the centroid, with an aspect ratio equal or larger than 0.1, defined as the ratio of the smaller to larger eigenvalue of the position covariance matrix (Essink et al., 2022). To ensure independent sampling, redundant clusters, i.e. those with overlapping drifter membership of more than 75 %, were removed. Figure 4 shows the vorticity scaled by the Coriolis parameter $f=\mathrm{9.4717}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹ derived from triplets (triangles) and from clusters (circles) for 23 February 2022. The results demonstrate consistency between the two very different methods, in line with method comparisons reported elsewhere (e.g. Berta et al., 2020).

Figure 4 Comparison of normalized vorticity estimates from drifter triplets and clusters (2–4 km scale) as derived from CaC drifters, for 23 February 2022. Top: Individual estimates from drifter triplets (triangles) and clusters (dots). Bottom: Bin averages obtained from estimates with each method, left – triplets, center – clusters. Bottom right: Scatter plot relating the bin averages with least-squares fit (gray, dashed line) and one-to-one line (black, solid line).

4.2 Vertical Velocity Estimates from Drifters

When horizontally colocated divergence estimates are available at different depths, vertical velocities can be approximated from the continuity equation $u_x+v_y+w_z=\mathrm{0}$, where u, v, and w are the eastward, northward, and upward velocity components, respectively. It follows that

$$\begin{array}{}\text{(9)}& w(-h)=w\left(h_{\mathrm{ref}}\right)+\underset{-h}{\overset{h_{\mathrm{ref}}}{\int }}\left(u_x+v_y\right)\mathrm{d}z.\end{array}$$

To obtain a vertical velocity estimate, both the vertical velocity at a reference depth h_ref and a vertical profile of horizontal divergence are needed. The reference depth is typically either taken to be the surface, where vertical velocities approximately vanish (Tarry et al., 2022), or a deep “depth of no motion” (Lodise et al., 2020). Here, we follow Tarry et al. (2022) and use the surface at z=0. In lieu of a continuous vertical profile of the horizontal divergence, we rely on a first-order approximation, interpolating linearly between the two known depths (Esposito et al., 2023). Thus, if we denote the known horizontal divergences at depths −h₀ and −h₁ as δ₀ and δ₁, then

$$\begin{array}{}\text{(10)}& w(-h)={\displaystyle \frac{\mathrm{2}({\mathit{\delta }}_{\mathrm{0}}{h}_{\mathrm{1}}-{\mathit{\delta }}_{\mathrm{1}}{h}_{\mathrm{0}})h-({\mathit{\delta }}_{\mathrm{0}}-{\mathit{\delta }}_{\mathrm{1}})h^{\mathrm{2}}}{\mathrm{2}(h_{\mathrm{1}}-h_{\mathrm{0}})}}.\end{array}$$

For our sampling depths of h₀=0 and h₁=15 m, the vertical velocity at h=15 m is then given by

$$\begin{array}{}\text{(11)}& w(-\mathrm{15})=\mathrm{15}{\displaystyle \frac{({\mathit{\delta }}_{\mathrm{0}}+{\mathit{\delta }}_{\mathrm{1}})}{\mathrm{2}}}.\end{array}$$

The uncertainty in the estimates for δ₀ and δ₁, in terms of median standard errors, are 0.08f and 0.07f, respectively. These values are typical for most bins, although for some of the bins near the edges of the computational domain, the standard errors reach 0.48f and 0.34f for δ₀ and δ₁, respectively (see Fig. S2). Propagating these standard errors to the w estimates yields a median standard error for w of 8 m d⁻¹.

5 Results

5.1 Daily surface eddy evolution

In this section, we describe the evolution of the eddy day-by-day during the splitting period, 23–28 February, skipping 24 and 27 February when insufficient data are available for a meaningful analysis. KP values are superimposed on satellite chlorophyll isolines, when available, for context and are normalized by the local Coriolis parameter f, taken to be $\mathrm{9.57}\times {\mathrm{10}}^{-\mathrm{5}}$, computed for a latitude of 40.5° N. Subsurface drifters are more sparse and primarily add insight to the eddy splitting analysis on 23 February, when the co-location of the two drifter datasets allows direct computation of the vertical velocity w. While available, we therefore do not show SVP results for other dates. The selected ecoCTD sections provide a three-dimensional context.

5.1.1 23 February

During 23 February (Fig. 5), the surface CaC drifters sample mostly the southern flank of the cyclonic eddy, characterized by northeastward flow, and the northern edge of an adjacent anticyclonic structure, indicated by the curvature in the trajectories (Fig. 5A). Velocities along the edge of the cyclone reach a maximum value of 0.94 m s⁻¹. The strong velocities sampled by the drifters straddle the boundary between the cyclone and its adjacent anticyclonic region, with vorticity ζ (Fig. 5B) showing high positive values, exceeding ζ=f, to the north, characterizing the cyclone core, and negative, less extreme values to the south. This is consistent with the UCTD section (Fig. 5C), showing the characteristic isopycnal doming within the cyclone, coincident with positive drifter-derived ζ (Fig. 5F) and a pronounced front separating the region with negative ζ. The strain rate (Fig. 5E) is especially pronounced in the northern part of the coverage, corresponding to the area of high positive ζ within the cyclone. As anticipated, divergence δ (Fig. 5D) shows a more complex pattern, that at first glance is not straightforwardly consistent with the eddy structure. It shows a prominent positive pattern, with values reaching f, roughly located at the boundary between positive and negative ζ, indicating upwelling, sandwiched between two smaller negative bands, suggesting downwelling. Along the section (Fig. 5F), δ is negative within the eddy and slightly positive at the front.

Figure 5 Surface and water column properties from in situ samplings in 23 February 2022: (A) surface currents derived from CaC (1m depth) drifter trajectories. (B, D, E) surface vorticity, divergence, and strain rate (normalized by local f) as derived from CaC drifter triplets; black thin contours show chlorophyll isolines (see Fig. 1). The black line indicates the location of the underway CTD potential density transect (kg m⁻³) shown in panel (C). Panel (F) shows the vorticity (upper panel) and divergence (lower panel) along the CTD transect.

The 15 m SVP drifters (Fig. 6) show similar patterns (Fig. 5), although the KP values are slightly attenuated, as can be expected (Esposito et al., 2021). Vertical velocities w, computed for the co-located bins (Fig. 6E), are mostly positive (upwelling), reaching values 𝒪(50–60 m d⁻¹), except for a small negative band to the north. The pattern mirrors the surface δ (Fig. 5D) in the overlapping region, confirming that surface divergence is a good indicator of upwelling and downwelling.

Figure 6 Sub-surface properties in 23 February 2022: (A) 15 m-depth currents derived from SVP drifter trajectories. (B, C, D) surface vorticity, divergence, and strain rate (normalized by local f) as derived from SVP drifter triplets; black thin contours show chlorophyll isolines (see Fig. 1). (E) Vertical velocity within the first 15 m depth layer, estimated through continuity equation from divergence at 1 and 15 m depth.

The presence of the anticyclonic structure to the southeast of the eddy and the high strain rate at the border of the cyclone suggest that large-scale straining could have played a role at least in the initial phase of the eddy elongation. The day of 23 February was also characterized by increased winds, that could have played a role in the observed upwelling and divergence patterns. This point will be discussed in Sect. 5.2.

5.1.2 25 February

For 25 February (Fig. 7), no satellite chlorophyll data are available to place the in situ observations into a larger context. However, the surface trajectories appear to sample the core of the cyclone, showing a recirculating elongated structure and reaching maximum velocity values of 0.75 m s⁻¹ (Fig. 7A). The drifters capture the strong positive vorticity ζ (Fig. 7B), exceeding 2f, in the center, coincident with the maximum strain rate (Fig. 7E), which also reaches values of 2f. The cyclone strengthening is accompanied by the strong isopycnal outcropping in the UCTD section (Fig. 7C), whereby the highest ζ values occur at the edge of the front, as expected (Fig. 7F). The values of δ along the section (Fig. 7F) are consistent with upwelling within the cyclone and downwelling at the southern front. The overall divergence pattern (Fig. 7D) shows mostly positive values northwest of the central core and negative values to the southeast, but exhibits additional complexity, with high small-scale variability.

Figure 7 Same as Fig. 5 but for 25 February 2022. Black thin contours show drifter trajectories (see Fig. 1). Chlorophyll concentrations are not available for this date.

5.1.3 26 February

During 26 February (Fig. 8), the surface drifters are mostly following a southwestward current, aligned with the wind, with maximum velocities of 0.79 m s⁻¹ (Fig. 8A). They do not show recirculation at the scale of the main eddy, but do capture recirculation in two smaller submesoscale eddies that have formed at the northern and southern extremes of their sampling region. The collapse of the original eddy structure can also be detected in the isopycnals in Fig. 8C, that have drastically deepened and flattened relative to the previous day. The vorticity (Fig. 8B) within the filament connecting the two smaller eddies is positive and strong, exceeding ζ=2f, with very high strain rate (Fig. 8E) in the same area. Positive vorticity maxima also occur within the two spun off eddies. The divergence (Fig. 8D) is mostly negative within the main filament and consistent with the observed downwelling associated with the collapsing of the isopycnals (Fig. 7C vs. Fig. 8C), with a drop of about 40 m in one day. Within the eddies the δ structure is more complex with alternating positive and negative values.

Figure 8 Same as Fig. 5 but for 26 February 2022. (Please note the reduced bin coverage around the ecoCTD transect in panel B compared to panel D) which explains the different length of the lines in panel (F).

5.1.4 28 February

By 28 February (Fig. 9), the original eddy structure has completely split into smaller cyclones. There are now three submesoscale eddies, with a tenuous connection between them, traced out by the surface drifters, which reach maximum velocity values of 0.70 m s⁻¹ (Fig. 9A). All the eddies are characterized by intense cyclonic vorticity reaching ζ=2f (Fig. 9B), with patches of high strain rate (Fig. 9E). Divergence (Fig. 9D) again shows high small-scale variability, although it is mostly positive, indicative of upwelling, inside the southernmost eddy. Upwelling is also consistent with the newly developing doming of the isopycnals seen in the section through this eddy (Fig. 9C). High positive vorticity coincides with the doming (Fig. 9F).

Figure 9 Same as Fig. 5 but for 28 February 2022.

5.1.5 Synopsis

The drifters illustrate the process of eddy elongation and permit the quantification of the strong intensification of vorticity and internal strain rate occurring on 25 February, followed by the collapse of the eddy center and the emergence of smaller eddies on 26 February, and terminating with the splitting into three submesoscale cyclones on 28 February. The overall picture of vorticity and strain rate evolution is consistent with the cyclo-geostrophic explanation proposed by Middleton et al. (2025) based on the analysis of subsurface data and idealized modeling , even though the values are significantly higher than reported there, likely due to a combination of surface intensification and less smoothing than in the previous study. This suggests that the process at the surface is also mostly dominated by internal dynamics, even though the drifters also provide evidence of a possible contribution of large-scale strain rate in the initial eddy elongation. The divergence patterns, on the other hand, are different, with the drifter data again showing significantly higher values but also more complicated structures than the smoothed quadrupole detected in the subsurface variational mapping of Middleton et al. (2025). The surface divergence field is characterized by bands of alternating divergence and convergence indicative of upwelling and downwelling at fine scales and in some cases opposite to what would be expected for cyclo-geostrophic effects. Since the wind was particularly strong on 26 February, at the height of the eddy splitting, we investigate its impact on the sea surface layer variability not resolved by subsurface hydrographic casts. These wind contributions to the observed divergence patterns and related vertical velocities are explored next.

5.2 Wind effects

Our investigation into possible wind effects considers contributions from nonlinear Ekman pumping w_Eknl (see Sect. 2.2) and NIOs (see Sect. 2.1), concentrating on the days of significant wind forcing, i.e., 23 and 25–27 February. As discussed in Sect. 2.1, even though the time scales are such that w_Eknl is not expected to be fully established, the tendency toward w_Eknl is expected to leave a signature in the observed patterns.

Recall that there are two components of nonlinear Ekman pumping (Eq. 7): w_ζ and w_∇ζ. For an isolated cyclone (Gaube et al., 2015), w_ζ is expected to lead to downwelling in the core, with elongation in the direction of the wind. Here, the observations (Sect. 5.1) clearly show the presence of a cyclonic eddy but are too sparse to delineate its complete spatial structure, let alone its surroundings. Consequently, we can only provide a qualitative description of the expected contribution of w_ζ. The observed pattern of isopycnal downwelling and increasing ellipticity observed during the eddy evolution (Sect. 5.1) is consistent with a tendency toward downwelling and eddy elongation along the wind direction. Its relevance, though, depends on the value of the expected w_ζ. An order of magnitude estimate, using $\left|{\mathit{u}}_{\mathrm{wind}}\right|\approx$ 7–10 m s⁻¹, $\left|\mathit{u}\right|\approx$ 0.10–0.5 m s⁻¹, $\left|\mathit{\zeta }\right|\approx$ 1 – 2f yields a corresponding range for w_ζ of ≈ 1–1.5 m d⁻¹, which is more than an order of magnitude smaller than the w derived from the drifter divergence calculations (Fig. 6) and implied by the observed isopycnal drop of more than 40 m on 26 February (Fig. 8). We therefore conclude that, even though w_ζ may contribute to the eddy evolution and its vertical velocity, its contribution is secondary and cannot explain the observed evolution.

The second term w_∇ζ depends on the vorticity gradient and on wind direction. As for w_ζ, we cannot provide a complete direct assessment of w_∇ζ since our measurements are spatially limited. Some considerations and estimates, though, can be provided combining the information from the wind data during the two events and the observed vorticity ζ. As a first step, we again perform an order of magnitude analysis, assuming the same ranges for $\left|{\mathit{u}}_{\mathrm{wind}}\right|\approx$ 7–10 and $\left|\mathit{u}\right|\approx$ 0.10–0.5 m s⁻¹, and for the vorticity gradients $\mathrm{\nabla }\mathit{\zeta }\approx \mathrm{\Delta }\mathit{\zeta }/\mathrm{\Delta }r\approx f/$(5–10) km, where r is the relevant spatial scale. The corresponding w_∇ζ is estimated in the range 5–20 m d⁻¹, which is comparable to the order of magnitude of the measured w.

A more precise estimate of w_∇ζ is possible where there is enough data coverage. The linear Ekman transport M_wind (Eq. 2) is shown in Fig. 10 (top panels) for 23, 25, and 26 February, respectively, with the nonlinear transport $\left|{\mathbf{M}}_{\mathrm{Eknl}}\right|$ (Eq. 5) superimposed wherever vorticity ζ is available from drifters (Figs. 5, 7, 8). As expected, $\left|{\mathbf{M}}_{\mathrm{Eknl}}\right|$ is lower than the background in the cyclonic eddy (positive ζ), but higher in the anticyclonic areas (negative ζ).

Figure 10 Top: Nonlinear Ekman transport magnitude and direction for (A) 23 February, (B) 25 February, and (C) 26 February as perturbed by surface vorticity estimated from drifter triplets coverage, superposed on the linear Ekman transport outside the drifter coverage. Bottom: Ekman pumping velocity derived from Ekman transport divergence for (D) 23 February, (E) 25 February, and (F) 26 February. In Panel (E) the thicker black line indicates the maximum vorticity magnitude.

A corresponding estimate of w_∇ζ for 23 February (Fig. 10D) shows a prominent upwelling pattern, reaching values of w_∇ζ≈ 20 m d⁻¹ in the transition area between positive and negative ζ (Fig. 5B). In this region, the transport $\left|{\mathbf{M}}_{\mathrm{Eknl}}\right|$ increases significantly in the direction of propagation of transport (Fig. 10A). Conceptually, this corresponds to the upwelling on the side of the cyclone to the right of the wind (Fig. 3). Some downwelling can be seen at the southern edge of the available data, where negative vorticity starts to decrease. The vertical velocity w computed from the drifters for the same day (Fig. 6E) shows a similar predominantly upwelling pattern in the same area. The surface divergence (Fig. 5D) is also consistently positive there, indicating upwelling.

A quantitative comparison between the estimated w_∇ζ and the actual w from the drifters is shown in Fig. 11. In grid cells where both values are available, the ratio $w/w_{\mathrm{\nabla }\mathit{\zeta }}$ is computed (panel a). It is almost entirely positive, reflecting that both estimates consistently show upwelling in this area. The few values greater than 1 indicate that a downwelling tendency from another source is softening the wind effect. Figure 11B compares the histograms of the full vertical velocity and the Ekman pumping component. It shows that w reaches values of 50–60 m d⁻¹, while w_∇ζ has a more peaked distribution, reaching only 20 m d⁻¹.

Figure 11 (A) Ratio of Ekman pumping to vertical velocity calculated from drifter divergence on 23 February. (B) Comparison of probability density functions for Ekman pumping estimates and for vertical velocity from drifters for 23 February.

For 25 February, the background transport direction (northwestward, Fig. 10B) and the cyclone vorticity structure (Fig. 7B) characterized by a maximum in the core of the eddy lead us to expect upwelling to occur north of the maximum and downwelling south of the maximum, consistent with Fig. 3. Because of the reduced drifter coverage on this day that limits the gradient calculation and because of the orientation of the vorticity maximum, the computation of w_∇ζ is limited to the downwelling side in the northwestern portion of the eddy and to the upwelling side in the southeastern portion (Fig. 10E), with maximum w_∇ζ≈10 m d⁻¹. For this day, the drifter-based w estimates have very limited coverage, but the surface divergence δ (Fig. 7D) exhibits a pattern of positive and negative δ to the northwest and to the southeast of the maximum ζ, respectively, suggestive of upwelling and downwelling as in Fig. 10E.

Finally, for 26 February, the w_∇ζ estimates (Fig. 10F) are mostly confined within the southernmost emerging eddy structure, since the elongated filament is too narrow to compute significant cross gradients. Due to the alternating negative and positive vorticity within the structure (Fig. 8B), and the transport direction (Fig. 10C), a pattern of positive and negative w_∇ζ to the east and west, respectively, emerges (Fig. 10F). For this day, no w estimates are available from drifters for comparison, but the δ estimates (Fig. 8D) show a consistent pattern of positive and negative divergence in the same area. We note that the main convergence/downwelling area observed in the central structure between the two eddies, and consistent with the isopycnal drop (Figs. 7C to 8C), does not have an obvious signature in terms of wind effects and likely results from internal dynamics, as suggested by the subsurface analysis of Middleton et al. (2025).

In addition to nonlinear Ekman pumping effects, the wind is also likely to induce NIOs that can affect the divergence and w distributions. As discussed in Sect. 2, we are not in a position to provide a quantitative estimate of these effects since we have only a partial picture of the eddy vorticity structure, which is moreover rapidly evolving, as is the wind. On the other hand, some general observations should be considered. For instance, in the presence of a vorticity maximum, we can expect that the NIOs tend to create two bands of opposite δ on either sides of the maximum (Esposito et al., 2023), with alternating sign in time at the inertial period (approximately 19 h). For 25 February, a maximum ζ is indeed recorded along the axis of the eddy (Fig. 7), with bands of opposite δ on each side. It is therefore possible that NIOs contribute to the observed pattern, even though we cannot verify whether or not the sign and the time dependence recorded by the drifters are indeed consistent with NIOs fluctuations. Similarly, NIOs could also contribute to the observed patterns in the submesoscale eddies on 26 and 28 February (Figs. 8, 9).

6 Summary, Discussion and Conclusions

The surface evolution of a splitting eddy under windy conditions has been investigated using drifter clusters that provide estimates of kinematic properties, hydrographic sections, and reanalysis wind data during the period 23–28 February 2022, aiming at several objectives.

The first main objective is to provide a fine scale description of the evolution of vorticity ζ, strain rate α and divergence δ in the mixed layer. The drifter sampling yields a high resolution (2 km binning) description of the eddy, even though the coverage is not always complete. Surface drifters with drogues centered in the first meter, CARTHE and CODE (CaC), together provide a detailed look initially at the flank of the original eddy, then at the elongated eddy, and finally at the constellation of smaller eddies after the splitting, while SVP drifters drogued at 15 m are more sparse but confirm that the flow is consistent within the surface layer. Results show high values of vorticity exceeding f, corresponding to a Rossby number of order 1 or higher, while divergence values also reach f, indicative of the submesoscale nature of the flow.

The patterns of vorticity and strain rate show the evolution of an eddy initially influenced by large-scale straining during 23 February, subsequently elongating and reaching maximum core values on 25 February, and then collapsing in the center on 26 February with the emergence of completely separated smaller eddies by 28 February. δ shows a more complex pattern, with bands of alternating sign, indicative of upwelling and downwelling. Positive divergence is prevalent in the southern flank of the eddy sampled during 23 February, while convergence is more prominent in the collapsing center on 26 February.

Vertical velocity w estimates are derived for 23 February, when surface drifters are co-located with SVP drifters. w values are found to be consistent with surface δ, showing mostly upwelling in the covered area, with values reaching 50–60 m d⁻¹. These values are consistent with previous findings in the Northwest Mediterranean Sea in the presence of frontal dynamics (Tintoré et al., 1991; Tarry et al., 2021).

Additional context for these results is obtained from Middleton et al. (2025), who analyzed a variational mapping of subsurface (below 20 m) hydrological data of the same eddy and an idealized model without external shear and wind stress. Overall, the surface evolution of ζ and α is found to be comparable to their subsurface expression, confirming that the main dynamics of the eddy splitting in the surface mixed layer are mostly controlled by internal processes. The δ pattern and evolution, on the other hand, are significantly different at the surface. While the subsurface divergence exhibits a smooth quadrupole pattern, indicative of cyclogeostrophic effects, here the presence of alternating bands suggests that other processes come into play. This is consistent with previous results that, while vorticity and strain are controlled by mesoscale and submesoscale processes, surface divergence is directly influenced by air-sea interaction and especially by wind forcing (Esposito et al., 2021).

The second main objective of the paper is to investigate the interaction between the eddy and the wind, with special emphasis on the resulting δ and vertical velocity w, derived from the vertical gradient of δ and the continuity equation. We estimate the nonlinear Ekman pumping velocity w_Eknl from wind data and ζ estimates, establishing the expected tendency of the surface flow, even though the high variability of the wind likely prevents its full realization. There are two components of w_Eknl, one that depends on vorticity w_ζ, resulting in downwelling within the cyclone, and one that depends on vorticity gradients w_∇ζ, resulting in a dipolar structure. The main contribution of Ekman pumping to the observed w pattern is found to be given by w_∇ζ, which causes alternating downwelling and upwelling on the two sides of the cyclone. The observed w and δ estimates reflect these expected w_∇ζ patterns. In addition, w_∇ζ is roughly of the same order of magnitude as the full w, indicating that wind effects on surface divergence and vertical velocity are significant. Overall, w was found to exhibit higher values (both negative and positive) than w_∇ζ, which is to be expected, both because other processes contribute to the vertical velocity and because the wind data are smoothed over much larger scales than the drifter observations, resulting in the smoothing of the w_∇ζ estimates. In addition to Ekman pumping, wind effects on divergence also include higher frequency processes such as NIOs. While a quantitative estimate of these effects is not possible, the observed vorticity and divergence patterns are consistent with an NIO contribution (Esposito et al., 2023).

The wind-induced fine-scale variability in vertical velocity found here is expected to impact the fine-scale pattern of biological tracers, interacting with the nutrient upwelling generally associated with cyclones and with the subduction of enriched waters during cyclone splitting (McGillicuddy et al., 2007; Mahadevan et al., 2008). Our findings indeed suggest that w_∇ζ could play an important role in the distribution of tracers such as chlorophyll and oxygen in the surface layer, even though its dipole structure tends to cancel upward and downward transport (Li et al., 2021), resulting in minimal net vertical exchange.

There are several other processes that could influence vertical transport, which were not analyzed due to insufficient available data. The Ekman effect has been found to lead to tilting in the vertical eddy axis (Stern, 1965; Tang et al., 2020), which could increase vertical velocity and transport (Li et al., 2022, 2023). Wind effects are also likely to modify turbulent processes in the surface layer, possibly resulting in the generation of submesoscale flows through the transient turbulent thermal wind balance (Dauhajre and McWilliams, 2018). Lastly, the combined effects of wind and frontal stratification could lead to submesocale instabilities that could result in subduction (Fox-Kemper et al., 2008; McWilliams, 2019).

In summary, the 2022 CALYPSO experiment yielded a rich observational dataset, specifically and fortuitously targeting a mesoscale eddy about to split into several submesoscale eddies, thereby providing an unusually detailed picture of the evolution. Yet, the resulting description provides only a partial 4D picture of the processes in action. A full accounting of the different processes, both steady and transient, however, remains out of reach. Describing all the facets and their interactions may well require further modeling work. Nevertheless, the present analysis definitively shows that cyclogeostrophic instabilities leading to eddy splitting are modified in the present case by strong shear and winds. Moreover, the vertical transport induced by the nonlinear interaction of the wind with the cyclone is a significant component of the total vertical exchange associated with the edges of the eddy. Other wind effects are likely to contribute to the high spatial variability of the divergence patterns observed by the drifters.